#OK to post
#Yuxuan Yang, April 30th, Assignment 26

with(combinat):
with(linalg):

#1
#I modified hkn ekn with negative subscripts are defined to be zero 
ScJT1:=proc(L,x,m) local i,j,n:
n:=nops(L):
det(matrix([seq([seq(hkn(L[i]+j-i,m,x),j=1..n)],i=1..n)])):
end:

ScJT2:=proc(L,x,m) local i,j,n,Lc:
Lc:=Conj(L):
n:=nops(Lc):
det(matrix([seq([seq(ekn(Lc[i]+j-i,m,x),j=1..n)],i=1..n)])):
end:
#evalb(ScJT1([3, 2, 1], x, 7) = ScY([3, 2, 1], x, 7))
#evalb(ScJT2([3, 2, 1], x, 7) = ScY([3, 2, 1], x, 7))


#Conj(L):Conjugate partition
Conj:=proc(L) local L1,k,i:

if L=[] then
  RETURN([]):
fi:

k:=nops(L):

L1:=L-[1$k]:

for i from 1 to nops(L1) while L1[i]<>0 do  od:
L1:=[op(1..i-1,L1)]:
[k, op(Conj(L1))]:

end:




#2
StartShape:=proc(Y,k) local i,j,l,S:
S:=[]:
for i from 1 to nops(Y) do
j:=0:
for l from 1 to nops(Y[i]) do
if Y[i][l]<=k then 
j:=j+1:
fi:
od:

if j>0 then
S:=[op(S),j]:
fi:
od:
S:
end:

EstStart:=proc(L,k,K) local i,pnk,fre,YT,S,j:
pnk:=Pn(k):
fre:=[0$nops(pnk)]:
for i from 1 to K do 
S:=StartShape(GNW(L),k):
for j from 1 while pnk[j]<>S do od:
fre[j]+=1:
od:
fre:
end:	
#EstStart([9$12],5,10000)
#[48, 1063, 1992, 3005, 2325, 1471, 96]












#Old
#Yuxuan Yang, April 24th, Assignment 25

with(combinat):

#1
HookSet:=proc(L,cell) local i,j,S:
S:={}:
for i from cell[1] to nops(L) while L[i]>=cell[2] do
 S:=S union {[i,cell[2]]}:
od:
for j from cell[2] to L[cell[1]] do
 S:=S union {[cell[1],j]}:
od:

S:
end:
 
#2
OneStepSNW:=proc(L) local m,n,i,j,H,k,allpos,pos:
allpos:=[seq(seq([n,m],m=1..L[n]),n=1..nops(L))]:
pos:=allpos[rand(1..nops(allpos))()]:
i:=pos[1]:
j:=pos[2]:
for k from 1 to add(L) do
H:=HookSet(L,[i,j]):
if nops(H)=1 then 
 RETURN([i,j]):
fi:
pos:=H[rand(2..nops(H))()]:
i:=pos[1]:
j:=pos[2]:
od:
pos:
end:

#3

GNW:=proc(L) local i,j,pos,YT,n,L1:
n:=add(L):
if nops(L)=1 then
RETURN([[seq(i,i=1..L[1])]]):
fi:
pos:=OneStepSNW(L):
L1:=L:
if L1[pos[1]]>1 then
L1[pos[1]]:=L1[pos[1]]-1:
YT:=GNW(L1):
YT[pos[1]]:=[op(YT[pos[1]]),n]:
RETURN(YT):
else
L1:=[op(1..nops(L1)-1,L1)]:
YT:=GNW(L1):
YT:=[op(YT),[n]]:
RETURN(YT):
fi:
end:
Pnk:=proc(n,k) local k1,L,L1:
option remember:

if not (type(n,integer) and type(k,integer) and n>=1 and k>=1 )then
 RETURN([]):
fi:
 
if k>n then
 RETURN([]):
fi:

if k=n then
 RETURN([[n] ]):
fi:

L:=[]:

for k1 from min(n-k,k) by -1 to 1 do
 L1:=Pnk(n-k,k1):
 L:=[op(L), seq([k, op(L1[j])],j=1..nops(L1))]:
od:

L:

end:


#Pn(n): The list of integer partitions of n in rev. lex. order
Pn:=proc(n) local k:option remember:[seq(op(Pnk(n,n-k+1)),k=1..n)]:end:



#Elementary symmetic polynomial

ekn:=proc(k,n,x) option remember:
if k=0 then
  RETURN(1):
fi:
if k>n or k<0 then
 RETURN(0):
fi:
expand(ekn(k,n-1,x)+x[n]*ekn(k-1,n-1,x)):
end:


#Full homog. symmetic polynomial
hkn:=proc(k,n,x) option remember:

if k=0 then
  RETURN(1):
fi:
if n=0 or k<0 then
  RETURN(0):
fi:
expand(hkn(k,n-1,x)+x[n]*hkn(k-1,n,x)):
end:




#5
delta:=proc(x,n) local j,i:
product(product(x[i]-x[j],j=i+1..n),i=1..n-1):
end:

guesssc:=proc(L,x,m) 
expand(Sc(L,x,m)*delta(x,m)):
end:

detpart:=proc(L,x,m) local i,j,M,degree,size,xnew:
size:=max(m,nops(L)):
degree:=[0$size]:
xnew:=[1$size]:
for i from 1 to m do
degree[i]:=degree[i]+m-i:
xnew[i]:=x[i]:
od:
for i from 1 to nops(L) do
degree[i]:=degree[i]+L[i]:
od:
M:=[seq([seq(xnew[i]^degree[j],j=1..size)],i=1..size)]:
expand(det(matrix(M))):
end:

checksc:=proc(L,x,m)
evalb(guesssc(L,x,m)=detpart(L,x,m)):
end:

#I checked several examples. The formula should work.
#The formula is in fact a generalized vandermonde det over the regular vandermonde det.
#6 no idea how to prove.



#ScY(L,x,m): The Schur polynomial in x[1],..., x[m]
#of the shape L, according to the determinant formula
#rediscovered by Yuxuan
ScY:=proc(L,x,m): normal(detpart(L,x,m)/delta(x,m)):end: