# OK to post homework
# Aurora Hiveley, 4/1/26, Assignment 18

Help:=proc(): print(`diceRolls(n), normalMom(n), AppxAve22(n,K), AppxJike(n,x,K)`):   end:


### Problem 1

# Using Moms, find the exact avergae and s.d. for rolling a fair die n times, and using SMos prove that 
# its scaled moments up to the 20th agree with those of the standard normal distribution
# Conjecture an expression, in n, for the averge of the [2,2] entry in a standard Young tableau of shape [n,n,n]

# probabilistic gen func for n dice rolls
diceRolls := proc(n) local k:
expand((add(x^k, k=1..6))^n):
end:

# Moms(diceRolls(n), x, 2);
# output: [7*n/2, sqrt(105)*sqrt(n)/6]


# the nth moment of the standard normal distribution 
normalMom:= proc(n) 
# recall only the even moments are nonzero, so:
if n mod 2 = 1 then
  0:
else
  (2*n/2)!/(2^(n/2)*(n/2)!):
fi:
end:

## to compare the results
# S := SMoms(diceRolls(n), x, 20); 
## pick some sufficiently large n, say n = 100
# S1 := [seq(subs(n=100,S[i]),i=1..nops(S))]:

# T := [seq(normalMom(n), n=1..20)];
# for i from 1 to 20 do 
#   print(evalf(S1[i]),T[i],evalf(abs(S1[i]-T[i])));
# od: 

## output: 
#                         350., 0, 350.
#                  17.07825129, 1, 16.07825129
#                           0., 0, 0.
#                 2.987314286, 3, 0.01268571429
#                           0., 0, 0.
#                 14.81046045, 15, 0.1895395452
#                           0., 0, 0.
#                 102.3624285, 105, 2.637571514
#                           0., 0, 0.
#                 905.7572293, 945, 39.24277074
#                           0., 0, 0.
#                9754.081800, 10395, 640.9181998
#                           0., 0, 0.            
#              1.236127613 10^5 , 135135, 11522.23866
#                           0., 0, 0.
#           1.799858310 10^6 , 2027025, 2.271666898 10^5
#                           0., 0, 0.
#           2.957465892 10^7 , 34459425, 4.884766079 10^6
#                           0., 0, 0.
#          5.408205888 10^8 , 654729075, 1.139084862 10^8




### Problem 2

# AppxAve22(n,K) that takes K random standard Young tableaux of shape [n,n,n], finds its [2,2] entry and takes the average

AppxAve22 := proc(n,K) local i:
add(RandSYT([n,n,n])[2,2], i=1..K)/K:
end:

# What is
# AppxAve22(30,1000)?
# output:  1237/200 =  6.185000000

# AppxAve22(30,10000)?
# output:  31177/5000 = 6.235400000


# Are these close? Do they agree with your conjecture for the exact value?
# pretty close together! 




### Problem 3 

# AppxJike(n,x,K): approximates the prob. gen. function of the weight pi->nops(RS(pi)[1]) by taking K random permutations (using randperm(n))

AppxJike := proc(n,x,K) local i:
add(x^nops(RS(randperm(n))[1]), i=1..K):
end:


# PlotDist(AppxJike(100,x,1000))
# PlotDist(AppxJike(200,x,1000))

# Do they look alike?
# they do look very similar!





### copied from C18.txt 


Help18:=proc(): print(` OT(L), RandSYT(L), Moms(f,x,r), WtE(S,f) , Jike(n,x)`):
print(`Moms(f,x,r), SMoms(f,x,r) `):
end:

with(combinat):

Jike:=proc(n,x) local S,pi:
S:=permute(n) :
add(x^nops(RS(pi)[1]), pi in S):
end:


#PlotDist(f,x): plots the prob. distribution inspired by the weight enumerator f in the variable x after it is
#turned to a prob. distribution (after dividing by subs(x=1,f)
PlotDist:=proc(f,x) local f1,i:
f1:=f/subs(x=1,f):
plot([seq([i,coeff(f1,x,i)],i=ldegree(f1,x)..degree(f1,x))]):
end:

#WtE(S,f,x): the weight-enumerator of the set S according to the statistic f(s)
#WtE(permute(5), pi->nops(RS(pi)[1]),x  );
WtE:=proc(S,f,x) local s:
add(x^f(s), s in S):
end:

#Moms(f,x,r): 
Moms:=proc(f,x,r) local mu,L,f1,sig,i:

f1:=f/subs(x=1,f):

mu:=subs(x=1,diff(f1,x)):

L:=[mu]:

f1:=f1/x^mu:
#g.f. of X-mu
#Centralized pgf

f1:=x*diff(f1,x):

f1:=x*diff(f1,x):
sig:=sqrt(subs(x=1,f1)):
L:=[mu,sig]:

for i from 3 to r do
f1:=x*diff(f1,x):
L:=[op(L), subs(x=1,f1)]:
od:
L:
end:

SMoms:=proc(f,x,r) local L,sig,i:
L:=Moms(f,x,r):
sig:=L[2]:
[op(1..2,L),seq(L[i]/sig^i    , i=3..r)]:
end:


OT:=proc(L) local n,c:
n:=convert(L,`+`):
c:=Cells(L)[rand(1..n)()]:
while HLc(L,c)>1 do
 c:=Hook(L,c)[rand(1..HLc(L,c))()]:
od:
c:
end:

#RandSYT(L): a random SYT of shape L
RandSYT:=proc(L) local k,c,i,L1,n,Y1:

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

n:=convert(L,`+`):

k:=nops(L):
c:=OT(L):
i:=c[1]:

L1:=[op(1..i-1,L),L[i]-1,op(i+1..k,L)]:
if L1[-1]=0 then
L1:=[op(1..nops(L1)-1,L1)]:
fi:

Y1:=RandSYT(L1):
if nops(L1)=k then
[op(1..i-1,Y1),[op(Y1[i]),n],op(i+1..k,Y1)]:
else
[op(1..k-1,Y1),[n]]:
fi:
end:


#old stuff
#C17.txt, March 26, 2006
Help17:=proc():  print(`Cells(L), Conj(L), NuSYT(L)`): end:

#Cells(L): the set of n (=sum(L)) cells [i,j]in the shape L
Cells:=proc(L) local k,i,j:
k:=nops(L):
{seq(seq([i,j]   , j=1..L[i]  ), i=1..k)}
end:

Conj:=proc(L) local k,C1,i1,L1,i:
option remember:
if L=[] then
RETURN([]):
fi:
k:=nops(L):
L1:=[seq(L[i]-1,i=1..k)]:

for i1 from 1 to nops(L1) while L1[i1]>0 do od:
i1:=i1-1:

L1:=[op(1..i1,L1)]:

C1:=Conj(L1):
[k,op(C1)]:
end:

#Hook(L,c): The set of the cells in the hoo corresponding to the cell c=[i,j]
#(i.e. the set of cells to the right and to the bottom of c
Hook:=proc(L,c) local k,i,j,C,i1,j1,i2:
k:=nops(L):
i:=c[1]:
j:=c[2]:
if not (i>=1 and i<=k) then
  RETURN(FAIL):
fi:

if not(j>=1 and j<=L[i]) then
  RETURN(FAIL):
fi:

C:={seq([i,j1],j1=j..L[i])}:

for i2 from i to k while j<=L[i2] do od:
i2:=i2-1:
C:=C union {seq([i1,j],i1=i..i2)}:
end:

#HL(L,c): the hook-length of cell c in the shape L
HL:=proc(L,c): nops(Hook(L,c)):end:

HLc:=proc(L,c) local i,j,L1:
L1:=Conj(L):
i:=c[1]:
j:=c[2]:
L[i]-j+L1[j]-i+1:

end:

#NuSYTc(L): implementing the Frame-Robinson-Thrall Hook Length Formula
NuSYTc:=proc(L) local n,C,i,c:
n:=add(L[i],i=1..nops(L)):
C:=Cells(L):
n!/mul(HL(L,c),c in C):
end:


#NuSYTcc(L): implementing the Frame-Robinson-Thrall Hook Length Formula Cleverly
NuSYTcc:=proc(L) local n,C,i,c:
n:=add(L[i],i=1..nops(L)):
C:=Cells(L):
n!/mul(HLc(L,c),c in C):
end:









#old stuff

#C16.txt, March 23, 2026
Help16:=proc(): print(`RSleft(pi), RS1(Y,i), RS(pi) `): end:

#RS1(Y,i): inputs a partial Young tableau and another integer i NOT yet in Y
#places it in the right place, by a bumping process it returns a tableau
#with one more box followed by the name of the row where it settled
RS1:=proc(Y,i) local k,NewY,lucy,bumpee,i1,j:
if Y=[] then
 RETURN([[i]],1):
fi:
k:=nops(Y):
lucy:=RS11(Y[1],i):
NewY[1]:=lucy[1]:
bumpee:=lucy[2]:

if bumpee=0 then
  RETURN([NewY[1],op(2..k,Y)],1):
fi:

for i1 from 2 to k  while bumpee<>0 do
 lucy:=RS11(Y[i1],bumpee):
  NewY[i1]:=lucy[1]:
 bumpee:=lucy[2]:
if bumpee=0 then
 RETURN([seq(NewY[j],j=1..i1),op(i1+1..k,Y)],i1):
fi:
od:

[seq(NewY[j],j=1..k),[bumpee]],k+1:
end:


with(combinat):

#RS(pi): inputs a permutation pi of ({1, ..., n:=nops(pi)) and outputs a pair of SYT of the SAME shape (with n boxes)
#The Robinson-Schenstead algorithm
RS:=proc(pi) local Yl, i, Yr,eaea, p:
if pi=[] then
 RETURN([[],[]]):
fi:
Yl:=[[pi[1]]]:
Yr:=[[1]]:
for i from 2 to nops(pi) do
eaea:=RS1(Yl,pi[i]):
Yl:=eaea[1]:
p:=eaea[2]:

if p<=nops(Yr) then
 Yr:=[op(1..p-1,Yr),[op(Yr[p]),i],op(p+1..nops(Yr),Yr)]:
else
Yr:=[op(Yr),[i]]:
fi:
od:

[Yl, Yr]:
end:






#RSeft(pi): inputs a permutation pi (of size nops(pi)) and outputs
#of whatever shape (with n boxes)
RSleft:=proc(pi) local Y,i:
Y:=[]:
for i from 1 to nops(pi) do
Y:=RS1(Y,pi[i])[1]:
od:
Y:
end:

#old stuff
#C14.txt; March 9, 2026

Help14:=proc(): print(` NuSYT(L), SYTpairs(n) , NuSYTpairs(n), RS11(a,i) `): end:

#RS11(a,i): inputs an INCREASING list of positive integers, a, and another positive integer i
#outputs a pair a1,j, where (usually a1 is of the same length as a) and i is put where it
#belongs and j is the entry that it bumped, unless i is larger than all the members of a
#(i.e. larger than a[-1]) then a1 is [op(a),i], and j is 0
RS11:=proc(a,i) local k,j:
k:=nops(a):
for j from 1 to k while a[j]<i do od:
if j=k+1 then
RETURN([op(a),i],0):
fi:
[op(1..j-1,a),i,op(j+1..k,a)],a[j]:
end:

NuSYTpairs:=proc(n) local S,P, p:
option remember:
P:=Par(n):
S:=0:
for p in P do
S:=S+NuSYT(p)^2:
od:
S:
end:

SYTpairs:=proc(n) local S,P, p,S1,s1,s2:
option remember:
P:=Par(n):
S:={}:
for p in P do
S1:=SYT(p):
S:= S union {seq(seq([s1,s2], s1 in S1), s2 in S1)}:
od:
S:
end:

#NuSYT(L): The NUMBER of Standard Young tableaux of shape L
NuSYT:=proc(L) local i, k,S,L1:
option remember:
k:=nops(L):

if k=0 then
  RETURN(1):
fi:
S:=0:
#Look for all legal rows where Mr. or Ms. n can place themselves
for i from 1 to k-1 do
if L[i]>L[i+1] then
 L1:=[op(1..i-1,L),L[i]-1,op(i+1..k,L)]:
 S:=S+NuSYT(L1): 
fi:
od:

if L[k]>1 then
L1:=[op(1..k-1,L),L[k]-1]:

S:=S+ NuSYT(L1):
else
L1:=[op(1..k-1,L)]:
S:=S+NuSYT(L1):
fi:
S:
end:



#old stuff
#C13.txt
Help13:=proc(): print(` PFG(L), SYT(L), PSYT(n) `): end:
Help12:=proc():
 print(`Park(n,k), Par(n), ParN(n,k)`):
end:
ParN:=proc(n,k) local s,S,T: S:=Par(n):T:={}:
for s in S do
 if s[1]=k then
  T:=T union {s}:
 fi:
od:
T:
end:
#Park(n,k): The set of partitions of n into exactly k parts
Park:=proc(n,k) local S,k1,S1,s1:
option remember:

if n<k then
   RETURN({}):
fi:

if k=0 then
 if n=0 then
   RETURN({[]}):
   else
    RETURN({}):
fi:
fi:

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)}:
od:
S:
end:


#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:

#PFG(L): The Ferrers graph of the partition L
PFG:=proc(L) local i:
for i from 1 to nops(L) do
lprint(1$L[i]):
od:
print(`----------------`):
end:




ParN:=proc(n,k) local s,S,T: S:=Par(n):T:={}:
for s in S do
 if s[1]=k then
  T:=T union {s}:
 fi:
od:
T:
end:

#Def. Standard Young tableau of shape L L=[3,3,2] a partition of 8
#
#111111
#111111
#11
#
#A way of filling {1, ..., 8} the empty boxes such that every row and every column are increasing
#
#L=[2,2]
#
#12
#34
#
#13
#24    Number of of YT of shape [2,2]
#
#Where can n (the tallest person) seat w/o disturbing shorter people
#

#SYT(L): The SET of Standard Young tableaux of shape L
SYT:=proc(L) local i, k,n,S,L1,S1,s1:
option remember:
k:=nops(L):
n:=add(L[i],i=1..k):
if k=0 then
  RETURN({[]}):
fi:
S:={}:
#Look for all legal rows where Mr. or Ms. n can place themselves
for i from 1 to k-1 do
if L[i]>L[i+1] then
 L1:=[op(1..i-1,L),L[i]-1,op(i+1..k,L)]:
 S1:=SYT(L1):

S:=S union {seq(
[op(1..i-1,s1),[op(s1[i]),n],op(i+1..k,s1)]
,s1  in S1)}:
fi:
od:

if L[k]>1 then
L1:=[op(1..k-1,L),L[k]-1]:

S1:=SYT(L1):
S:=S union {seq(
[op(1..k-1,s1),[op(s1[k]),n]]
,s1  in S1)}:
else
L1:=[op(1..k-1,L)]:
S1:=SYT(L1):
S:=S union {seq(
[op(1..k-1,s1), [n]]
,s1  in S1)}:
fi:
S:
end:

#PSYT(Y): prints the SYT Y
PSYT:=proc(Y) local i:
for i from 1 to nops(Y) do
lprint(op(Y[i])):
od:
end: