# OK to post homework
# Lucy Martinez, 03-31-2026, Assignment 18

with(combinat):


# Question 1:
# Using Moms, find the exact average 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
# ANSWER: [Recall that the first 2 terms of Moms returns [mu,sig] ]
# [First part] 
# Running Moms(((x+x^2+x^3+x^4+x^5+x^6))^n,x,2);
#  returns the following information:
#  The average is (7*n)/2 and the s.d. is 1.707825129*sqrt(n)
# [Second part]
# For this part, I did the following:
#  fdie:=x+x^2+x^3+x^4+x^5+x^6:
#  n:=100:
#  fn:=expand(fdie^n):
#  Sdie:=SMoms(fn, x, 20):
# The following loop compares the difference between the scaled moments for
#  rolling a fair die and the scaled moments of the standard normal distribution
# for i from 3 to 20 do
#  die:=evalf(Sdie[i]):
#  Ndis:= (NormalSMom(i)):
#  dif:=abs(die - Ndis): #this is the difference of the i-th scaled moments for both cases
#  print(die, Ndis,dif):
# od:
# Here are the numbers I got.
# Of course, we would need to increase the number of times we roll the die
#  to obtain more accurate results.
# i=3: 0., 0, 0.
# i=4: 2.987314277, 3, 0.012685723
# i=5: 0., 0, 0.
# i=6: 14.81046039, 15, 0.18953961
# i=7: 0., 0, 0.
# i=8: 102.3624279, 105, 2.6375721
# i=9: 0., 0, 0.
# i=10: 905.7572219, 945, 39.2427781
# i=11: 0., 0, 0.
# i=12: 9754.081708, 10395, 640.918292
# i=13: 0., 0, 0.
# i=14: 123612.7600, 135135, 11522.2400
# i=15: 0., 0, 0.
# i=16: 1.799858287*10^6, 2027025, 227166.713
# i=17: 0., 0, 0.
# i=18: 2.957465851*10^7, 34459425, 4.88476649*10^6
# i=19: 0., 0, 0.
# i=20: 5.408205803*10^8, 654729075, 1.139084947*10^8


# Scaled moments of standard Normal N(0,1)
# Odd  moments = 0
# Even moment of order 2k = (2k-1)!! = 1*3*5*...*(2k-1)
NormalSMom:=proc(r)
 if r mod 2=1 then
  return(0):
 else
  return((2*r/2)!/(2^(r/2)*(r/2)!)):
 ##doublefactorial(r-1)   # (r-1)!! 
 fi:
end:


# Conjecture an expression, in n, for the average of the [2,2] entry 
#  in a standard Young tableau of shape [n,n,n]
STYPGF22:=proc(n,x) local L, f, T, allT:
L:=[n,n,n]:
allT:=SYT(L):
f:=add(x^(T[2][2]), T in allT):
f:
end:



# Question 2:
# Write a procedure AppxAve22(n,K) that takes K random standard Young tableaux
#  of shape [n,n,n], finds its [2,2] entry and takes the average
# What is AppxAve22(30,1000)?
#  ANSWER: 6221/1000 which is roughly 6.221
# AppxAve22(30,10000)?
#  ANSWER: 31177/5000 which is roughly 6.2354
# Are these close? Do they agree with your conjecture for the exact value?
#  ANSWER: They are close (off on the second, third and fourth digits).
AppxAve22:=proc(n,K) local i,Y,s22,s,L:
s:=0:
for i from 1 to K do
 L:=[n,n,n]:
 Y:=RandSYT(L):
 s22:=Y[2][2]:
 s:=s+s22:
od:
s/K:
end:



# Question 3: Write a procedure AppxJike(n,x,K)
#  that approximates the prob. gen. function of the weight pi->nops(RS(pi)[1])
#  by taking K random permutations (using randperm(n))
# Do PlotDist(AppxJike(100,x,1000),x)
# and PlotDist(AppxJike(200,x,1000),x)
# Do they look alike?
# ANSWER: Yes, they roughly look alike.
AppxJike:=proc(n,x,K) local i,s,pi:
s:=0:
for i from 1 to K do
 pi:=randperm(n):
 s:=s+x^nops(RS(pi)[1]):
od:
s/K:
end:


#########################From previous class:
# C18.txt, March 30, 2026

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

#Notes:
#Want to see how many rows does the Tableaux from the RSK give young
# For example if pi:=[1,2,4,3] then RS(pi)[2]=[[1,2,5],[3,4]]
#
#Turns out the number of rows of RS(pi)[1] is equal to the number
# of the largest decreasing subsequence

#Statistics Notes:
# If f=Sum(Pr(X=i), all i) is the probability generating function then
#  the average is defined as
#  Sum(Pr(X=i)*i, all i) 
# The r-th moment:
#  m_r(X)=add((X(s)-av)^r, s in S)/nops(S)
# The scaled moments called alpha is defined as
#  alpha_r=m_r/((s.d)^r) where s.d=standard deviation
#
# alpha_4=m_4/sig^4 is called the Kurtosis
# THE STANDARD DISTRIBUTION:
# m_r=0 if r is odd AND
# m_(2r)= 1*3*5*7*...*(2*r-1)=(2*r)!/(2^r*r!) (the even ones)


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

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

#f1 derived twice:
f1:=x*diff(f1,x):
f1:=x*diff(f1,x):
sig:=evalf(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(f,x,r): The scaled moments
SMoms:=proc(f,x,r) local L,sig,i,mu:
L:=Moms(f,x,r):
mu:=L[1]:
sig:=L[2]:
#L is the list of moments
# L[i]/sig^i are the scaled moments
[mu,sig,seq(L[i]/sig^i, i=3..r)]:
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:


Jike:=proc(n,x) local S,pi:
S:=permute(n):
add(x^nops(RS(pi)[1]), pi in S):
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:



##################################from before
# C17.txt, March 23, 2026

Help17:=proc():
 print(`Cells(L), Conj(L), Hook(L,c), HL(L,c), NuSYTc(L)`):
 print(`HLc(L,c), NuSYTcc(L)`):
end:

#S:=SYT([3,3,3]):
# nops(S)=42
#One way that we could get a random Young Tableaux is to call in Maple: 
# S[rand(1..nops(S))()];
#However, SYT(L) is constructing all Young Tableaux - not efficient to do this
#We will set up for next week so we can construct a way to generate random Young Tableaux
#
#Robinson-Frame-Thrall
#Suppose we have placed the following boxes:
# ***
# ***
# ***
#NuSYT(L)=n!/mul(nops(Hook(c)), c a cell in the shape L)
# Let L be of shape:
# ***
# ***
# ***
# Here, n=9 (9 cells)
# To count the hook of any cell: count the the number of 
#  stars in the same column and to the right of it (row)
# So now, the hook length for cell=(1,1) is 5 
#  the hook length for cell (1,2) is 4
#  the hook length for cell (2,2) is 3
# If we do it all, we get
#  9!/(5*4*3*4*3*2*3*2*1)=42


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

#Def: Conjugate of a partition: each column (left to right) becomes a row
# from top to bottom
# Example: 
# Say we have the following diagram
# 331
# ***
# ***
# *
# Then the conjugate is 322
# ***
# **
# **

#Conj(L): Given a partition L outputs the conjugate of L
Conj:=proc(L) local k,L1,i1,C1,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 hook corresponding
# to the cell c=[i,j]
# i.e. the set of cells to the right and to the bottom of call
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)}:
C:
end:


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

#HLc(L,c): the hook-length of the cell c in the shape L
# done cleverly
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 lenght 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 lenght formula
# done 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:


#######################################
Help16:=proc():
 print(`RS1(Y,i), RS1left(pi), RS(pi)`):
 print(` `):
end:

#We will finish the RSK algorithm
# Example:
# Given the following partial STY 
# 4 5 7
# 6 8
# Want to insert 2:
# 2 5 7  and bumpee=4 so now we have to insert 4 next
# 6 8    
# Want to insert 4 (after second row):
# 2 5 7
# 4 8    and bumpee=6 so now we have to insert 6 next
# Want to insert 6 (after third row):
# 2 5 7
# 4 8
# 6
#The previous SYT is the output  

#Now given a permutation pi=31452
#Begin STY:
# 3
# insert 1 and bump 3
# 1
# 3
# insert 4
# 1
# 3 4
# insert 5
# 1
# 3 4 5
# insert 2
# 1 2
# 3 4 5 


#RS1(Y,i): Inputs a partial Young tableaux and another integer i NOT yet in Y
# and places it in the right place, by a bumping process
# it returns a tableaux with one more box followed by the row 
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:



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



#RS(pi): inputs a permutation pi of ({1,2,...,n:=nops(pi)) and outputs
# a pair of SYT of the SAME shape (with n boxes)
#The Robinson–Schensted algorithm
RS:=proc(pi) local Y1, i,Yr,eaea,p:
if pi=[] then
 return([[],[]]):
fi:
Y1:=[[pi[1]]]:
Yr:=[[1]]:

for i from 2 to nops(pi) do
 eaea:=RS1(Y1,pi[i]):
 Y1:=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:
[Y1,Yr]:
#Yr is usually what they call Q
end:



############################
# C14.txt, March 09, 2026

Help14:=proc():
 print(`NuSYT(L), SYTpairs(n), NuSYTpairs(n)`):
 print(`RS11(a,i)`):
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:

#now we look for all the legal rows where the element n can be placed
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:


#Pairs of SYT of the same shape
# We are interested in pairs (Y1,Y2) such that
# shape of Y1 equals shape of Y2
# For n=2:
#  SYT([1,1])={[ [1], [2] ] } so there is only 1 of these STY of shape [1,1]
#
#  SYTpairs(2)= { [[[1, 2]], [[1, 2]]], [ [[1], [2]], [[1], [2]] ]}

#SYTpairs(n):
SYTpairs:=proc(n) local S,P,p,S1,s1,s2:
option remember:
P:=Par(n):
S:={}:
for p in P do
 S1:=SYT(p): #the set of Young Tableaux of the partition p (shape p)
 S:=S union {seq(seq([s1,s2], s1 in S1),s2 in S1)}:
od:
S:
end:

#NuSYTpairs(n):
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:



#RSK: Find a mapping from the set of permutations to the set of STYpairs(n)
#Input: permutation of {1,2,...,n}
#Output: a pair of SYTs of the SAME shape (with n boxes)
#Example: Given [1,4,5,9] then where can we place 10 so that
# it is still increasing? --> put it at the end
# What if I want to place the number 3 using RSK such that if I place it
#  somewhere other than at the end, we must bump the number that was in it
#  So, [1,4,5,9] becomes [1,3,5,9]
# 

#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 placed 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
# For example: RS1([11,12,13],4)= [4,12,13], 11
#  where the second output is the element that got bumped out
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:



#Now that we have the helper function RS11(a,i), we do the following:
#
#Input a Young Tableaux (YT) with distinct entries and an integer i
# NOT yet in the YT, outputs a larger YT with ONE more box, obtained as follows:
#  Try to insert i in the first row if it fits w/o bumping it
#  otherwise put it in its rightful place and bump somebody
#  now the one that got bumped has to bump someone else
#  and place the bumpee in the second row


#Example of the algorithm above
#Y= [ [1,4,8], [2,7], [3] ]
#i=5

#148 bumpee is 8
#27
#3

#becomes 

#145
#278
#3

#######################################
# C13.txt, March 05, 2026

Help13:=proc():
 print(`ParN(n,k), PFG(L), SYT(L), PSYT(Y)`):
end:


#ParN(n,k): outputs the partitions of n for which the first part equal to k
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:


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



#SYT(L): the set of standard young tableaux of shape L
SYT:=proc(L) local i,s1,k,n,S,L1,S1:
option remember:
k:=nops(L):
n:=add(L[i],i=1..k):

if k=0 then
 return({[]}):
fi:

S:={}:

#now we look for all the legal rows where the element n can be placed
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:



#Def: We are going to construct the Standard Young tableaux
# of shape L=[3,3,2], which is a partition of 8

# The following is the Ferrers graph for [3,3,2]
# 111
# 111
# 11

# Now, imagine each number "1" is an empty box. You fill each box
#  with the numbers {1,...,8} such that every row and every column
#  are increasing from left to right and top to bottom, respectively
# One way is the following:
# 123
# 456
# 78

# Example 2:
# Let L=[2,2] so n=4
# This is the shape:
# 11
# 11
#
# So the Young Tableaux are the following 2:
#
# 12
# 34
#
# AND
#
# 13
# 24

# Where can n (the tallest person) sit w/o disturbing shorter people?
#

# Below is a recap from last week:
#
# Recall that the number of partitions of n, denoted by p(n) satisfies
# sum(p(n)*q^n,n=0..infinity)=Prod(1/(1-q^k),k=1..infinity)
#
# G.F. of distinct partitions Prod((1+q^k),k=1..infinity)
# sum(d(n)*q^n,n=0..infinity)=(1+q)*(1+q^2)*(1+q^3)*...
#  where d(n) is the number of partitions with distinct parts
#
# sum(o(n)*q^n,n=0..infinity)=Prod(1/(1-q^(2*k-1)),k=1..infinity)=
#  1/((1-q)*(1-q^3)*(1-q^5)*...)
#  where o(n) is the number of partitions where all the parts are odd 
#
# Recall that (1+z)=(1-z^2)/(1-z)
# SO,
# 1/((1-q)*(1-q^3)*(1-q^5)*...) = (1-q^2)*(1-q^4)*(1-q^6)*(1-q^8)*.../((1-q)*(1-q^2)*(1-q^3)*(1-q^4)...)
#                           = ( (1-q^2)/(1-q) )*( (1-q^4)/(1-q^2) )*( (1-q^6)/(1-q^3) )*...
#                           =  (1+q)*(1+q^2)*(1+q^3)*...
# On the left hand side we have the G.F. of the partitions with odd parts
# On the right hand side we have the G.F. of the partitions with distinct parts
# For more detail, see: https://sites.math.rutgers.edu/~zeilberg/numtheory/L22.pdf




###################################
# C12.txt, March 02, 2026

Help12:=proc():
 print(`Park(n,k), Par(n)`):
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: