#Please do not post homework
#AJ Bu, April 12 2021, Assignment 22

########################	PROBLEM 1	########################
#randomSYT(L):  inputs an integer partition L (aka shape) 
#outputs, uniformly at random, one of the syt(L) Standard Young tableaux of that shape

randomSYT:=proc(L) local n,S,die1,s,m:
	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
			S:=[op(S),[op(1..i-1,L),L[i]-1,op(i+1..nops(L),L)]]:
		fi:
	od: 
		if L[-1]>1 then
			S:=[op(S),[op(1..nops(L)-1,L),L[nops(L)]-1]]:
		else	
			S:=[op(S),[op(1..nops(L)-1,L)]]:
		fi:
		
	die1:=[seq(syt(s), s in S)]:  
	L1:=S[LD(die1)]: 
	if nops(L)=nops(L1) then
		m:=max[index](L-L1):
	else
		m:=nops(L):
	fi:
	S1:=randomSYT(L1):  
	pi:=AP(S1,m,n):  
	
end:


########################	PROBLEM 2	########################
#TotalSYT(n): inputs a positive integer n
#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,L:
	P:=Pn(n):
	add(syt(L), L in P):
end:

#Is that sequence in the OEIS? What is its A-number?

##seq(TotalSYT(i),i=1..10)
###        1, 2, 4, 10, 26, 76, 232, 764, 2620, 9496
###	This is on OEIS: A000085 (except ours gives 0 for n=0 instead of 1)

########################	PROBLEM 3	########################
#TotalSYTpower(n,k): inputs a positive integer n
#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,L:
	P:=Pn(n):
	add(syt(L)^k, L in P):
end:

##seq(TotalSYTpower(i,2),i=1..10)
###      1, 2, 6, 24, 120, 720, 5040, 40320, 362880, 3628800
### This is on OEIS: A000142	


##seq(TotalSYTpower(i,3),i=1..10)
###      1, 2, 10, 64, 596, 8056, 130432, 2534960, 59822884, 1718480368
### This is on OEIS: A130721	

##seq(TotalSYTpower(i,4),i=1..10)
###      1, 2, 18, 180, 3060, 101160, 3807720, 174986280, 10699554600, 927701102160
### This is on OEIS: A129627	

##seq(TotalSYTpower(i,10),i=1..10)
###      1, 2, 1026, 119124, 82094580, 1126524259080, 5563004909321160, 43453047082604239080, 620787527477497337506920, 82539616591562766578923554000
### This is on OEIS: A218437	
	
########################	PROBLEM 4 (OPTIONAL CHALLENGES)	########################
#a. Conjecture a closed-form formula, in terms of a1, a2, for the number of Standard Young 
#tableaux of two rows of lengths a1,a2. In other words, an algebraic formula for syt([a1,a2]) 
#(where of course a1>= a2 >= 0).

#Some observations:
#For "a lot" of them have syt([a1,a2])= syt([a1-1,a2])+syt([a1,a2-1])
#syt([a1,a2])=syt([a1-1,a2])+syt([a1,a2-1]) 
#syt([a1,a1-1])=syt([a1,a1])= the a1-th Catalan number = binomial(2a1,a1)/(a1+1)
#So, syt([a1,a1])=binomial(a1+a1,a1)* 1/(a1+1)
#syt([a1,1])=a1=binomial(a1+1,a1)/(a1+1)*(a1+1-1)
#binomial(2a1,a1)=(a1+a1)!/(a1!a1!)=(a1+a1-1)!(a1+a1)/(a1!(a1-1)!a1)=2*binomial(a1+a1-1,a1)
#So, syt([a1,a1-1])=binomial(a1+a1-1,a1)* 2/(a1+1)

#CONJECTURE: syt([a1,a2])=binomial(a1+a2,a1)*(a1+1-a2)/(a1+1).



########################	FROM C22.txt 	########################
#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




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