#OK to post
#Yuxuan Yang, April 24th, Assignment 24

with(combinat):

#2

TabPairs:=proc(n,m) local i,p,S,T,L:
p:={}:
for L in Pn(n) do
 for S in TAB(L,m) do
  for T in SYT(L) do
   p:=p union {[S,T]}:
  od:
 od:
od:
p:
end:

#3
#I think George's InvRS() works generally
#The key point is in procedure invert1()

# for j from 2 to nops(T[i]) while T[i][j] < a do od:

#this line find the position T[i][j] which is the first one >= a
#then T[i][j-1] is the position where we did the operation.































#C11.txt: Maple Code for Lecture 11 of Dr.Z's Experimental Mathematics class, Spring 2021, taught ar Rutgers University

Help11:=proc(): print(`PnClass(n), Pnk(n,k), Pn(n),pnk(n,k), pn(n) `): 
print(`pnSeq(N), pnSeqGF(N) `):
end:

#PnClass(n): THE SET of integer partitions of n
PnClass:=proc(n) local n1,L,L1,j:
option remember:
if n=0 then 
 RETURN({[]}):
elif n<0 then
 RETURN({}):
else
L:={}:

  for n1 from 1 to n do
   L1:=PnClass(n-n1):
   for j from 1 to nops(L1) do
     L:= L union {sort([op(L1[j]),n1],`>`)}:
   od:
 od:
   RETURN(L):
fi:
end:


#Pnk(n,k): The LIST of integer partitions of n with largest part k
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:

#pnk(n,k):
pnk:=proc(n,k) local k1:
option remember:
if not (type(n,integer) and type(k,integer) and n>=1 and k>=1 )then
 RETURN(0):
fi:
if n=k then
 RETURN(1):
else
 RETURN(add(pnk(n-k,k1),k1=1..min(k,n-k))):
fi:
end:

pn:=proc(n) local k: option remember: 
if n=0 then
 RETURN(1):
fi:

add(pnk(n,k),k=1..n):end:

pnSeq:=proc(N) local n: [seq(pn(n),n=0..N)]:end:

#pnSeqGF(N):
pnSeqGF:=proc(N) local i,q,f:
f:=taylor(mul(1/(1-q^i),i=1..N),q=0, N+1 ):
[seq(coeff(f,q,i),i=0..N)]:
end:

#From C10.txt
#C10.txt: Maple code for Lecture 10
Help10:=proc(): print(` SAc(f,x,K), LD(L) ,PlotDist(f,x,K) , OneTrial(p,n)`): 
print(`ManyTrials(p,n,K,x)`):
end:

#Taken from Blair Seidler's homework (who adapted a previous
#code from Dr. Z.)
#SAc(f,x,K): The SA(X,x,K) of the rv X such whose Pgf is f of x
SAc:=proc(f,x,K) local i,mu,sig2,M,f1:
if subs(x=1,f)=0 then
 RETURN(0):
fi:

f1:=f/subs(x=1,f):
mu:=subs(x=1,diff(f1,x)):
f1:=f1/x^mu:
f1:=x*diff(f1,x):
f1:=x*diff(f1,x):
sig2:=subs(x=1,f1):

M:=[mu,sig2]:

if sig2=0 then RETURN(M): fi:
for i from 3 to K do
  f1:=x*diff(f1,x):
#  M:=[op(M), evalf(subs(x=1,f1)/sig2^(i/2))]:
  M:=[op(M), subs(x=1,f1)/sig2^(i/2)]:
od:
M:
end:

#LD(L): Inputs a list of positive integers L (of n:=nops(L) members)
#outputs an integer i from 1 to n with the prob. of i being
#proportional to L[i]
#For example LD([1,2,3]) should output 1 with prob. 1/6
#output 2 with prob. 1/3
#output 3 with prob. 3/6=1/2
LD:=proc(L) local n,i,su,r:
n:=nops(L):
r:=rand(1..convert(L,`+`))():
su:=0:

for i from 1 to n do
 su:=su+L[i]:
  if r<=su then
    RETURN(i):
  fi:
od:
end:

#PlotDist(f,x,K): Inputs a prob. generating function of
#some r.v. X outputs a graph of its SCALED distribution
#A graph of the pdf of (X-mu)/sig. It draws it
#from -K*sig to K*sig

PlotDist:=proc(f,x,K) local i,mu,sig,M,f1:
if subs(x=1,f)=0 then
 RETURN(0):
fi:

f1:=f/subs(x=1,f):
mu:=subs(x=1,diff(f1,x)):
f1:=f1/x^mu:
f1:=x*diff(f1,x):
f1:=x*diff(f1,x):
sig:=sqrt(subs(x=1,f1)):
if sig=0 then RETURN(FAIL): fi:

#corrected after class, following Blair's suggestion
plot([seq([(i-mu)/sig,coeff(f,x,i)/subs(x=1,f)*sig],i=trunc(mu-K*sig)..trunc(mu+K*sig))]):

end:

#OneTrial(p,n): The outcome of ONE trial of tossing a coin n times
#where Pr(Head)=p. Returns the number of Heads
#p=1/3
OneTrial:=proc(p,n) local L,p1:
p1:=convert(p,rational):
L:=[denom(p1)-numer(p1), numer(p1)]:
add(LD(L)-1,i=1..n):
end:

#ManyTrials(p,n,K,x): the EMPIRICAL PGF  in the variable x, of repeating K times
#OneTrial(p,n)
ManyTrials:=proc(p,n,K,x) local i:
add( x^OneTrial(p,n) , i=1..K)/K:
end:



#C22.txt
Help22:=proc(): print(` Pn(n), AP(Y,i,m), SYT(L)  `): end:

#Ap(Y,i,m):write a procedure that inputs a Standard Young Tableau
#and an integer i between 1 and nops(Y)+1 and inserts the
#entry m at  the end of the i-th row, or if i=nops(Y)+1 make
#a new row with 1 cell occupied with m
AP:=proc(Y,i,m):

if i=nops(Y)+1 then
 RETURN([op(Y),[m]]):
fi:


if i>1 and nops(Y[i])=nops(Y[i-1]) then
  RETURN(FAIL):
fi:

[op(1..i-1,Y),[op(Y[i]),m],op(i+1..nops(Y),Y)]:


end:

#SYT(L): inputs an integer partion (given as a weakly decreasing
#list of POSITIVE integers OUTPUTS the LIST of
#Standard Young Tableax of shape L (n=Number of boxes of L)
#(A way of filling 1, ..., n in the boxes such that horiz. and
#vertically they are increasing
SYT:=proc(L) local n,L1,S,S1,i,j:
option remember:
n:=add(L):

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

S:=[]:

for i from 1 to nops(L)-1 do
 if L[i]>L[i+1] then
  L1:=[op(1..i-1,L),L[i]-1,op(i+1..nops(L),L)]:
  S1:=SYT(L1):
 S:=[op(S), seq(AP(S1[j],i,n),j=1..nops(S1))]:
 fi:
od:


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

 S1:=SYT(L1):
 S:=[op(S), seq(AP(S1[j],nops(L),n),j=1..nops(S1))]:

S:

end:


#syt(L): inputs an integer partion (given as a weakly decreasing
#list of POSITIVE integers OUTPUTS the NUMBER of
#Standard Young Tableax of shape L (n=Number of boxes of L)
#(A way of filling 1, ..., n in the boxes such that horiz. and
#vertically they are increasing
syt:=proc(L) local n,L1,S,S1,i,j:
option remember:
n:=add(L):

if n=1 then 
  RETURN(1):
fi:

S:=0:

for i from 1 to nops(L)-1 do
 if L[i]>L[i+1] then
  L1:=[op(1..i-1,L),L[i]-1,op(i+1..nops(L),L)]:
  S:=S+syt(L1):
 
fi:
od:


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

 S:=S+syt(L1):
S:
end:

##FROM C11.txt
#Pnk(n,k): The LIST of integer partitions of n with largest part k
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:
#END FROM C11.txt

#C24.txt: Maple code for Lecture 24#OK to post homework

Help24:=proc(): print(`SYTpairs(n), RS(pi), InvRS(P) , Words(m,n)`): 
print(`APP(Y,b,m), Box(L) , NewBox(B) TAB(L,m) `):
end:

#APP(Y,b,m): inputs a tableaux Y with entires {1,...,m-1} and
#and a vector b of either the same length of Y or one longer and appends
#b[i] m's to the end of row i for i=1..nops(Y) and if nops(b)=nops(Y)+1 a new row of b[nops(b)] m's. For example
#APP([[1,1,2],[2,2]],[1,1,1],3) should be
#[[1,1,2,3],[2,2,3],[3]]
APP:=proc(Y,b,m) local i,Y1:
Y1:=[seq([op(Y[i]),m$b[i]],i=1..nops(Y))]:
if nops(b)=nops(Y)+1 then
 Y1:=[op(Y1),[m$b[-1]]]:
fi:
Y1:
end:

#TAB(L,m): the set of all tableaux of shape L with entries
#{1,2,...,m}
TAB:=proc(L,m)  local L1,B,i,S,b,y,S1:
option remember:

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

if m<nops(L) then
  RETURN({}):
fi:

if nops(L)=1 and m=1 then
 RETURN({[[1$L[1]]]}):
fi:


B:=NewBox([seq(L[i]-L[i+1],i=1..nops(L)-1),L[nops(L)]]):

S:={}:

for b in B do
 L1:=[seq(L[i]-b[i],i=1..nops(L))]:
 if L1[-1]=0 then
  L1:=[op(1..nops(L1)-1,L1)]:
fi:

 S1:=TAB(L1,m-1):

S:=S union {seq(APP(y,b,m),y in S1)}:

od:
S:

end:

#NewBox(B): inputs a list of NON-NEGATIVE integers and outputs
#the set of all lists [i1,i2,...] i1<=B[1], i2<=B[2] ,...
NewBox:=proc(B) local k,B1,S1,i,v:
option remember:
if B=[]  then
 RETURN({[]}):
fi:

k:=nops(B):
B1:=[op(1..k-1,B)]:

S1:=NewBox(B1):
{seq(seq([op(v),i], v in S1),i=0..B[k])}:

end:

#Box(L): Inputs a list L (where k=nops(L)) of non-negative integers L outputs
#the  LIST of all LISTS [a1,a2,..., ak] in LEX ORDER such that ai is in {0,1,.., L[i])
#In particular the n-dim UNIT cube would be Box([1$n]);
Box:=proc(L) local k,i,L1,B1,j:
option remember:

if not type(L,list) then
  print(L, `should be a list `):
 RETURN(FAIL):
fi:

k:=nops(L):
if k=0 then 
 RETURN([[]]):
fi:
 if not ({seq(type(L[i],integer),i=1..nops(L))}={true} and min(op(L))>=0) then
   print(`bad input`):
  RETURN(FAIL):
 fi:




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

B1:=Box(L1):
 [ seq(seq([op(B1[j]),i],i=1..L[k]),j=1.. nops(B1))]:
end:



#From George Spahn, 4/18/2021, Assignment 23

#2

# SYTpairs(n) inputs a positive integer n and outputs the set of pairs of 
# standard Young tableaux of the same shape [S,T] each with n cells

SYTpairs:=proc(n) local s,L,T,i,j:

L:={}:
for s in Pn(n) do:
 T:=SYT(s):
 L := L union {seq(seq([T[i],T[j]] ,j=1..nops(T)) ,i=1..nops(T))}:
od:
L:
end: 

###############################################################################
#3


#RS1m(Y,m): ONE STEP OF THE RS algorithm
# modified to also return the row of the insertion
RS1m:=proc(Y,m) local i,Y1:

Y1:=Ins(Y,1,m):
if Y1[2]=0 then
  RETURN(Y1[1],1):
fi:

#continued after class
for i from 2 to nops(Y)+1 while Y1[2]<>0 do
 Y1:=Ins(Y1[1],i,Y1[2]):
od:
Y1[1],i-1:
end:

#RS1mOld(Y,m): ONE STEP OF THE RS algorithm
# modified to also return the row of the insertion
RS1mOld:=proc(Y,m) local i,Y1:

Y1:=Ins(Y,1,m):
if Y1[2]=0 then
  RETURN(Y1[1],1):
fi:

#continued after class
for i from 2 to nops(Y)+1 while Y1[2]<>0 do
 Y1:=InsOld(Y1[1],i,Y1[2]):
od:
Y1[1],i-1:
end:

#RS(pi):  the pair [S,T] that is the output
#of the Robinson-Scheastead mapping
RS:=proc(pi) local Y1,Y2,i,row:
Y1:=[[pi[1]]]:
Y2:=[[1]]:
for i from 2 to nops(pi) do
 Y1,row:=RS1m(Y1,pi[i]):
 if row = nops(Y2)+1 then 
  Y2:=[op(Y2),[i]]:
 else 
  Y2[row]:=[op(Y2[row]),i]:
 fi:
od:
[Y1,Y2]:
end:

#############################################################################
#4

invert1 := proc(T1,row) local T,a,i,j,c:
T:=T1:
i:=row:
a := T[i][-1]:
if nops(T[i]) = 1 then 
 T:=[op(1..nops(T)-1, T)]:
else
 T[i]:=[op(1..nops(T[i])-1, T[i])]:
fi:
while i>1 do 
 i:=i-1:
 for j from 2 to nops(T[i]) while T[i][j] < a do od:
 c:=T[i][j-1]:
 T[i][j-1]:=a:
 a:=c:
od:
T,a:
end:

# InvRS(P) inputs a pair [S,T] of standard Young tableaux of the same shape 
# and outputs a permutation of 1,2,...n
InvRS:=proc(P) local S,T,perm,n,x,j,a,i:
S:=P[1]:
T:=P[2]:
perm := []:

n:= max( seq(max(T[i]),i=1..nops(T)) ):
for x from n to 1 by -1 do

 for j from 1 to nops(T) while not(x in T[j]) do od:
 S,a:=invert1(S,j):
 perm:=[a,op(perm)]:  

od:
perm:
end:

##############################################################################
#old stuff

with(combinat):

#AllSYT(n): The set of all Standard Young Tableau with n cells
AllSYT:=proc(n) local S,i:
S:=Pn(n):
{seq(op(SYT(S[i])),i=1..nops(S))}:
end:

#InsOld(Y,i,m): Inputs a (partially filled tableux Y,
#a row i, between 1 and nops(Y)+1 and an integer m
#and tries to place it LEGALLY at row i
#(if it is larger than all the current entries of row i
#and does not stick out, i=1 OR nops(Y[i])<nops(Y[i-1])
#and outputs the new tableau and either 0 or the new bumpbee
InsOld:=proc(Y,i,m) local j:
if not (1<=i and i<=nops(Y)+1) then
 RETURN(FAIL):
fi:

if i=nops(Y)+1 then
  RETURN([op(Y),[m]],0):
fi:

if i=1 then 
  if m>=Y[1][-1] then
    RETURN([[op(Y[1]),m], op(2..nops(Y),Y)],0):
   fi:
fi:

for j from 1 to nops(Y[i]) while m>Y[i][j] do od:

if j=nops(Y[i])+1  and (i=1 or nops(Y[i])<nops(Y[i-1])) then

#corrected after class , Maple does not like op(m..n,L) if m-n>=2
if j+1<=nops(Y[i]) then
  RETURN ([op(1..i-1,Y),[op(1..j-1,Y[i]),m,op(j+1..nops(Y[i]),Y[i])],
op(i+1..nops(Y),Y)],0):
else

RETURN ([op(1..i-1,Y),[op(1..j-1,Y[i]),m],
op(i+1..nops(Y),Y)],0):
fi:


else
RETURN([op(1..i-1,Y),[op(1..j-1,Y[i]),m,op(j+1..nops(Y[i]),Y[i])],
op(i+1..nops(Y),Y)],Y[i][j]):
fi:

end:

#Ins(Y,i,m): Inputs a (partially filled tableux Y,
#a row i, between 1 and nops(Y)+1 and an integer m
#and tries to place it LEGALLY at row i
#(if it is larger than all the current entries of row i
#and does not stick out, i=1 OR nops(Y[i])<nops(Y[i-1])
#and outputs the new tableau and either 0 or the new bumpbee
Ins:=proc(Y,i,m) local j:
if not (1<=i and i<=nops(Y)+1) then
 RETURN(FAIL):
fi:

if i=nops(Y)+1 then
  RETURN([op(Y),[m]],0):
fi:

if i=1 then 
  if m>=Y[1][-1] then
    RETURN([[op(Y[1]),m], op(2..nops(Y),Y)],0):
   fi:
fi:

for j from 1 to nops(Y[i]) while m>=Y[i][j] do od:

if j=nops(Y[i])+1  and (i=1 or nops(Y[i])<nops(Y[i-1])) then

#corrected after class , Maple does not like op(m..n,L) if m-n>=2
if j+1<=nops(Y[i]) then
  RETURN ([op(1..i-1,Y),[op(1..j-1,Y[i]),m,op(j+1..nops(Y[i]),Y[i])],
op(i+1..nops(Y),Y)],0):
else

RETURN ([op(1..i-1,Y),[op(1..j-1,Y[i]),m],
op(i+1..nops(Y),Y)],0):
fi:


else
RETURN([op(1..i-1,Y),[op(1..j-1,Y[i]),m,op(j+1..nops(Y[i]),Y[i])],
op(i+1..nops(Y),Y)],Y[i][j]):
fi:

end:

#SYT(L): inputs an integer partion (given as a weakly decreasing
#list of POSITIVE integers OUTPUTS the LIST of
#Standard Young Tableax of shape L (n=Number of boxes of L)
#(A way of filling 1, ..., n in the boxes such that horiz. and
#vertically they are increasing
SYT:=proc(L) local n,L1,S,S1,i,j:
option remember:
n:=add(L):

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

S:=[]:

for i from 1 to nops(L)-1 do
 if L[i]>L[i+1] then
  L1:=[op(1..i-1,L),L[i]-1,op(i+1..nops(L),L)]:
  S1:=SYT(L1):
 S:=[op(S), seq(AP(S1[j],i,n),j=1..nops(S1))]:
 fi:
od:


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

 S1:=SYT(L1):
 S:=[op(S), seq(AP(S1[j],nops(L),n),j=1..nops(S1))]:

S:

end:


##FROM C11.txt
#Pnk(n,k): The LIST of integer partitions of n with largest part k
Pnk:=proc(n,k) local k1,L,L1,j:
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:
#END FROM C11.txt