#OK to post homework
#Austin DeCicco, 4/18/26, Assignment 22

#BEGIN REFORMATTED C21
############################################################################################################################################################

with(combinat):

#SubsPi(f,x,pi): inputs a polynomial f(x[1], ..., x[n]) (n:=nops(pi)) and ouputs the new polynomial obtained by substituting x[i]->x[pi[i]]
SubsPi:=proc(f,x,pi): local n, i: n:=nops(pi): subs({seq(x[i]=x[pi[i]],i=1..n)},f); end proc:

#IsSymS(f,x,n): checks that all the images of f under the symmetric group coincide with f (that is a symmetric polynomial)
IsSymS:=proc(f,x,n): local Sn, pi: Sn:=permute(n): evalb({seq(SubsPi(f,x,pi), pi in Sn)}={f}); end proc:

#IsSym(f,x,n): checks that f(x[1], ..., x[n]) is symmetric the Pablo way
IsSym:=proc(f,x,n): local i: evalb({seq(subs({x[i]=x[i+1],x[i+1]=x[i]},f),i=1..n-1)}={f}); end proc:

#Symm(f,x,n): inputs an ARBITRARY polynomial f(x[1],...,x[n]) and outputs the sum of all the images of f under S_n
Symm:=proc(f,x,n): local Sn, pi: Sn:=permute(n): add(SubsPi(f,x,pi), pi in Sn); end proc:

#rSymm(f,x,n): inputs an ARBITRARY polynomial f(x[1],...,x[n]) and outputs the sum of all the distinct images
rSymm:=proc(f,x,n): local Sn, pi: Sn:=permute(n): convert({seq(SubsPi(f,x,pi), pi in Sn)}, `+`); end proc:

#Park(n,k): The set of partitions of n into exactly k parts
Park:=proc(n,k) option remember: local S, k1, S1, s1: if n<k then RETURN({}); end if; if k=0 then if n=0 then RETURN({[]}); else RETURN({}); end if; end if;
S:={}: for k1 from 0 to k do S1:=Park(n-k,k1): S:= S union {seq([op(s1+[1$k1]),1$(k-k1)], s1 in S1)}: end do; S; end proc:

#Par(n): The set of partitions of n, the clever way.
Par:=proc(n): local k: {seq(op(Park(n,k)), k=1..n)}; end proc:

#mL(x,n,L): inputs a variable name x a positive integer n and an integer partition L outputs the symmetrizer of x[1]^L[1]*...*x[nops(L)]^L[nops(L)], m_L(x[1], ..., x[n]) MONOMIAL symmetric polymomial
mL:=proc(x,n,L): local k, i: k:=nops(L): rSymm(mul(x[i]^L[i],i=1..k),x,n); end proc:

#ekn(k,n,x): The elementary symmetric polynomial of degree k in x[1], ..., x[n]
ekn:=proc(k,n,x): rSymm(mL(x,n,[1$k]),x,n); end proc:

#pkn(k,n,x): The power symmetric polynomial
pkn:=proc(k,n,x): rSymm(mL(x,n,[k]),x,n); end proc:

#eknC(k,n,x): A cleverer way to output e_k(x[1], ..., x[n])
eknC:=proc(k,n,x) option remember: if k>n then RETURN(0); end if; if n=0 then if k=0 then RETURN(1); else RETURN(0); end if; end if; expand(eknC(k,n-1,x)+ x[n]*eknC(k-1,n-1,x)); end proc:

############################################################################################################################################################
#END REFORMATTED C21

#BEGIN REFORMATTED C22
############################################################################################################################################################

with(linalg):

#mBase: The basis of Monomial Symmetric polynomial of total degree n
mBase:=proc(n,x): local P, i: P:=Par(n): [seq(mL(x,n,P[i]),i=1..nops(P))]; end proc:

#eBase(n,x): The e-base of the algebra of symmetric polynomials of degree n in x[1],...,x[n]
eBase:=proc(n,x): local P, i, j: P:=Par(n): expand([seq(mul(ekn(P[i][j],n,x),j=1..nops(P[i])),i=1..nops(P))]); end proc:

#pBase(n,x): The p-base of the algebra of symmetric polynomials of degree n in x[1],...,x[n]
pBase:=proc(n,x): local P, i, j: P:=Par(n): expand([seq(mul(pkn(P[i][j],n,x),j=1..nops(P[i])),i=1..nops(P))]); end proc:

#hkn(k,n,x): A cleverer way to output h_k(x[1], ..., x[n])
hkn:=proc(k,n,x) option remember: if k=0 then RETURN(1); end if; if n=0 then if k=0 then RETURN(1); else RETURN(0); end if; end if; expand(hkn(k,n-1,x)+ x[n]*hkn(k-1,n,x)); end proc:

#hBase(n,x): The e-base of the algebra of symmetric polynomials of degree n in x[1],...,x[n]
hBase:=proc(n,x): local P, i, j: P:=Par(n): expand([seq(mul(hkn(P[i][j],n,x),j=1..nops(P[i])),i=1..nops(P))]); end proc:

#Coe(f,x,L): inputs any polynomial f of x[1], ..., x[n], and a list of integers L output the coeff. of x[1]^L[1]*x[2]^L[2]*...
Coe:=proc(f,x,L): local c, i: c:=f: for i from 1 to nops(L) do c:=coeff(c,x[i],L[i]): end do; c; end proc:

#mVec(f,x,n): inputs a symmetric poly of x[1], ..., x[n] of (total) degree n outputs the coeffi vector (of length nops(Par(n)) when you express it in terms of m-bases
mVec:=proc(f,x,n): local P, i: P:=Par(n): [seq(Coe(f,x,P[i]),i=1..nops(P))]; end proc:

#eTOmMat(n): The nops(Par(n)) by nops(Par(n)) matrix whose [i,j] entry is the coefficient of  mBase(n)[j] in eBase(n)[i], i.e. when eBase(n)[i] is expressed as a linear combination of m_L).
eTOmMat:=proc(n): local EB, x, i: EB:=eBase(n,x): matrix([seq(mVec(EB[i],x,n),i=1..nops(EB))]); end proc:

############################################################################################################################################################
#END REFORMATTED C22

#1. By using the command "inverse" in the linalg package, and procedure eTOmMat(n), Write a procedure mTOeMat(n).
#   Using this express m_[6](x[1], ..., x[6]), alias pkn(6,6,x), as a polynomial in ekn(1,6,x),... ekn(6,6,x).

#Verfiying understanding of E2M transformation
EvalE2M:=proc(n):
local tr, S, i, j, m:
tr:=eTOmMat(n): m:=nops(Par(n)):
S:={seq(evalb(expand(eBase(n,x)[i])=expand(add(mBase(n,x)[j]*tr[i,j],j=1..m))),i=1..m)}:
if nops(S) > 1 then return(false); end if; S[1]; end proc:

mTOeMat:=proc(n): local Matrix: Matrix:=eTOmMat(n): inverse(Matrix); end proc:

#Verfiying validity of M2E transformation
EvalM2E:=proc(n):
local tr, S, i, j, m:
tr:=mTOeMat(n): m:=nops(Par(n)):
S:={seq(evalb(expand(mBase(n,x)[i])=expand(add(eBase(n,x)[j]*tr[i,j],j=1..m))),i=1..m)}:
if nops(S) > 1 then return(false); end if; S[1]; end proc:

#2. Write a procedure pBase(n,x), analogous to eBase(n,x) and hBase(n,x)

#This was provided in class, pBase:=proc(n,x): local P, i, j: P:=Par(n): expand([seq(mul(pkn(P[i][j],n,x),j=1..nops(P[i])),i=1..nops(P))]); end proc:

#3. Write procedures pTOmMat(n) and mTOpMat(n) that input a positive integer n and outputs the nops(Par(n)) by nops(Par(n)) transition matrices from the p basis to the m basis, and vice versa.

pTOmMat:=proc(n): local EB, x, i: EB:=pBase(n,x): matrix([seq(mVec(EB[i],x,n),i=1..nops(EB))]): end proc:
mTOpMat:=proc(n): local Matrix: Matrix:=pTOmMat(n): inverse(Matrix); end proc:

#Verfiying validity of P2M transformation
EvalP2M:=proc(n):
local tr, S, i, j, m:
tr:=pTOmMat(n): m:=nops(Par(n)):
S:={seq(evalb(expand(pBase(n,x)[i])=expand(add(mBase(n,x)[j]*tr[i,j],j=1..m))),i=1..m)}:
if nops(S) > 1 then return(false); end if; S[1]; end proc:

#Verfiying validity of M2P transformation
EvalM2P:=proc(n):
local tr, S, i, j, m:
tr:=mTOpMat(n): m:=nops(Par(n)):
S:={seq(evalb(expand(mBase(n,x)[i])=expand(add(pBase(n,x)[j]*tr[i,j],j=1..m))),i=1..m)}:
if nops(S) > 1 then return(false); end if; S[1]; end proc:

#Verfiying validity of all transformations
Eval22:=proc(n):
local S, i: S:={}:
for i from 1 to n do
S:=S union {EvalE2M(i),EvalM2P(i),EvalP2M(i),EvalM2E(i)}:
end do; if nops(S) > 1 then return(false); end if; S[1]; end proc:

if Eval22(4) then print("All transformation matricies coded correctly and verified."); end if;