#core partitions, etc
#Anthony Zaleski, Spring 2016

print(`CORE PARTITIONS`):
print(`-----`):
print(`Execute Help(); to begin.`):

with(gfun):
with(combinat):

################
#HELP PROCEDURES
################


Help:=proc(): 
if nargs=0 then
	print(`core.txt contains the following primary procedures.`):
	print(`Type Help(<procedure name>); for help on a procedure.`):
	print(`Type Help1(); to list secondary procedures.`):
	#print(`DyckPaths, GoodWalks, GoodParts, NewGoodParts, NewGoodWalks`):
	#print(`GF1, GFNK, GFgw`):
	print(`Gs, CoreThms`):#,Gsg, GsTheo, GsRec`): 
	#print(`qGuessRec, qGuessRecOpt, qRec2Seq, GuessPol`):
else 
	LookupProc(args[1]):
fi:
end:


Help1:=proc(): 
if nargs=0 then
	print(`core.txt contains the following secondary procedures.`):
	print(`Type Help(<procedure name>); for help on a procedure.`):
	print(`Type Help(); to list primary procedures.`):
	print(`Gsg, GsTheo, GsRec`):
	print(`GFgw, GFgwN, GuessGFgwN, Gkl, Gklg, Fk`):
	print(`Nds, Gds`):
	print(`qGuessRec, qGuessRecOpt, qRec2Seq, GuessPol`):
	print(`GuessFibPol1, GuessFibPol, GuessFibPolList`):
	print(`Moms, gwMomList, gwStraightMoms, Straight2Central`):
	print(`qGuessRec1, GuessPol1`):
	print(`RepParts, Path2Part, OddParts, Idealize, NewGoodWalks1, FactorialMoment`):
	print(`qBin,F`):
	print(`abMoms `):
else 
	LookupProc(args[1]):
fi:
end:


LookupProc:=proc(a) local EG:
EG:=`For example, try`:
if a=Help then
	print(`Help(): The primary help procedure.`):
elif a=Help1 then
	print(`Help1(): The secondary help procedure.`):
elif a=Fk then 
	print(`Fk(k,t): the g.f. whose t^n coeff. is the number of `):
	print(`(n+1,n+2)-core partitions with parts repeated <=k times`):
	print(EG):
	print(`Fk(2,t);`):
elif a=RepParts then
	print(`RepParts(a,b,k): All (a,b)-core partitions where each part is`):
	print(`repeated at most k (=1 by default) times.`):
	print(EG):
	print(`RepParts(4,5,1);`):
elif a=DyckPaths then
	print(`DyckPaths(r,s,a,b): set of all walks from`):
	print(`(0,0) to (a,b) staying >= y=x*s/r.`):
	print(EG):
	print(`DyckPaths(4,5,4,5);`):
elif a=Path2Part then
	print(`Path2Part(W): converts a walk to (a,b), where`):
	print(`a,b are coprime, to an (a,b)-core partition.`):
	print(EG):
	print(`Path2Part([[0,0],[0,1],[1,1],[1,2],[2,2],[2,3],[2,4],[3,4],[4,4],[4,5],[5,5],[6,5],[7,5],[8,5]]);`):
elif a=GoodWalks then
	print(`GoodWalks(a,b): All (a,b) Dyck paths mapping to partitions with distinct parts.`):
	print(EG):
	print(`seq(nops(GoodWalks(n,n+1)),n=1..5);`):
elif a=GoodParts then
	print(`GoodParts(a,b): All (a,b)-core partitions with distinct parts.`):
	print(EG):
	print(`GoodParts(4,5);`):
elif a=NewGoodParts then 
	print(`NewGoodParts(a,b): A faster way to list all (a,b)-core partitions with distinct parts.`):
	print(EG):
	print(`NewGoodParts(4,5);`):
elif a=NewGoodWalks then
	print(`NewGoodWalks(a,b): A faster way to find all (a,b) Dyck paths mapping to partitions with distinct parts.`):
	print(EG):
	print(`NewGoodWalks(4,5);`):
elif a=GF1 then
	print(`GF1(a,b,t): sum of t^size(p) for p ranging over`):
	print(`(a,b)-core partitions with distinct parts.`):
	print(EG):
	print(`GF1(4,5,t);`):
elif a=GFgw then
	print(` GFgw(a,b,t): Like GF1(a,b,t): only faster.`):
	print(EG):
	print(`GFgw(4,5,t);`):
elif a=GFgwN then
	print(` GFgwN(a,b): Enumeration version of GFgw.`):
	print(EG):
	print(`GFgwN(4,5);`):

elif a=GFNK then
	print(`GFNK(N,K,t): The sum of the Nth factorial moment of GF1(k,k+1,t) for k=0..K.`):
	print(EG):
	print(`GFNK(0,10,t);`):
elif a=Poly2GF then
	print(`Poly2GF(p,t): uses built-in guessgf to guess a`):
	print(`generating function whose series begins like p(t).`):
	print(EG):
	print(`Poly2GF(GFNK(0,10,t),t);`):
elif a=F then
	print(`F(n): the nth Fibonacci number.`):
	print(EG):
	print(`F(0);`):
elif a=GuessFibPol1 then
	print(`GuessFibPol1(L,n,d): Inputs list L, outputs [A,B], where A,B are polys`):
	print(`of degree <=d s.t. L[n]=F(n)A(n)+F(n+1)B(n).`):
	print(EG):
	print(`GuessFibPol1([1,1,2,3,5,8,13],n,1);`):
elif a=GuessFibPol then
	print(`GuessFibPol(L,n): Inputs list L, outputs [A,B], where A,B are polys,`):
	print(`s.t. L[n]=F(n)A(n)+F(n+1)B(n).`):
	print(EG):
	print(`GuessFibPol([1,1,2,3,5,8,13],n);`):
elif a=GuessFibPolList then
	print(`A list of [A,B] pairs where the kth pair interpolates`):
	print(`kth column of L. (See GuessFibPol.)`):
	print(EG):
	print(`GuessFibPolList(gwMomList(10,1),n);`):
elif a=Moms then
	print(`Moms(f,t,K,div:=0): straight 1-Kth moments of f.`):
	print(`if div=1, divides by subs(t=1,f).`):
	print(EG):
	print(`Moms(GFgw(4,5,t),t,3);`):
elif a=gwMomList then
	print(`gwMomList(N,K): outputs NxK list s.t. L[n][k] is the kth straight`):
	print(`(and before division) moment of GFgw(n,n+1,t).`):
elif a=gwStraightMoms then 
	print(`gwStraightMoms(N,K,n,r): list of K expressions in n and r describing the`):
	print(`straight (after dividing) moments of GFgw(n,n+1,t).  r denotes F(n+1)/F(n).`):
	print(EG):
	print(`gwStraightMoms(13,2,n,r);`):
elif a=Straight2Central then 
	print(`Straight2Central(L): converts list of straight moments (starting with first)`):
	print(`to central ones.`):
	print(EG):
	print(`Straight2Central(gwStraightMoms(13,2,n,r));`):
elif a=CoreThms then
	print(`CoreThms(N,K,story:=1,longprint:=0): using data from Gs(n,n+1,t), n=1..N, computes the limits of`):
	print(`standard central moments 1..K of X_s as n->infinity.`):
	print(`If story=1, outputs theorems; otherwise returns a list.  Use longprint=1 to`):
	print(`print expressions in Maple-friendly format.`):
	print(EG):
	print(`CoreThms(17,3);`):
elif a=Gs then
	print(`Gs(q,s): Sum of q^|p|, where p ranges over (s,s+1)-core`):
	print(`partitions with distinct parts.`):
	print(EG):
	print(`Gs(q,4)`):
elif a=qGuessRec1 then
	print(`qGuessRec1(L,q,n,k,d): inputs list L of polys in q.`):
	print(`Outputs [ini,rec] s.t. ini=L[1..k] and`):
	print(`rec is a list of k+1 polys in q^n of degree <=d`):
	print(`satisfying rec[1](n)*L[n-k]+...+rec[k+1](n)*L[n]=0`):
	print(`for all n.  If not possible, returns FAIL.`):
	print(EG):
	print(`qGuessRec1([seq(q^i,i=1..50)],q,n,1,0);`):
elif a=qGuessRec then 
	print(`qGuessRec(L,q,n,dmax:=4): inputs list L of polys in q.`):
	print(`Outputs [ini,rec] s.t. ini=L[1..k] and`):
	print(`rec is a list of k+1 polys in q^n of degree <=dmax`):
	print(`satisfying rec[1](n)*L[n-k]+...+rec[k+1](n)*L[n]=0`):
	print(`for all n.  If not possible, returns FAIL.`):
	print(EG):
	print(`qGuessRec([seq(q^i,i=1..50)],q,n);`): 
elif a=qRec2Seq then
	print(`qRec2Seq(Rec,q,n,N): inputs Rec=[ini,rec] specifying`):
	print(`a q-recurrence and outputs the first N terms.`):
	print(EG):
	print(`qRec2Seq([[q],[q,-1]],q,n,5);`):
elif a=Gsg then 
	print(`Gsg(q,s,d): sum of q^size(P), where P ranges over`):
	print(`partitions of perimeter<=s with parts distance >=d apart.`):
	print(EG):
	print(`Gsg(q,3,1)`):
elif a=Gkl then
	print(`Gkl(q,k,l): sum of q^size(P), where P ranges over`):
	print(`distinct part partitions with l parts and largest part k.`):
	print(EG):
	print(`Gkl(q,3,2)`):
elif a=Gklg then
	print(`Gklg(q,k,l,d): sum of q^size(P), where P ranges over`):
	print(`partitions with l parts distance >=d from each other and largest part k.`):
	print(EG):
	print(`Gklg(q,3,2,1)`):
elif a=GuessPol1 then
	print(`GuessPol1(X,Y,x,Deg): tries to find the UNIQUE degree Deg polynomial`):
	print(`in x that interpolates the lists the lists X,Y.  Returns FAIL if not possible.`):
elif a=GuessPol then
	print(`GuessPol(L,n,K:=3,MaxDeg:=10,MinDeg:=0): finds the polynomial in n with smallest degree between`):
	print(`MinDeg and MaxDeg which interpolates L. For safety, if nops(L)<=degree+K, returns FAIL.`):
	print(EG):
	print(`GuessPol([seq(n^2,n=1..20)],n);`):
elif a=qGuessRecOpt then
	print(`qGuessRecOpt(L,q,n,R:=4): like qGuessRec, but uses a different`):
	print(`search pattern, letting k+d=1..R where k is the order and d is the`):
	print(`max degree of the polynomials.`):
	print(EG):
	print(`qGuessRecOpt([seq(q^i,i=1..50)],q,n,2);`):
elif a=qBin then
	print(`qBin(q,m,n): the q-binomial coefficient {m \choose n}_q.`):
	print(EG):
	print(`qBin(q,3,2);`):
elif a=GsTheo then
	print( `GsTheo(q,s): computes Gs(q,s) using the formula in Theorem 2.2.`):
	print(EG):
	print(`Gs(q,5)-GsTheo(q,5);`):
elif a=GsRec then
	print(`GsRec(q,N): computes [Gs(q,1),..Gs(q,N)] using the recursion.`):
	print(EG):
	print(`GsRec(q,6)-[seq(Gs(q,s),s=1..6)];`):
elif a=GuessGFgwN then
	print(`GuessGFgwN(a,b,n,MaxD:=10,K:=20): inputs expressions a,b`):
	print(`depending on n, guesses polynomial for GFgwN(a,b)`):
	print(`of max degree MaxD.  Uses K data points.`):
	print(EG):
	print(`GuessGFgwN(3,3*n+1,n);`):
elif a=Nds then
	print(`Nds(d,s): The number of (s,ds-1)-core partitions with distinct parts.`):
	print(EG):
	print(`Nds(3,7);`):
elif a=Gds then
	print(`Gds(q,d,s): The g.f. (according to size) of (s,ds-1)-core partitions with distinct parts.`):
	print(EG):
	print(`Gds(q,3,7);`):
elif a=abMoms then 
	print(`abMoms(a,b,s,N,K): inputs numeric a,b with gcd(a,b)=1; symbolic s, large integer N; `):
	print(`outputs number, expectation, and up to Kth central moments of`):
	print(`size of a,as+b - core partitions as functions of s. Also outputs limits of standardized moments.`):
	print(EG):
	print(`abMoms(4,1,s,20,3);`):
elif a=OddParts then
	print(`OddParts(a,b): the set of (a,b) core partitions into odd parts.`):
	print(EG):
	print(`OddParts(5,6);`):
elif a=Idealize then 
	print(`Idealize(p): converts a partition (set of partitions)`):
	print(`to an order ideal (set of ideals).`):
	print(EG):
	print(`Idealize([3,2,1]); `):
else
	print(`No such procedure exists.`):
fi:
end:


#############################################
#LISTING (a,b) DYCK PATHS AND CORE PARTITIONS
#############################################

DyckPaths:=proc(r,s,a,b) local S,S1,p: option remember:
if a=0 and b=0 then return {[[0,0]]}: fi:

if a<0 or b<0 or b*r<a*s then return {}: fi:

S1:=DyckPaths(r,s,a-1,b):
S:={seq([op(p),[a,b]],p in S1)}:

if (b-1)*r>=a*s then
	S1:=DyckPaths(r,s,a,b-1):
	S:=S union {seq([op(p),[a,b]],p in S1)}:
fi:
S:
end:


Path2Part:=proc(W) local S,row,i,j,a,b:
a:=W[-1][1]:
b:=W[-1][2]:
S:={}:
row:=0:
for i from 2 to nops(W) do
   if W[i]-W[i-1]=[0,1] then
      S:=S union 
         {seq(row*a-j*b,j=W[i][1]+1..floor(row*a/b))}:
      row:=row+1:
      if row=b then break: fi:
   fi:
od:
[seq(S[-i]-nops(S)+i,i=1..nops(S))]:
end:


GoodWalks:=proc(a,b) local S,S1,w,P:
S:=DyckPaths(a,b,a,b):
S1:={}:
for w in S do
   P:=Path2Part(w):
   if nops({op(P)})=nops(P) then S1:=S1 union {w}: fi:
od:
S1:
end:


GoodParts:=proc(a,b) local S,S1,w,P:
S:=DyckPaths(a,b,a,b):
S1:={}:
for w in S do
   P:=Path2Part(w):
   if nops({op(P)})=nops(P) then S1:=S1 union {P}: fi:
od:
S1:
end:


OddParts:=proc(a,b) local S,S1,w,P,p:
S:=DyckPaths(a,b,a,b):
S1:={}:
for w in S do
   P:=Path2Part(w):
   if `and`(seq(p mod 2=1,p in P)) then 
     S1:=S1 union {P}:
   fi:
od:
S1:
end:

Idealize:=proc(p) local i,p1:
if type(p,`list`) then
	{seq(p[i]+nops(p)-i,i=1..nops(p))}:
else
	{seq(Idealize(p1),p1 in p)}:
fi:
end:


RepParts:=proc(a,b,k:=1) local S,S1,w,P,M,m,isGood:
S:=DyckPaths(a,b,a,b):
S1:={}:
for w in S do
   P:=Path2Part(w):
   M:=convert(P,multiset):
   isGood:=true:
   for m in M do
      if m[2]>k then isGood:=false: break: fi:
   od:
   if isGood then S1:=S1 union {P}: fi:
od:
S1:
end:



NewGoodParts:=proc(a,b) local w:
{seq(Path2Part(w),w in NewGoodWalks(a,b))}:
end:

Label:=proc(r,s,a,b)
(b-1)*r-(a+1)*s:
end:



Labels:=proc(r,s,a,b) local j:
 {seq((b-1)*r-j*s,j=a+1..floor((b-1)*r/s))}:
end:

SetPlus:=proc(S,t) local s: 
{seq(s+t,s in S)}:
end:


NewGoodWalks1:=proc(r,s,a,b,S,checked) local L,L1,W,W1,w: option remember:
if a=0 and b=0 then return {[[0,0]]}: fi:
if a<0 or b<0 or b*r<a*s then return {}: fi:

L1:=Labels(r,s,a,b):

if checked then
  L:={Label(r,s,a,b)}:
else
  L:=L1:
fi:

if (SetPlus(L,1) union SetPlus(L,-1)) intersect S <> {} then return {}: fi:

W:=NewGoodWalks1(r,s,a-1,b,S,true):
W:={seq([op(w),[a,b]],w in W)}:
W1:=NewGoodWalks1(r,s,a,b-1,S union L1,false):
W union {seq([op(w),[a,b]],w in W1)}:
end:


NewGoodWalks:=proc(a,b)
NewGoodWalks1(a,b,a,b,{},true):
end:






#####################
#GENERATING FUNCTIONS
#####################

#(n+1,n+2), <=k-repetition of parts (at most k consecutive labels)


Ids:=proc(k) local ab,m,S: option remember:
ab:=(a,b)->{seq(j,j=a..b)}:
{seq(seq([m,S union ab(2,m)], S in choose(ab(max(2,m+1),k))), m=0..k)}:
end:

AmSk:=proc(m,S,k) local a,S1,ab,ids,x: option remember:
if k=0 then 
	if m=0 and S={} then return 1: 
	else return 0: fi:
fi:
ids:=Ids(k-1):
a:=0:
for x in ids do
	S1:=SetPlus(x[2],1):
	if x[1]>=m-1 and evalb(x[1]>0)=evalb(member(2,S)) and S1=S minus {2} then
		a:=a+AmSk(op(x),k-1):
	fi:
od:
a:
end:

(* 
#SLOW version of Fk
Fk:=proc(k,t) local ids,f,eq,var,ab,n,j,id,x,y,rhs,S1,S2:
ab:=(a,b)->{seq(j,j=a..b)}:
ids:={seq(seq([m,S union ab(2,m)], S in choose(ab(max(2,m+1),k))), m=0..k)}:
var:={seq(f[id], id in ids)}:
eq:={}:
for x in ids do
	rhs:=AmSk(op(x),k)*t^k:
	S1:=x[2]:
	for y in ids do
		S2:=SetPlus(y[2],1):
		if y[1]>=x[1]-1 and  (y[1]=0 or member(2,S1) or k=1) and (evalb(y[1]>0)=evalb(member(2,S1)) or k=1) and S2 minus {k+1}=S1 minus {2} and not (x[1]>0 and S1 union S2=ab(2,k+1)) and not (k=1 and x[1]>0 and y[1]>0) then
			rhs:=rhs+t*f[y]:
		fi:
	od:
	eq:=eq union {f[x]=rhs}:
od:
var:=solve(eq,var):
normal(convert(taylor(-t+1/2-(1/2)*sqrt(-4*t+1),t=0,k+2),polynom)/t^2 + subs(var, add(f[x],x in ids))):
end:

*)

#FAST version of Fk
Fk:=proc(k,t) local g:
g:=1+expand(t*convert(taylor(-t+1/2-(1/2)*sqrt(-4*t+1),t=0,k+2),polynom)/t^2):
normal(g/(1-g*t));
end:




########



GF1:=proc(a,b,t) local p:
if gcd(a,b)<>1 then return FAIL: fi:
add(t^convert(p,`+`),p in NewGoodParts(a,b)):
end:


FactorialMoment:=proc(p,t,n) option remember:
if n=0 then return subs(t=1,p): fi:
FactorialMoment(diff(p,t),t,n-1):
end:


GFNK:=proc(N,K,t) local k:
add(t^k*FactorialMoment(GF1(k,k+1,t),t,N),k=0..K):
end:


Poly2GF:=proc(p,t) local d:
guessgf([seq(coeff(p,t,d),d=0..degree(p,t))],t):
end:


GFgw1:=proc(r,s,a,b,S,checked,t,x) local L,L1,W,W1: option remember:
if a=0 and b=0 then return(1): fi:
if a<0 or b<0 or b*r<a*s then return(0): fi:

L1:=Labels(r,s,a,b):

if checked then
  L:={Label(r,s,a,b)}:
else
  L:=L1:
fi:

if (SetPlus(L,1) union SetPlus(L,-1)) intersect S <> {} then return(0): fi:

W:=GFgw1(r,s,a-1,b,S,true,t,x):
#W:={seq([op(w),[a,b]],w in W)}:
W1:=GFgw1(r,s,a,b-1,S union L1,false,t,x):
W1:=W1*x^nops(L1)*t^convert(L1, `+`):
W+W1:
end:



#GFgw(a,b,t): an efficient version of GF1
GFgw:=proc(a,b,t) local x,f,i: 
if gcd(a,b)<>1 then return FAIL: fi:
f:=GFgw1(a,b,a,b,{},true,t,x):

expand(add(coeff(f,x,i)/t^binomial(i,2),i=0..degree(f,x))):

end:


GFgw1N:=proc(r,s,a,b,S,checked) local L,L1,W,W1: option remember:
if a=0 and b=0 then return(1): fi:
if a<0 or b<0 or b*r<a*s then return(0): fi:

L1:=Labels(r,s,a,b):

if checked then
  L:={Label(r,s,a,b)}:
else
  L:=L1:
fi:

if (SetPlus(L,1) union SetPlus(L,-1)) intersect S <> {} then return(0): fi:

GFgw1N(r,s,a-1,b,S,true)+GFgw1N(r,s,a,b-1,S union L1,false):
end:


#GFgwN(a,b,t):
GFgwN:=proc(a,b) local x,f,i: 
if gcd(a,b)<>1 then return FAIL: fi:
GFgw1N(a,b,a,b,{},true):
end:


Gkl:=proc(q,k,l) local k1: option remember:
if l=1 then return q^k: fi:
expand(q^k*add(Gkl(q,k1,l-1),k1=1..k-1)):
end:


Gs:=proc(q,s) local k,l:
1+add(add(Gkl(q,k,l),l=1..s-k),k=1..s):
end:


Gklg:=proc(q,k,l,d) local k1: option remember:
if l<0 then return 0: fi:
if l=0 then return 1: fi:
if l=1 then return q^k: fi:
expand(q^k*add(Gklg(q,k1,l-1,d),k1=1..k-d)):
end:


Gsg:=proc(q,s,d) local k,l:
1+add(add(Gklg(q,k,l,d),l=1..s-k),k=1..s):
end:


#################
#GUESSING MOMENTS
#################


F:=proc(n)  option remember:
if n=0 or n=1 then
 n:
else
 F(n-1)+F(n-2):
fi:
end:

GuessFibPol1:=proc(L,n,d) local A,B,a,b,f,i,eq,var:
if nops(L)<2*d+4 then
  return FAIL:
fi:
A:=add(a[i]*n^i,i=0..d):
B:=add(b[i]*n^i,i=0..d):
var:={seq(a[i],i=0..d),seq(b[i],i=0..d)}:
eq:={seq(L[i]=subs(n=i,A)*F(i)+subs(n=i,B)*F(i+1),i=1..2*d+4)}:
var:=solve(eq,var):
A:=subs(var,A):
B:=subs(var,B):
if var=NULL 
  or {seq(L[i]-subs(n=i,A)*F(i)-subs(n=i,B)*F(i+1),i=2*d+5..nops(L))}<>{0}
  then return FAIL: fi:

[A,B]:
end:


GuessFibPol:=proc(L,n) local d,AB:
for d from 0 to nops(L)/2-2 do
  AB:=GuessFibPol1(L,n,d):
  if AB<>FAIL then return AB: fi:
od:
FAIL:
end:


Moms:=proc(f,t,K,div:=0) local L,tD,g,k:
if K=0 then return [subs(t=1,f)]: fi:
tD:=g->t*diff(g,t):
g:=f:
L:=[]:
for k from 1 to K do
  g:=tD(g):
  L:=[op(L),subs(t=1,g)]:
od:

if div=1 then return L/subs(t=1,f): fi:
L:
end:


GuessFibPolList:=proc(L,n) local k,nn,LL,AB:
LL:=[]:
for k from 1 to nops(L[1]) do
  AB:=GuessFibPol([seq(L[nn][k],nn=1..nops(L))],n):
  if AB=FAIL then return LL: fi:
  LL:=[op(LL),AB]:
od:
LL:
end:


gwMomList:=proc(N,K) local q,n,L,t: option remember:
#[seq(Moms(Gs(t,n),t,K),n=1..N)]:
L:=qRec2Seq([[1,1+q,q^2+q+1,2*q^3+q^2+q+1],[-q^n/q,-q^n/q,0, -1, 1]],q,n,N):
[seq(Moms(L[n],q,K),n=1..N)]
end:


gwStraightMoms:=proc(N,K,n,r) local L,k:
L:=GuessFibPolList(gwMomList(N,K),n):
[seq(L[k][1]/r+L[k][2],k=1..nops(L))]:
end:


Straight2Central:=proc(L) local k,j:
[seq(normal(add(L[k-j]*L[1]^j*(-1)^j*binomial(k,j),j=0..k-1)+(-L[1])^k),k=1..nops(L))]:
end:


CheckNormal:=proc(L) local i,j:
for i from 1 by 2 to nops(L) do
	if L[i]<>0 then
		return false:
	fi:
od:

for i from 2 by 2 to nops(L) do
	if L[i]<>mul(2*j+1,j=1..i/2-1) then
		return false:
	fi:
od:
true:
end:


CoreThms:=proc(N,K,story:=1,longprint:=0) local gr,L,L1,L2,n,r,k,F,st,Print:
st:=time():
gr:=(1+sqrt(5))/2:
L:=gwStraightMoms(N,K,n,r):

if K<2 then 
  print(`ERROR in CoreThms: need K>=2.`):
  return FAIL:
fi:

if nops(L)<K then 
	print(`ERROR in CoreThms: increase N.`):
	return FAIL: 
fi:



L1:=Straight2Central(L):
L2:=[seq(limit(subs(r=gr,L1[k]/L1[2]^(k/2)),n=infinity),k=1..nops(L1))]:

if story=1 then
  if longprint=0 then Print:=s->print(s):
  else Print:=s->lprint(s):
  fi:

  print(`Let X_n be the random variable "size of an (n,n+1)-core partition`):
  print(`into distinct parts," and F(n) be the nth Fibonacci number.`):
  print(``):
  print(`Theorem 0: The number of (n,n+1) core partitions into distinct parts`):
  print(`Is F(n+1).`):
  print(``):
  
  print(`Theorem 1:`):
  print(`E[X_n]=`):
  Print(normal(subs(r=F(n+1)/F(n),L[1]))):
  print(``):
  print(`which is asymptotic to`):
  Print(normal(L[1])):
  print(``):
  print(`where r is the Golden ratio.`):
  print(``): print(``):
  print(`Theorem 2:`):
  print(`Var(X_n)=`):
  Print(normal(subs(r=F(n+1)/F(n),L1[2]))):
  print(``):
  print(`which is asymptotic to`):
  Print(normal(L1[2])):
  print(``):
  if limit( normal( subs(r=gr,L1[2]) / subs(r=gr,L[1])^2) , n=infinity) = 0 then
  	print(`Thus, X_n is concentrated about the mean.`):
  else
 	print(`Thus, X_n is not concentrated about the mean.`):
  fi:

  for k from 3 to nops(L1) do
    print(``):
    print(`Theorem`,k):
    print(`Central moment`,k,`of X_n is`):
    Print(normal(subs(r=F(n+1)/F(n),L1[k]))):
    print(``):
    print(`which is asymptotic to`):
    Print(normal(L1[k])):
    print(``):
    print(`The limit of the standardized moment is`,L2[k]):
    print(``):
  od:

  print(``):
  print(`In short, the first`,nops(L2), `standardized central moments of X_n limit to`):
  print(L2):

  if CheckNormal(L2) then
	print(`And hence, X_n appears asymptotically normal.`):
  else 
	print(`And hence, X_n is not asymptotically normal.`):
  fi:
  print(``):
  print(`The computation time was`,time()-st,`seconds.`):
  return NULL:
fi:

L2:
end:


######################
#GUESSING Q RECURSIONS
######################

qGuessRec1:=proc(L,q,n,k,d) local a,c,i,j,p,rec,var,var1,N,eq:
N:=k*(d+1)+4:

if nops(L)<N then return FAIL: fi:

var:={seq(seq(c[i][j],j=0..d),i=0..k)}:
rec:=[seq(add(q^(n*j)*c[i][j],j=0..d),i=0..k)]:
eq:={seq(add(subs(n=j,rec[-i])*L[j+1-i],i=1..k+1),j=k+1..nops(L))}:

N:=min(nops(eq),N):

var1:=solve(subs(q=3,eq[1..N]),var):
if var1=NULL then return FAIL: fi:

var1:=solve(eq[1..N],var):
var:={seq(a=1,a in var)}:

rec:=subs(var,subs(var1,rec)):
if var1=NULL or rec[-1]=0  then return FAIL: fi:

for i from N+1 to nops(eq) do
  if expand(subs(var1,eq[i]))<>0 then return FAIL: fi:
od:

[L[1..k],rec]:
end:


qGuessRec:=proc(L,q,n,dmax:=4) local k,Rec,d:
for k from 1 to nops(L) do
  for d from 0 to dmax do
    Rec:=qGuessRec1(L,q,n,k,d):
    if Rec<>FAIL then return Rec: fi:
  od:
od:
FAIL:
end:


qRec2Seq:=proc(Rec,q,n,N) local L,rec,i,k:
L:=Rec[1]:
if N<=nops(L) then return L[1..N]: fi:
rec:=Rec[2]:
k:=nops(rec)-1:
for i from nops(Rec[1])+1 to N do
  L:=[op(L),
    expand(-add(subs(n=i,rec[-j-1])*L[-j],j=1..k)/subs(n=i,rec[-1])) ]:
od:
L:
end:

GsRec:=proc(q,N) local n:
qRec2Seq([[1,1+q,q^2+q+1,2*q^3+q^2+q+1],[-q^n/q,-q^n/q,0, -1, 1]],q,n,N):
end:


qGuessRecOpt:=proc(L,q,n,R:=4) local k,Rec,d,r:

for r from 1 to R do
  for k from 1 to r do
    Rec:=qGuessRec1(L,q,n,k,r-k):
    if Rec<>FAIL then return Rec: fi:
  od:
od:

FAIL:
end:



#############################
#POLYNOMIAL ANSATZ
#############################


CheckAnsatz:=proc(X,Y,f,var) local neq,j,nvar,i:

neq:=nops(X):
nvar:=nops(var):

if not type(X[1],`list`) then
	return CheckAnsatz([seq([X[i]],i=1..neq)],Y,f,[var]):
fi:

for i from 1 to neq do
	if Y[i]<>subs({seq(var[j]=X[i][j],j=1..nvar)},f) then
		return false:
	fi:
od:
true:
end:


GuessPol1:=proc(X,Y,x,Deg) local eq,var,a,P,i:
P:=add(x^i*a[i],i=0..Deg):
var:={seq(a[i],i=0..Deg)}:
eq:={seq(subs(x=X[i],P)=Y[i], i=1..Deg+1)}:

var:=solve(eq,var):
P:=subs(var,P):

if var=NULL or CheckAnsatz(X,Y,P,x)=false then 
	return FAIL: 
fi:
P:
end:


GuessPol:=proc(Y,x,K:=3,MaxDeg:=10,MinDeg:=0) local d,P,X:
X:=[seq(i,i=1..nops(Y))]:

for d from MinDeg to MaxDeg do
	if nops(X)<=d+K then
		#print(`ERROR in GuessPol: the lists X and Y are too short`):
		#print(`to make a reliable guess.  Try more data points.`):
		return FAIL:
	fi:

	P:=GuessPol1(X,Y,x,d):
	if P<>FAIL then
		return P:
	fi:
od:
FAIL:
end:



#######################
#Q-BINOMIAL COEFFICIENT
#######################

qBin:=(q,m,r)->mul(1-q^(m-i),i=0..r-1)/mul(1-q^i,i=1..r):

GsTheo:=(q,s)->expand(normal(add(q^binomial(m+1,2)*qBin(q,s-m,m),m=0..s))):


#########################################
#GUESSING NUMBER OF (a,b)-CORE PARTITIONS
#########################################

GuessGFgwN:=proc(a,b,n,MaxD:=10,K:=20) local i:
GuessPol([seq(GFgwN(subs(n=i,a),subs(n=i,b)),i=1..K)],n,3,MaxD):
end:


################
#STRAUB (s,ds-1)
################


Nds:=proc(d,s) option remember:
if s=1 then return 1: fi:
if s=2 then return d: fi:
Nds(d,s-1)+d*Nds(d,s-2):
end:


Gds1:=proc(q,t,d,s,a) local i,j: option remember:
if a=s-1 then return add(q^(s*i^2/2-s*i/2+a*i)*t^i,i=0..d-1): fi:
if a>=s then return 1: fi:
expand( Gds1(q,t,d,s,a+1)+add(q^(s*i^2/2-s*i/2+a*i)*t^i,i=1..d)*Gds1(q,t,d,s,a+2) ):
end:


Gds:=proc(q,d,s) local g,t:
g:=Gds1(q,t,d,s,1):
expand(add(q^(k*(1-k)/2)*coeff(g,t,k),k=0..degree(g,t))):
end:


############################
#(a,a*s-b), where gcd(a,b)=1
############################

Gabs:=proc(q,a,b,s):
if b=1 then return Gds(q,s,a): fi:
if gcd(a,b)<>1 then return FAIL: fi:
GFgw(a,a*s-b,q):
end:


GuessPoly1:=proc(X,Y,x,Deg) local eq,var,a,P,i:
P:=add(x^i*a[i],i=0..Deg):
var:={seq(a[i],i=0..Deg)}:
eq:={seq(subs(x=X[i],P)=Y[i], i=1..Deg+1)}:

var:=solve(eq,var):
P:=subs(var,P):

if var=NULL or CheckAnsatz(X,Y,P,x)=false then 
	return FAIL: 
fi:
P:
end:


GuessPoly:=proc(X,Y,x,MaxDeg:=10,MinDeg:=0) local d,P:
if nops(X)<>nops(Y) then
	print(`ERROR in GuessPoly: X,Y must be lists of the same length.`):
	return(FAIL):
fi:

for d from MinDeg to MaxDeg do
	if nops(X)<d+4 then
		#print(`ERROR in GuessPoly: the lists X and Y are too short`):
		#print(`to make a reliable guess.  Try more data points.`):
		return FAIL:
	fi:

	P:=GuessPoly1(X,Y,x,d):
	if P<>FAIL then
		return P:
	fi:
od:
FAIL:
end:



abMoms:=proc(a,b,s,N,K) local q,L,M0,M,M1,n,k:
L:=[seq(Gabs(q,a,b,n),n=1..N)]:
M0:=GuessPoly([seq(n,n=1..N)],subs(q=1,L),s):
if M0=FAIL then return FAIL: fi:

L:=[seq(Moms(L[n],q,K),n=1..N)]:
for k from 1 to K do
	M[k]:=GuessPoly([seq(n,n=1..N)],[seq(L[n][k],n=1..N)],s):
	if M[k]=FAIL then break: fi:
	M[k]:=normal(M[k]/M0):
od:
M:=convert(M,`list`):
M1:=Straight2Central(M):
M0,[M[1],op(2..nops(M),M1)],
[seq(limit(M1[k]/M1[2]^(k/2),s=infinity),k=2..nops(M))]:
end: