#OK to post homework
#Austin DeCicco, 4/3/26, Assignment 17

#BEGIN REFORMATTED C12
############################################################################################################################################################

with(combinat):

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

#A41s(N): OEIS sequence A41 done stupidly (construct the actual set) and take the cardinality
A41s:=proc(N): local n: [seq(nops(Par(n)),n=1..N)]; end proc:

#A41ok(N): OEIS sequence A41 done much more efficiently using Euler's generating function Sum(p(n)*q^n,n=0..infinity)=Prod(1/(1-q^k), k=1..infinity)
A41ok:=proc(N): local k, f, q: f:=mul(1/(1-q^k),k=1..N): f:=taylor(f,q=0,N+4): [seq(coeff(f,q,k),k=1..N)]; end proc:

#EulerM(N,q): The truncation to 1^N of Prod(1-q^k, k=1..infinity)
EulerM:=proc(N,q): local f, k: f:=mul(1-q^k,k=1..N): f:=taylor(f,q=0,N+1); end proc:

#ConjP(N): The conjectured powers of the polynomial in EulerM
ConjP:=proc(N): local j, L: L:=[]: for j from 1 to N while (3*j^2-j)/2<=N do L:=[op(L),(3*j^2-j)/2]: end do; for j from 1 to N while (3*j^2+j)/2<=N do
L:=[op(L),(3*j^2+j)/2]: end do; sort(L); end proc:

#EPTtest(N): testing Euler's pentagonal numbers (3*j^2+j)/2
EPTtest:=proc(N): local q, j, f, g, k: f:=mul(1-q^k,k=1..N):
g:=add((-1)^j*q^((3*j^2+j)/2), j=1..trunc(2*sqrt(N/3)+1))+add((-1)^j*q^((3*j^2-j)/2), j=1..trunc(2*sqrt(N/3)+1)): taylor(f-g,q=0,N+1); end proc:

############################################################################################################################################################
#END REFORMATTED C12

#BEGIN REFORMATTED C13
############################################################################################################################################################

#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]): end do; print(`----------------`): end proc:

#SYT(L): The SET of Standard Young tableaux of shape L
SYT:=proc(L) option remember: local i, k, n, S, L1, S1, s1: k:=nops(L): n:=add(L[i],i=1..k): if k=0 then RETURN({[]}); end if; S:={}:
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)}; end if; end do;
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)}; end if; S; end proc:

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

############################################################################################################################################################
#END REFORMATTED C13

#BEGIN REFORMATTED C14
############################################################################################################################################################

#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 end do; if j=k+1 then RETURN([op(a),i],0); end if; [op(1..j-1,a),i,op(j+1..k,a)],a[j]; end proc:

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

SYTpairs:=proc(n) option remember: local S, P, p, S1, s1, s2: 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)}: end do; S; end proc:

#NuSYT(L): The NUMBER of Standard Young tableaux of shape L
NuSYT:=proc(L) option remember: local i, k, S, L1: k:=nops(L): if k=0 then RETURN(1); end if; S:=0: 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); end if; end do; 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); end if; S; end proc:

############################################################################################################################################################
#END REFORMATTED C14

#BEGIN REFORMATTED C16
############################################################################################################################################################

#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); end if; 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);
end if; 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); end if; end do; [seq(NewY[j],j=1..k),[bumpee]],k+1; end proc:

#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([[],[]]); end if; 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]]; end if; end do; [Yl, Yr]; end proc:

#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]: end do; Y; end proc:

############################################################################################################################################################
#END REFORMATTED C16

#BEGIN REFORMATTED C17
############################################################################################################################################################

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

Conj:=proc(L) option remember: local k, C1, i1, L1, i: if L=[] then RETURN([]); end if; k:=nops(L): L1:=[seq(L[i]-1,i=1..k)]:
for i1 from 1 to nops(L1) while L1[i1]>0 do end do; i1:=i1-1: L1:=[op(1..i1,L1)]: C1:=Conj(L1): [k,op(C1)]; end proc:

#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 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); end if;
if not(j>=1 and j<=L[i]) then RETURN(FAIL); end if; C:={seq([i,j1],j1=j..L[i])}: for i2 from i to k while j<=L[i2] do end do;
i2:=i2-1: C:=C union {seq([i1,j],i1=i..i2)}: end proc:

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

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

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

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

############################################################################################################################################################
#END REFORMATTED C17

#1. [A little bit challenging, do your best] Write a procedure RandSYT(L) that inputs a partition L and outputs a uniformly-at-random standard Young tableau of shape L,
#   by implementing the The Greene-Nijenhuis-Wilf proof of the Hook-Lenght formula
#   Hints: 1. If you have a set S, to get a uniformly-at-random member of S do S[rand(1..nops(S))()].
#          2. Write a procedure OneTrial, that inputs a partition L and picks a random corner-cell, following the "kick the buck" algorithm of the paper.
#             Let L1 be the shape obtained by removing that corner, then recursively find RandSYT(L1) and put "n" in the removed corner. You may have to add a new row with only n in it.

#Picks a random corner of L
RCorner:=proc(L):
local C, rcell, hookset, rhook:
C:=Cells(L):
rcell:=C[rand(1..nops(C))()]:
hookset:=Hook(L,rcell) minus {rcell}:
if nops(hookset) < 1 then return(rcell); end if;
while nops(hookset) > 0 do
rhook:=hookset[rand(1..nops(hookset))()]:
hookset:=Hook(L,rhook) minus {rhook}:
end do;
rhook;
end proc:

#Returns L without a random corner along with the corner that was removed
OneTrial:=proc(L):
local corner, row, L1:
corner:=RCorner(L): row:=corner[1]:
if op(L[row])-1 > 0 then
L1:=[op(L[1..row-1]),op(L[row])-1,op(L[row+1..nops(L)])];
else
L1:=[op(L[1..row-1])];
end if;
return([L1,corner]);
end proc:

RandSYT:=proc(L):
local L1, sum, i, T, hold, n:
n:=nops(L):
T:=[seq([0$L[i]],i=1..n)]:
sum:=0: for i from 1 to n do
sum:=sum + op(L[i]): end do; L1:=L:
while nops(L1) > 0 do
hold:=OneTrial(L1):
T[hold[2][1]][hold[2][2]]:=sum:
sum:=sum - 1:
L1:=hold[1]:
end do;
T;
end proc:

#RandSYT([3,3,3]);
#Verified working

#2. By running RandSYT([5$5]) 10000 times, estimate the probability that "2" is in location [1,2] (rather than [2,1]).

Experiment17_2:=proc():
local i, count, L:
count:=0:
for i from 1 to 10000 do
L:=RandSYT([5$5]);
if op(L[1][2])=2 then count:=count + 1; end if;
end do;
count/10000;
end proc:

#Experiment17_2();
#Outputs: 2489/5000

#3. [Optional challenge, 5 dollars] Find the exact value of that probability.

#The exact probability is trivially 1/2, 2 can only be placed in [1,2] or [2,1] as 1 is always placed in [1,1]
#and placing 2 anywhere else would imply the existence of a natural number between 1 and 2.
#Since 2 is being placed uniformly at random under the construction algorithm the probability it's in any
#of it's possible places is 1/#of places = 1/2. QED