Help := proc();
print(`Mamas(L), YT(L), YTixi(i),YTaxb(a,b),YTaxbStrict(a,b),Cat(n)`);
print(`Phi1(YT),Phi2(vec),Phi(YT),Phi1Inv(vec)`);
print(`Par(n), Par1(n,k), ParBox(n,i,j), ParBox1(n,k,l)`);
print(`F(L), Fset(shapes)`);
print(`AllixiParts(i), AllaxbParts(a,b),AllaxbPartsStrict(a,b)`);
print(`GuessRat(L,x), GuessRat1(L,x,d), GuessPol(L,n), GuessPol1(L,d,n)`);
end:

with(plots):
with(combinat):

#========================================================
#Cat(n): the nth Catalan number

Cat := proc(n)

RETURN(binomial(2*n,n)/(n+1));

end:

#========================================================
#Mamas(L): all the shapes obtained from L by
#nibbling a corner
Mamas:=proc(L) local i,k,S:

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

S:={}:
k:=nops(L):

for i from 1 to k-1 do
 if L[i]>L[i+1] then
  S:=S union {[op(1..i-1,L),L[i]-1,op(i+1..k,L)]}:
 fi:
od:

if L[k]>=2 then
   S:=S union {[op(1..k-1,L),L[k]-1]}:
else
  S:=S union {[op(1..k-1,L)]}:
fi:

S:

end:

#========================================================
#YT(L): inputs a partition shape L and outputs the number
#of standard Young tableaux of shape L using the "going
#down" formula

YT:=proc(L) local S,s:

option remember:

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

S:=Mamas(L):

add(YT(s), s in S):
end:

#========================================================
#F(L):the set of standard YT of shape L

F:=proc(L) 
local A,L1,n,A1,i,m:

option remember:

m:=nops(L):
n:=convert(L,`+`):

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

A:={}:

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)]:
		A1:=F(L1):
		A1:={seq( [ op(1..i-1,a1), [op(a1[i]),n],op(i+1..nops(a1),a1)], a1 in A1)}:
		A:=A union A1:
	fi:
od:

if L[m]>1 then
	L1:=[op(1..m-1,L),L[m]-1]:
	A1:=F(L1):
	A1:= {seq( [ op(1..m-1,a1), [op(a1[m]),n]], a1 in A1)}:
	A:=A union A1:
else
	L1:=[op(1..m-1,L)]:
	A1:=F(L1):
	A1:= {seq( [ op(1..m-1,a1),[n] ], a1 in A1)}:
	A:=A union A1:
fi:

RETURN(A):

end:

#========================================================
#Fset(shapes): the set of standard young tableaux of shape
#L for all L in shapes

Fset := proc(shapes);

RETURN({seq(op(F(L)), L in shapes)});

end:

#========================================================
#Par(n): the set of integer-partitions of n
Par:=proc(n) local k:
option remember:

{seq( op(Par1(n,k)),k=1..n)}:

end:

#========================================================
#Par1(n,k): the set of partitions of n whose
#largest part is k
Par1:=proc(n,k) local k1,beth,b,S:
option remember:

if n=1 then
 if k=1 then
   RETURN({[1]}):
 else
   RETURN({}):
 fi:
fi:

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

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

S:={}:

for  k1 from 1 to k do
 beth:=Par1(n-k,k1):
 S:=S union {seq([k,op(b)], b in beth)}:
od:
 
S:
end:

#========================================================
#ParBox(n,i,j): the set of partitions of n whose largest
#part is <=i and number of parts is <=j, i.e- they fit in
#and ixj (i=cols, j=rows) box

ParBox := proc(n,i,j) option remember;
local k,l;

RETURN({seq(seq(op(ParBox1(n,k,l)), l=1..j),k=1..i)});

end:

#========================================================
#ParBox1(n,k,l): the set of partitions of n whose largest
#part is k and has l parts

ParBox1 := proc(n,k,l) option remember;
local smaller,RET,L;

if n=1 then
	if k=1 then
		if l=1 then
			RETURN({[1]});
		fi;
	fi;
fi;

if k=1 then
	if l=n then
		RETURN({[seq(1,i=1..n)]});
	fi;
fi;

if l=1 then
	if k=n then
		RETURN({[n]});
	fi;
fi;

if k*l<n or k>n or l>n or k<=0 or l<=0 or n<=0 then
	RETURN({});
fi;

smaller := {seq(op(ParBox1(n-k,k1,l-1)), k1=1..k)};
RET := {};

for L in smaller do
	RET := RET union {[k,op(L)]};
od;

RETURN(RET);

end:

#========================================================
#AllixiParts(i): outputs the list of partitions fitting
#inside an i by i box.

AllixiParts := proc(i)
local n;

RETURN({seq(op(ParBox(n,i,i)),n=1..i^2)});

end:

#========================================================
#AllaxbParts(a,b): outputs the list of partitions fitting
#inside an a (cols) by b (rows) box.

AllaxbParts:=proc(a,b)
local n;

RETURN({seq(op(ParBox(n,a,b)),n=1..a*b)});

end:

#========================================================
#AllaxbPartsStrict(a,b): outputs the list of partitions fitting
#inside an a (cols) by b (rows) box and no smaller box.
#i.e.- the largest part is exactly a and there are exactly
#b parts

AllaxbPartsStrict := proc(a,b)
local parts,L;

parts := AllaxbParts(a,b);
for L in parts do
	if L[1]<>a or nops(L)<>b then
		parts := parts minus {L};
	fi;
od;

RETURN(parts);

end:

#========================================================
#YTixi(i): the number of Young Tableaux fitting inside an
#i by i box

YTixi := proc(i)
local shapes,poss,L,n;

shapes :=AllixiParts(i);
poss := 0;
for L in shapes do
	poss := poss+YT(L);
od;

RETURN(poss);

end:

#========================================================
#YTaxb(a,b): the number of Young Tableaux fitting inside an
#a (columns) by b (rows) box
YTaxb := proc(a,b)
local shapes,poss,L,n,EX;

shapes := AllaxbParts(a,b);
poss := 0;
for L in shapes do
	poss := poss + YT(L);
od;

RETURN(poss);

end:

#========================================================
#YTaxbStrict(a,b): the number of Young tableaux fitting 
#inside an a (columns) by b (rows) box and not any smaller
#box. So the number of rows must equal b and the number of
#columns must equal a. Whatever happens in between is fine

YTaxbStrict := proc(a,b)
local shapes, poss,L,n,shapes2;

shapes := AllaxbPartsStrict(a,b);
poss := 0;
#shapes2 := {};
for L in shapes do
	poss := poss + YT(L);
od;

RETURN(poss);
end:



#============================================================
#Recall that a rational function f(x) is a quotient of 
#polynomials f(x)=P(x)/Q(x). Write a program GuessRat1(L,x,d) 
#that inputs a list L, a variable x, and a non.neg. 
#integer d such that 2*d+5 <= nops(L), and outputs a 
#rational function P(x)/Q(x) with degree(P,x), 
#degree(Q,x) <= d, such that P(i)/Q(i)=L[i] for all 
#i from 1 to nops(L), if one exist, and FAIL otherwise. 

GuessRat1 := proc(L,x,d)
local P,Q,var,eq,i,Ans;

if 2*d+5>nops(L) then
	print(`List not big enough`);
	RETURN(FAIL);
fi;

P := add(a[i]*x^i, i=0..d);
Q := add(b[i]*x^i, i=0..d);

var := {seq(a[i], i=0..d),seq(b[j], j=0..d)};

eq := {};

for i from 1 to nops(L) do
	if not subs(x=i,Q)=0 then
		eq := eq union {subs(x=i,P/Q)-L[i]};
	else
		print(`denominator equals 0 for i=`,i);
		RETURN(FAIL);
	fi;
od;

Ans := solve(eq,var);

if Ans=NULL then
	RETURN(FAIL);
fi;

RETURN(factor(subs(Ans,P/Q)));

end:

#============================================================
#By trying d=1..(nops(L)-5)/2, write a program GuessRat(L,x) 
#that guesses a rational function P(x)/Q(x) that fits L.

GuessRat := proc(L,x)
local d,PoQ;

for d from 1 to (nops(L)-5)/2 do
	PoQ := GuessRat1(L,x,d);
	if not PoQ = FAIL then
		RETURN(PoQ);
	fi;
od;

RETURN(FAIL);

end:

#============================================================
#GuessPol1(L,d,n): guesses a polynomial of degree d in n for
# the list L, such that L[i]=P[i] for i>=s0 for example, try: 
#GuessPol1([(seq(i,i=1..10),1,n,1);
GuessPol1:=proc(L,d,n) local P,i,a,eq,var:
if d>=nops(L)-3 then
 ERROR(`the list is too small`):
fi:

P:=add(a[i]*n^i,i=0..d):
var:={seq(a[i],i=0..d)}:
eq:={seq(subs(n=i,P)-L[i],i=1..d+4)}:

var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
fi:

subs(var,P):

end:


#============================================================
#GuessPol(L,n): guesses a polynomial of degree d in n for
# the list L, such that L[i]=P[i] for i>=1 for example, try: 
#GuessPol([(seq(i,i=1..10),n);
GuessPol:=proc(L,n) local d,gu:

for d from 0 to nops(L)-4 do
 gu:=GuessPol1(L,d,n):
 if gu<>FAIL then
    RETURN(gu):
 fi:
od:

FAIL:

end:

GuessExp := proc(L,n)
local var,eq,ans;

var := {a,b};
eq := {seq(a*b^i-L[i], i=1..nops(L))};

ans := solve(eq,var);
if ans = NULL then
	RETURN(FAIL);
else
	RETURN(ans);
fi;

end:

#========================================================
#Phi1(YT): inputs a young tableaux YT of 2 rows and 
#outputs a vector of length 2*(nops(YT[1])+1) 1 and -1 
#where no partial sum is <0. This is described
#in the write up.

Phi1 := proc(YT)
local n,m,vec,i;

n := nops(YT[1]);
m := nops(YT[2]);

vec := [];
for i from 1 to n+m do
	if i in YT[1] then
		vec := [op(vec),1];
	elif i in YT[2] then
		vec := [op(vec),-1];
	fi;
od;

vec := [op(vec),1,seq(-1,i=n+m+2..2*n+2)];

RETURN(vec);

end:

#========================================================
#Phi2(vec): inputs a vector of length 2*n and plots the
#corresponding Dyck path

Phi2 := proc(vec)
local i,points,n1;

points := [[0,0]];

for i from 1 to nops(vec) do
	n1 := nops(points);
	if vec[i]=1 then
		points := [op(points),[points[n1][1],points[n1][2]+1]];
	elif vec[i]=-1 then
		points := [op(points),[points[n1][1]+1,points[n1][2]]];
	fi;
od;

RETURN(PLOT(POINTS(op(points)),CURVES(points),AXESSTYLE(BOXED),SYMBOL(CIRCLE,15)));

end:

#========================================================
#Phi1Inv(vec): inputs a vector of length 2*n and outputs
#a young tableaux of 2 rows with top row length n-1.
#The algorithm is described in the write up.

Phi1Inv := proc(vec)
local YT,oneFlag,loc,i,FLAG;

YT:=[[1],[]];
oneFlag := 1;
FLAG:=true;
loc := 1;

for i from 2 to nops(vec) while FLAG do
	if vec[i]=1 then
		oneFlag := oneFlag+1;
		if oneFlag > nops(vec)/2-1 then
			FLAG:=false;
		else
			YT:=[[op(YT[1]),i],YT[2]];
		fi;
	elif vec[i]=-1 then
		YT:=[YT[1],[op(YT[2]),i]];
	fi;
od;

RETURN(YT);

end:

#========================================================
#Phi(YT): Phi2(Phi1(YT)). Inputs a YT and outputs its 
#corresponding Dyck path as described in the write up.

Phi := proc(YT);

RETURN(Phi2(Phi1(YT)));

end: