#OK to post homework
#Blair Seidler, 2021-04-14 Assignment 22

with(combinat):

Help:=proc(): print(`randomSYT(L),TotalSYT(n),TotalSYTpower(n,k)`): end:

# 1.

#randomSYT(L): inputs an integer partition L (aka shape) and outputs, uniformly at random,
#one of the syt(L) Standard Young tableaux of that shape
randomSYT:=proc(L) local L1,die,S1,n,i:
n:=add(L):
if n=1 then RETURN([[1]]): fi:
die:=[]:
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)]:
    die:=[op(die),syt(L1)]:
  else
    die:=[op(die),0]:
  fi:
od:

if L[-1]>1 then
  L1:=[op(1..nops(L)-1,L),L[i]-1]:
else
  L1:=[op(1..nops(L)-1,L)]:
fi:
die:=[op(die),syt(L1)]:
i:=LD(die):
L1:=[op(1..i-1,L),L[i]-1,op(i+1..nops(L),L)]:
if L1[-1]=0 then L1:=[op(1..nops(L1)-1,L1)]: fi:
S1:=randomSYT(L1):
Ap(S1,i,n):
#S:

end:

# 2.

#TotalSYT(n): inputs a positive integer n, and outputs the total number of standard Young
#tableaux with n cells [Hint: sum-up syt(L) over all members of Pn(n)
TotalSYT:=proc(n) local P,p:
if n=0 then RETURN(1): fi:
P:=Pn(n):
add(syt(p),p in P):
end:

# 3.

#TotalSYTpower(n,k): inputs a positive integer n, and outputs the total number of k-tuples of
#standard Young tableaux of the same shape with n cells#[Hint: sum-up syt(L)k over all members of Pn(n)
TotalSYTpower:=proc(n,k) local P,p:
if n=0 then RETURN(1): fi:
P:=Pn(n):
add(syt(p)^k,p in P):
end:


#Is that sequence in the OEIS? What is its A-number
#Yes: 1, 2, 4, 10, 26, 76, 232, 764, 2620, 9496, 35696, 140152, 568504, 2390480,
#     10349536, 46206736, 211799312, 997313824, 4809701440, 23758664096
# Is A000085 in the OEIS

#### From C22.txt ####
#Ap(Y,i,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:

#k=2 is A000142 in OEIS
#k=3 is A130721
#k=4 is A129627
#k=5 is A218432
#k=6 is A218433

#SYT(L) inputs an integer partition given as a weakly decr. list of pos integers
#outputs the list of Standard Young Tableaux of shape L
SYT:=proc(L) local n,L1,S1,S,i:
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[i]-1]:
else
  L1:=[op(1..nops(L)-1,L)]:
fi:
S1:=SYT(L1):
S:=[op(S), seq(Ap(S1[j],i,n),j=1..nops(S1))]:

S:

end:

#syt(L) inputs an integer partition given as a weakly decr. list of pos integers
#outputs the number of Standard Young Tableaux of shape L
syt:=proc(L) local n,L1,S,i:
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[i]-1]:
else
  L1:=[op(1..nops(L)-1,L)]:
fi:
S:=S+syt(L1):

end:

#### From C11.txt ####
Help11:=proc(): print(`Pnk(n,k), Pn(n)`): 
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:

#### From C4.txt ####
#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: