#OK to post
#Yuxuan Yang, April 12th, Assignment 22

with(combinat):

with(GroupTheory):

randomSYT:=proc(L) local n,L1,die1,i,SYT1:
n:=add(L):
if n=1 then 
  RETURN([[1]]):
fi:
die1:=[]:
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)]:
  die1:=[op(die1),syt(L1)]:
 else
  die1:=[op(die1),0]:
 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:
die1:=[op(die1),syt(L1)]:
i:=LD(die1):
if i=nops(L) and L[-1]=1 then 
 L1:=[op(1..nops(L)-1,L)]:
else
 L1:=[op(1..i-1,L),L[i]-1,op(i+1..nops(L),L)]:
fi:
SYT1:=randomSYT(L1):
AP(SYT1,i,n):
end:





TotalSYT:=proc(n) add(syt(L),L in Pn(n)): end:
#A000085 in OEIS



TotalSYTpower:=proc(n,k) add(syt(L) ^ k , L in Pn(n)): end:

#k=2 n! A000142
#k=3    A130721
#k=3    A129627





#########################################################
#syt([i,j])=binomial(i + j, j) - binomial(i + j, j - 1)##
#########################################################
#the proof is induction. It is just the same as the idea of proc syt.
#syt([i,j])=syt([i-1,j])+syt([i,j-1]) if i>=j+1
#syt([i,j])=syt([i,j-1]) if i=j
#these two are obvious by checking algebra.






#for 3, by some random experiments, I got
#syt([10, 8, 5])=
#((multinomial(23, 10, 8, 5) - multinomial(23, 11, 7, 5) - multinomial(23, 10, 9, 4) + multinomial(23, 11, 8, 4)) - multinomial(23, 10, 9, 4)) + multinomial(23, 11, 8, 4)
#conj: syt([a,b,c])=multinomial(a+b+c,a,b,c)-multinomial(a+b+c,a+1,b-1,c)+2*multinomial(a+b+c,a+1,b,c-1)-2*multinomial(a+b+c,a,b+1,c-1)
#when we have this, the proof is induction naturally.
#firstly,multinomial(a+b+c,a,b,c)=multinomial(a+b+c-1,a-1,b,c)+multinomial(a+b+c-1,a,b-1,c)+multinomial(a+b+c-1,a,b,c-1)
#then it is about boundary things.
#syt([a,a,c])=syt([a,a-1,c])+syt([a,a,c-1]) if a>c
#syt([a,b,b])=syt([a,b,b-1])+syt([a-1,b,b]) if a>b
#syt([a,a,a])=syt([a,a,a-1])
#but maple tells me it is INCORRECT
#The key point is, we need several multinomials and we need "syt"=0 when b=a+1 and c=b+1
#Then, I got a new formula
conj3syt:=proc(a,b,c) 
#############################################################################################
multinomial(a+b+c,a,b,c)-multinomial(a+b+c,a+1,b-1,c)-multinomial(a+b+c,a,b+1,c-1)+        ##
multinomial(a+b+c,a+1,b+1,c-2)+multinomial(a+b+c,a+2,b-1,c-1)-multinomial(a+b+c,a+2,b,c-2):## 
#############################################################################################
end:
#By checking algebra on the boundaries, this is correct. The proof is induction.


mymultinomial:=proc(L) local i:
for i from 1 to nops(L) do
 if L[i]<0 then 
  RETURN(0):
 fi:
od:
multinomial(op(L)):
end:

##############################################################################################
sytk:=proc(L) local i,k,pi:                                                                 ##
k:=nops(L):                                                                                 ##
add(PermParity(Perm(pi))*mymultinomial([add(L),seq(pi[i]-i+L[i],i=1..k)]),pi in permute(k)):##
end:                                                                                        ##
##############################################################################################


#the proof is checking all recursion situations.
#some of the ideas are in the pdf attached.






















#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



#C4.txt:
Help4:=proc(): print(` Snk(n,k), NuSnk(n,k) , LD(L), Rnk(n,k)  `): end:

#Snk(n,k): Inputs NON-NEG. integer n and an integer k and OUTPUTS
#THE SET OF k-element subsets of {1,...,n}
#FOR EXAMPLE
#Snk(3,2) should be {{1,2},{1,3},{2,3}}
Snk:=proc(n,k) local S1,S2,s:
option remember:

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

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

S1:=Snk(n-1,k):
S2:=Snk(n-1,k-1):
S1 union {seq( s union {n}, s in S2)}:
end:

#NuSnk(n,k): Inputs NON-NEG. integer n and an integer k and OUTPUTS
#THE NUMBER OF  k-element subsets of {1,...,n}
#FOR EXAMPLE
#NuSnk(3,2) should be 3
NuSnk:=proc(n,k) option remember:

if k<0 or k>n then
 RETURN(0):
fi:

if n=0 then
  if k=0 then
    RETURN(1):
  else
    RETURN(0):
  fi:
fi:

NuSnk(n-1,k)+ NuSnk(n-1,k-1):

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:

#Rnk(n,k): Inputs NON-NEG. integer n and an integer k and OUTPUTS
#A (UNIFORMLY AT-RANDOM)  k-element subsets of {1,...,n}
#FOR EXAMPLE
#Rnk(3,2) should output each of {1,2},{1,3},{2,3} with prob. 1/3

Rnk:=proc(n,k)  local c,S:
if k<0 or k>n then
 RETURN(FAIL):
fi:

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

c:=LD([NuSnk(n-1,k), NuSnk(n-1,k-1)]):

if c=1 then
  RETURN(Rnk(n-1,k)):
else
 RETURN(Rnk(n-1,k-1) union {n}):
fi:

end: