#######################################################################
##OddArmin.txt: Save this file as  OddArmin.txt                       #
## To use it, stay in the                                             #
##same directory, get into Maple (by typing: maple <Enter> )          #
##and then type:  read  OddArmin.txt<Enter>                           #
##Then follow the instructions given there                            #
##                                                                    #
##Written by Anthony Zaleski and Doron Zeilberger, Rutgers University #
#DoronZeil at gmail dot com                                           #
#######################################################################
 
#Created: Aug.-Nov. 2017
 
print(`Created:  Aug.-Nov. 2017. `):
print(` This is OddArmin.txt `):
print(`It is one of the packages that accompany the article `):
print(`On the Intriguing Problem of Counting (n+1,n+2)-core partitions into Odd Parts`):
print(`by Anthony Zaleski and Doron Zeilberger`):

print(` available from the authors' websites`):
print(``):
print(`Please report bugs to DoronZeil at gmail dot com `):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/  .`):

print(`---------------------------------------`):
 print(`For a list of the Supporting procedures type ezra1();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the Story procedures type ezraS();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):

print(`---------------------------------------`):
 print(`For a list of the Checking procedures type ezraC();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

print(`---------------------------------------`):
 print(`For a list of the MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

with(combinat):

ezraS:=proc()

if args=NULL then
 print(` The story procedures are:  OAseqS, , StoryAD, StoryAll, StoryO `): 
 print(``):

else
ezra(args):
fi:

end:

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are:  AnSG, AnSGPdb, AnSGP, AnSGPa, CtoR, Dec, GoodP1, GoodP, GuessRec, JoP, KidsProf, Label, NuOI, NuOIP,`):
 print(` OIbuW, OIgf, OIP, Prof, ProfDB, SpP111, SpP11, SpP1, SpP , Wid  `):
 print(``):

else
ezra(args):
fi:

end:

ezraC:=proc()

if args=NULL then
 print(` The Checking procedures are:  AnG, AnGr, CheckDec, CheckJoP1, CheckJoP, CheckOIP, CheckOIP1, CheckSpP, CheckSpP1, FkSlow, GkSlow  `):
else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are:    Ak, Fk, Gk, NuOIr, NuOIG, NuOIGr, NuOIGr1, OAseq, OAseqR, OAseqR1, OI, OIr, OIG   `):


elif nargs=1 and args[1]=AnG then
print(`AnG(n,c): The set of good order ideals in OI(n) with spacing c. Try: `):
print(` AnG(5,1); `):

elif nargs=1 and args[1]=AnGr then
print(`AnGr(n,r,c): The set of good order ideals in the union of OIr(n,r1) for r1 from 0 to r, with spacing c. Done by brute force.`):
print(`Try: `):
print(`AnGr(5,2,1);`):

elif nargs=1 and args[1]=AnGr1 then
print(`AnGr1(n,r,c): The set of good order ideals in the union of OIr(n,r) with spacing c. Done by brute force.`):
print(`Try: `):
print(`AnGr1(5,2,1);`):

elif nargs=1 and args[1]=AnSG then
print(`AnSG(n,c): The set of good order ideals in OI(n) with spacing c. Try: `):
print(` AnSG(5,1); `):

elif nargs=1 and args[1]=AnSGP then
print(`AnSGP(n,c,P): The set of semi-good order-ideals in OIn(n) with spacing c with profile P. Try: `):
print(` AnSGP(5,1,[[0,1]]); `):

elif nargs=1 and args[1]=AnSGPa then
print(`AnSGPa(n,c,P,a): The set of semi-good order-ideals in OIn(n) with spacing c with profile P and whose main diagonal has a elements. Try:`):
print(`AnSGPa(7,1,[[1,1]],1);`):

elif nargs=1 and args[1]=AnSGPdb then
print(`AnSGPdb(n,c,P):  a database for AnSGPdb(n,c,P), according to its a. Try:`):
print(`AnSGdb(7,1,[[0,0]]);`):


elif nargs=1 and args[1]=CheckDec then
print(`CheckDec(n): checks procedure Dec(S,n) . Try:`):
print(`CheckDec(4);`):


elif nargs=1 and args[1]=CheckJoP then
print(`CheckJoP(n,c): checks procedure JoP(a,P1,P2,c) for all members of AnSG(n,c):`):
print(` Try: `):
print(` CheckJoP(7,1); `):

elif nargs=1 and args[1]=CheckJoP1 then
print(`CheckJoP1(n,c,P,a): checks procedure JoP(a,P1,P2,c) for all the members of KidsProf(n,c,P,a).`):
print(`Try:`):
print(`CheckJoP1(7,1,[[1,1]],0);`):

elif nargs=1 and args[1]=CheckOIP then
print(`CheckOIP(n,c): checks whether OIP(n,c,P)=AnSGP(n,c,P) for all profiles P. Try:`):
print(`CheckOIP(3,1);`):

elif nargs=1 and args[1]=CheckOIP1 then
print(`CheckOIP1(n,c,P): checks whether OIP(n,c,P)=AnSGP(n,c,P). Try:`):
print(`CheckOIP1(3,1,[[1,1]]);`):

elif nargs=1 and args[1]=CheckSpP1 then
print(`CheckSpP1(n,c,P,a): checks procedure SpP(n,c,P,a) `):
print(` Try: `):
print(` CheckSpP1(7,1,[[1,1]],3); `):

elif nargs=1 and args[1]=CheckSpP then
print(`CheckSpP(n,P): Checks CheckSpP1(n,c,P,a) for all a and c. Try:`):
print(`CheckSpP(8,[[1,1]]);`):

elif nargs=1 and args[1]=CtoR then
print(`CtoR(S,t), the rational function, in t, whose coefficients`):
print(`are the members of the C-finite sequence S. For example, try:`):
print(`CtoR([[1,1],[1,1]],t);`):

elif nargs=1 and args[1]=Dec then
print(`Dec(S,n): Given an order-ideal S of P_{n+1,n+2} (i.e. a member of OI(n)), outputs the`):
print(` pair [a,I1,I2] where a is the number of continuous solid dots in the outer diagonal (0<=a<=n) `):
print(` i1 is a member of OI(a-1) (it is {} if a=0) and I2 is a member of OI(n-a-1). `):
print(` Try: `):
print(` Dec({[0,0],[1,0],[0,1]},2); `):

elif nargs=1 and args[1]=Ak then
print(`Ak(k,t): the generating function for NuOIr(n,k) for a fixed k. Try:`):
print(`Ak(3,t);`):

elif nargs=1 and args[1]=Fk then
print(`Fk(k,t): inputs a pos. integer k and a variable name t, outputs the generating function`):
print(` whose Maclaurin coeff. of t^n is the number of (n+1,n+2)-core partitions, where`):
print(`no part is repeated more than k times. Try:`):
print(`Fk(3,t);`):


elif nargs=1 and args[1]=FkSlow then
print(`FkSlow(k,t): Like Fk(k,t), but much slower. For checking only.  Try:`):
print(`FkSlow(3,t);`):

elif nargs=1 and args[1]=GkSlow then
print(`GkSlow(r,t,N): the generating function for odd (n+1,n+2)-core partitions whose corresponding order-ideals are supported `):
print(`in the outer-most k diagonals. BY GUESSING. Try:`):
print(`GkSlow(2,t,30);`):

elif nargs=1 and args[1]=Gk then
print(`Gk(r,t): the generating function, in the variable t, for odd (n+1,n+2)-core partitions whose corresponding order-ideals are supported `):
print(`in the outer-most r diagonals. BY GUESSING. Try:`):
print(`Gk(2,t);`):

elif nargs=1 and args[1]=GoodP then
print(`GoodP(r): all the good profiles of length <=r. Try:`):
print(`GoodP(4);`):

elif nargs=1 and args[1]=GoodP1 then
print(`GoodP1(r): all the good profiles of length r. Try:`):
print(`GoodP1(4);`):

elif nargs=1 and args[1]=GuessRec then
print(`GuessRec(L): inputs a sequence L and tries to guess`):
print(`a recurrence operator with constant cofficients `):
print(`satisfying it. It returns the initial  values and the operator`):
print(` as a list. For example try:`):
print(`GuessRec([1,1,1,1,1,1]);`):

elif nargs=1 and args[1]=JoP then
print(`JoP(a,P1,P2,c): inputs a non-negative integer a and two profiles P1 and P2, and labelling paramter c (0 or 1)`):
print(`Finds the profile of the parent. Try`):
print(`JoP(3,[[0,1]],[[0,1],[1,0]],0);`):

elif nargs=1 and args[1]=KidsProf then
print(`KidsProf(n,c,P,a): all the children of the  profile P for AnSGPa(n,c,P,a). Try:`):
print(`KidsProf(7,1,[[0,0]],1);`):

elif nargs=1 and args[1]=Label then
print(`Label(pt,a,c): The label of the lattice point pt in OI(n) with spacing c. Try:`):
print(`Label([0,0],5,c);`):


elif nargs=1 and args[1]=NuOI then
print(`NuOI(n): The number of  order-ideals of P_{n+1,n+2}, should be the Catalan number.`):
print(`NuOI(3);`):

elif nargs=1 and args[1]=NuOIr then
print(`NuOIr(n,r): The number of order-ideals of P_{n+1,n+2} of width <=r. Should be equal to`):
print(`nops(OIr(n,r)), but, of course, much faster.`):
print(`Try: NuOIr(3,2); `):

elif nargs=1 and args[1]=NuOIG then
print(`NuOIG(n,c): The number of good order ideals in OI(n) with spacing c, done recursively, via dynamical programming. `):
print(`and a scheme. `):
print(` Try: `):
print(` NuOIG(5,1); `):

elif nargs=1 and args[1]=NuOIGr then
print(`NuOIGr(n,r,c): The number of good order ideals in OI(n) with spacing c, and width<=r. done recursively, via dynamical programming. `):
print(`and a scheme. `):
print(` Try: `):
print(` NuOIGr(5,1,1); `):

elif nargs=1 and args[1]=NuOIGr1 then
print(`NuOIGr1(n,r,c): The number of good order ideals in OI(n) with spacing c, and width=r, done recursively, via dynamical programming. `):
print(`and a scheme. `):
print(` Try: `):
print(` NuOIGr1(5,1,1); `):

elif nargs=1 and args[1]=NuOIP then
print(`NuOIP(n,c,P): The number of semi-good order-ideals of P_{n+1,n+2} with labelling parameter c, and profile P.`):
print(`Try:`):
print(`NuOIP(3,1,[[1,1]]);`):


elif nargs=1 and args[1]=OAseq then
print(`OAseq(N): the first N terms of the two sequences related to the Odd-Armin sequence and  OEIS sequence A047749, and their`):
print(`interleavings. The output is  a list of these four sequences`):
print(`Try: `):
print(` OAseq(11); `):

elif nargs=1 and args[1]=OAseqR then
print(`OAseqR(N,r): The first N terms of two sequences related to the restricted Odd-Armin sequence whose support is the r outmost diagonals`):
print(`and their `):
print(` interleavings. The output is  a list of these four sequences `):
print(` Try: `):
print(` OAseqR(11,2); `):

elif nargs=1 and args[1]=OAseqR1 then
print(`OAseqR1(N,r): The first N terms of two sequences related to the restricted Odd-Armin sequence whose width is exactly r. `):
print(`and their `):
print(` interleavings. The output is  a list of these four sequences `):
print(` Try: `):
print(` OAseqR1(11,2); `):


elif nargs=1 and args[1]=OAseqS then
print(`OAseqS(N): the story of the two sequences related to the Odd-Armin sequence and  OEIS sequence A047749. Try:`):
print(`OAseqS(11);`):

elif nargs=1 and args[1]=OI then
print(`#OI(n): The set of order-ideals of P_{n+1,n+2}. Try:`):
print(`OI(3);`):

elif nargs=1 and args[1]=OIbuW then
print(`OIbuW(n): inputs a positive integer n and outputs a list of length n+1 whose i+1-th entry is the set of`):
print(`order ideals with width i. Try: `):
print(`OIbuW(n): `):

elif nargs=1 and args[1]=OIr then
print(`OIr(n,r): The set of order-ideals of P_{n+1,n+2} with width <=r, done recursively. `):
print(`OIr(3,1);`):


elif nargs=1 and args[1]=OIG then
print(`OIG(n,c): The set of good order ideals in OI(n) with spacing c, done recursively, via dynamical programming. `):
print(`and a scheme. `):
print(` Try: `):
print(` OIG(5,1); `):

elif nargs=1 and args[1]=OIP then
print(`OIP(n,c,P): The subset of OI(n) that are semi-good and have profile P with labelling parameter c. Try:`):
print(`OIP(3,1,[[1,1]]);`):

elif nargs=1 and args[1]=OIgf then
print(`OIgf(t): the generating function for NuOI(n)`):

elif nargs=1 and args[1]=Prof then
print(`Prof(S,n,c):The profile of the order ideal S with labelling parameter c. It is only useful for members of AnG(n,c)`):
print(`and AnSG(n,c) but makes sense for any order-ideal. Try:`):
print(`Prof({[2,0]},3,1);`):

elif nargs=1 and args[1]=ProfDB then
print(`ProfDB(n,c): Inputs a positive integer n. Creates a database, according to profiles of the members of AnSG(n,c):`):
print(`ProfDB(5,1);`):

elif nargs=1 and args[1]=SpP111 then
print(`SpP111(Z): Inputs a pair [a,b] and outputs the set {[[a,0],[1,b]],[[a,1],[0,b]]`):


elif nargs=1 and args[1]=SpP11 then
print(`SpP11(L): Inputs a lists of pairs L and outputs the set of pairs of lists of the same length that gives it. Try:`):
print(`SpP11([[1,1],[0,0]])`);

elif nargs=1 and args[1]=SpP1 then
print(`SpP1(L,r,s): Given a list of pairs, finds all their splittings [L1,L2] with nops(L1)=r and nops(L2)=s  Try:`):
print(`SpP1([[1,1],[0,0],[0,1]],1,3):`):

elif nargs=1 and args[1]=SpP then
print(` SpP(n,c,P,a): Inputs a profile in OI(n) and positive integers n and a, (0<=a<=n) and labelling parameter c`):
print(`(0 or 1). Outputs`):
print(`the set of pairs {[P1,P2]} where P1 is a profile of OI(a-1) and P2 is a profile of OI(n-a-1).`):
print(`Try: `):
print(`SpP(8,1,[[1,1]],3); `):


elif nargs=1 and args[1]=StoryAD then
print(`StoryAD(R,t,N): tells the story of  the generating functions enumerating (n+1,n+2)-core partitions where each part`):
print(`is repeated at most  k times for k from 1 to R `):
print(`N is a pos. integer indicating how many terms in the sequence to display.`):
print(` Try:`):
print(`StoryAD(5,t,30);`):

elif nargs=1 and args[1]=StoryAll then
print(`StoryAll(R,t,N): tells the story of enumerating order ideals of P_{n+1,n+2}  supported in the first k outermost diagonals`):
print(`for k from 1 to R `):
print(`it out puts the sequences up to the N-th term, Try:`):
print(`StoryAll(5,t,30);`):

elif nargs=1 and args[1]=StoryO then
print(`StoryO(N,R,t): tells the story of restriced odd (n+1,n+2)-core partitions up to width R, computing N terms. Try: `):
print(`StoryO(30,2,t);`):


elif nargs=1 and args[1]=Wid then
print(`Wid(S,n): Given an order ideal, S, or  P_{n+1,n+2} (i.e. a member of OI(n)) outputs its width. Try`):
print(`Wid({[3,0]},4);`):

else
print(`There is no ezra for`,args):
fi:

end:


##From Maple package Cfinite

#SeqFromRec(S,N): Inputs S=[INI,ope]
#where INI is the list of initial conditions, a ope a list of
#size L, say, and a recurrence operator ope, codes a list of
#size L, finds the first N0 terms of the sequence satisfying
#the recurrence f(n)=ope[1]*f(n-1)+...ope[L]*f(n-L).
#For example, for the first 20 Fibonacci numbers,
#try: SeqFromRec([[1,1],[1,1]],20);
SeqFromRec:=proc(S,N) local gu,L,n,i,INI,ope:
INI:=S[1]:ope:=S[2]:
if not type(INI,list) or not type(ope,list) then
 print(`The first two arguments must be lists `):
 RETURN(FAIL):
fi:
L:=nops(INI):
if nops(ope)<>L then
 print(`The first two arguments must be lists of the same size`):
RETURN(FAIL):
fi:

if not type(N,integer) then
 print(`The third argument must be an integer`, L):
 RETURN(FAIL):
fi:

if N<L then
 RETURN([op(1..N,INI)]):
fi:

gu:=INI:

for n from 1 to N-L do
 gu:=[op(gu),expand(add(gu[nops(gu)-i+1]*ope[i],i=1..L))]:
od:

gu:
end:

#CtoR(S,t), the rational function, in t, whose coefficients
#are the members of the C-finite sequence S. For example, try:
#CtoR([[1,1],[1,1]],t);

CtoR:=proc(S,t) local D1,i,N1,L1,f,f1,L:
if not (type(S,list) and  nops(S)=2 and type(S[1],list) and type(S[2],list)
        and nops(S[1])=nops(S[2]) and type(t, symbol) ) then
   print(`Bad input`):
   RETURN(FAIL):
fi:

D1:=1-add(S[2][i]*t^i,i=1..nops(S[2])):
N1:=add(S[1][i]*t^(i-1),i=1..nops(S[1])):
L1:=expand(D1*N1):
L1:=add(coeff(L1,t,i)*t^i,i=0..nops(S[1])-1):
f:=L1/D1:
L:=degree(D1,t)+10:
f1:=taylor(f,t=0,L+1):
if expand([seq(coeff(f1,t,i),i=0..L)])<>expand(SeqFromRec(S,L+1)) then
print([seq(coeff(f1,t,i),i=0..L)],SeqFromRec(S,L+1)):
 RETURN(FAIL):
else
 RETURN(factor(f)):
fi:
end:

#GuessRec1(L,d): inputs a sequence L and tries to guess
#a recurrence operator with constant cofficients of order d
#satisfying it. It returns the initial d values and the operator
# as a list. For example try:
#GuessRec1([1,1,1,1,1,1],1);
GuessRec1:=proc(L,d) local eq,var,a,i,n:

if nops(L)<=2*d+2 then
 print(`The list must be of size >=`, 2*d+3 ):
 RETURN(FAIL):
fi:

var:={seq(a[i],i=1..d)}:

eq:={seq(L[n]-add(a[i]*L[n-i],i=1..d),n=d+1..nops(L))}:

var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
else
 RETURN([[op(1..d,L)],[seq(subs(var,a[i]),i=1..d)]]):
fi:

end:




#GuessRec(L): inputs a sequence L and tries to guess
#a recurrence operator with constant cofficients 
#satisfying it. It returns the initial  values and the operator
# as a list. For example try:
#GuessRec([1,1,1,1,1,1]);
GuessRec:=proc(L) local gu,d:

for d from 1 to trunc(nops(L)/2)-2 do
 gu:=GuessRec1(L,d):
 if gu<>FAIL then
   RETURN(gu):
 fi:
od:
FAIL:

end:


#End  Maple package Cfinite




#IsSGood(S,n,c): Is the member S of OI(n) , with spacing c semi-good? . Try
#S:=OI(6)[100]; IsSGood(S,6,1);
IsSGood:=proc(S,n,c) local gu,gu1, i,j:

if S={} then
 RETURN(true):
fi:

gu:=SigL(S,n,c):


for i from 1 to nops(gu) do
 gu1:=gu[i]:
 
 for j from 1 to nops(gu1)-1 do
 
 if gu1[j]=gu1[j+1] then
  RETURN(false):
 fi:
 
 od:

od:


true:

end:

#SigL(S,n,c): The signature of the order ideal S,  as a list of lists,  in OI(n) with spacing parity c. Try
#S:=OI(5)[100]; SigL(S,5,0);
SigL:=proc(S,n,c) local gu,i,j:
gu:=StoLL(S):
gu:=[seq([seq (Label(gu[i][j],n,c) mod 2, j=1..nops(gu[i]))],i=1..nops(gu))]:
end:


#Label(pt,a,c): The label of the lattice point pt in Tr(a,a) or Tr(a,b) for any b. Try
#Label([0,0],5,c);
Label:=proc(pt,a,c):
1+(a-1)*c-pt[1]*c-pt[2]*(c-1) mod 2 :
end:

#StoL1(S): inputs a set of 2D lattice points with the same sum-of-coordinates, outputs a list where they
#are ordered in increasing order of their second coodinate. Try:
#StoL1({[5,0],[4,1],[3,2]});
StoL1:=proc(S) local i,pt,lu,m:

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

lu:={seq(pt[1]+pt[2],pt in S)}:

if nops(lu)<>1 then
 RETURN(FAIL):
fi:

m:=lu[1]:

lu:={seq(pt[2],pt in S)}:

lu:=sort(convert(lu,list)):

[seq([m-lu[i],lu[i]],i=1..nops(lu))]:

end:

#StoLL(S): inputs a list of lattice points, converts it to a list of lists where the first entry
#is list of lattice points with largest sum, the second entry is with one less, etc.
#Try
#StoLL({[7,0],[6,1],[0,7],[6,0]}):
StoLL:=proc(S) local m,pt,i,T,gu:

m:=max(seq(pt[1]+pt[2],pt in S)):

for i from 0 to m do
T[i]:={}:
od:

for pt in S do
T[pt[1]+pt[2]]:=T[pt[1]+pt[2]] union {pt}:
od:

gu:=[seq(T[i],i=0..m)]:

for i from 1 to nops(gu) while gu[i]={} do od:

gu:=[op(i..nops(gu),gu)]:

gu:=[seq(gu[nops(gu)+1-i],i=1..nops(gu))]:

[seq(StoL1(gu[i]),i=1..nops(gu))]:

end:

#CF(S): the canonical form of the order-ideal S. Try
#CF({});
#CF(OI(4)[40]);
CF:=proc(S) local gu,i:
gu:=StoLL(S):
[seq(op(gu[i]),i=1..nops(gu))]:
end:

#Sig(S,a,c): The signature of the order ideal S in OI(a) with spacing parity c. Try
#S:=OI(5)[100]; Sig(S,5,0);
Sig:=proc(S,a,c) local gu,i:
gu:=CF(S):
gu:=[seq (Label(gu[i],a,c) mod 2, i=1..nops(gu))]:
end:

#IsGood(S,a,c): Is the order-ideal S of OI(a), with spacing c good? . Try
#S:=OI(5)[100]; IsGood(S,5,1);
IsGood:=proc(S,a,c) local gu,i:

if S={} then
 RETURN(true):
fi:

gu:=Sig(S,a,c):

if gu[1]=0 then
 RETURN(false):
fi:

for i from 1 to nops(gu)-1 do
 if gu[i]=gu[i+1] then
  RETURN(false):
 fi:
od:


true:

end:


#AnG(n,c): The set of good order ideals in OI(n) with spacing c. Try
#AnG(5,1);
AnG:=proc(n,c) local gu,mu,mu1:
option remember:

mu:=OI(n):

gu:={}:

for mu1 in mu do
 if IsGood(mu1,n,c) then
    gu:=gu union {mu1}:
 fi:
od:

gu:

end:


#AnGr1(n,r,c): The set of good order ideals in OIr(n,r) with spacing c. Try
#AnGr1(5,2,1);
AnGr1:=proc(n,r,c) local gu,mu,mu1:
option remember:

mu:=OIr(n,r):

gu:={}:

for mu1 in mu do
 if IsGood(mu1,n,c) then
    gu:=gu union {mu1}:
 fi:
od:

gu:

end:

#AnGr(n,r,c): The set of good order ideals in the union of OIr(n,r) for r1<=r, with spacing c. Try
#AnGr(5,2,1);
AnGr:=proc(n,r,c) local gu,mu,mu1,r1:
option remember:

mu:=OIr(n,r):

gu:={}:

for mu1 in mu do
 if IsGood(mu1,n,c) then
    gu:=gu union {mu1}:
 fi:
od:

gu:

end:

#AnSG(n,c): The set of semi-good order ideals in OI(n) with spacing c. Try
#AnSG(5,1);
AnSG:=proc(n,c) local gu,mu,mu1:
option remember:

mu:=OI(n):

gu:={}:

for mu1 in mu do
 if IsSGood(mu1,n,c) then
    gu:=gu union {mu1}:
 fi:
od:

gu:

end:

###########################End From previous Ks

#Hag(kv,v): adds v to each point of kv
Hag:=proc(kv,v) local pt:
{seq([pt[1]+v[1],pt[2]+v[2]],pt in kv)}:
end:

#Label(pt,a,c): The label of the lattice point pt in OI(n) with spacing c. Try:
#Label([0,0],5,c);
Label:=proc(pt,a,c):
1+(a-1)*c-pt[1]*c-pt[2]*(c-1)  :
end:



#CV1(S,n): given a subest of {[n-1,0], [n-2,1], ..., [0,n-1]} outputs its characteristic vector
CV1:=proc(S,n) local i,S1,lu:
S1:={seq([n-1-i,i],i=0..n-1)}:

if S minus S1<>{} then
 RETURN(FAIL):
fi:

lu:=[]:

for i from 0 to n-1 do
if member([n-1-i,i],S) then
 lu:=[op(lu),1]:
else
 lu:=[op(lu),0]:
fi:
od:
lu:
end:


#CV2(L): inputs a 0-1 vector and outputs a vector [b0,a1,b1,a2,b2, ..., ak,b_{k+1}] where the first and last components
#are non-negative and the rest are positive such that the vector is 0^b0 1^a1 .... 1^ak 0^b_{k+1}. Try:
#CV2([0,1,1,1,0];
CV2:=proc(L) local i,n,L1,b0,a1,gu:

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

n:=nops(L):
for i from 1 to n while L[i]=0 do od:
if i=n+1 then
 RETURN([n]):
fi:

b0:=i-1:
for i from b0+1 to n while L[i]=1 do od:
if i=n+1 then
 RETURN([b0,n-b0,0]):
fi:
a1:=i-1-b0:

L1:=[op(i..n,L)]:

gu:=CV2(L1):
[b0,a1,op(gu)]:
end:



#Dec(S,n): Given an order-ideal S of P_{n+1,n+2} (i.e. a member of OI(n)), outputs the
#pair [a,I1,I2] where a is the number of continuous solid dots in the outer diagonal (0<=a<=n)
#i1 is a member of OI(a-1) (it is {} if a=0) and I2 is a member of OI(n-a-1).
#Try:
#Dec({[0,0],[1,0],[1,1]},2);
Dec:=proc(S,n) local  a,i1,j1,S1,S2,ku:
option remember:

for i1 from 1 to n+1 while member([n-i1,i1-1],S) do od:
a:=i1-1:


if a=0 then
 RETURN([0,{},Hag(S,[0,-1])]):
fi:

ku:={seq(seq([i1,j1], i1=n-a..n-1-j1 ),j1=0..a-1)  }:


S1:=S intersect ku:
S2:=S minus S1:


if {seq([n-i1,i1-1],i1=1..a)} minus S1<>{} then
   RETURN(FAIL):
fi:

S1:=S1 minus {seq([n-i1,i1-1],i1=1..a)}:
S1:=Hag(S1,[-(n-a),0]):


S2:=Hag(S2,[0,-a-1]):

[a,S1,S2]:

end:

#CheckDec(n): checks procedure Dec(S,n) . Try:
#CheckDec(4);
CheckDec:=proc(n) local gu,gu1,T1,T2,a,lu:

gu:=OI(n):

for a from 0 to n do
 T1[a]:={}:
 T2[a]:={}:
od:

for gu1 in gu do
lu:=Dec(gu1,n):
a:=lu[1]:
T1[a]:=T1[a] union {lu[2]}:
T2[a]:=T2[a] union {lu[3]}:
od:

if T1[0]<>{{}} then
 RETURN(false):
fi:

for a from 1 to n do
 if T1[a]<>OI(a-1) then
  RETURN(false):
 fi:
od:

if T2[n]<>{{}} then
  RETURN(false):
fi:

for a from 0 to n-1 do
 if T2[a]<>OI(n-1-a) then
  RETURN(false):
 fi:
od:

true:

end:


#Prof(S,n,c):The profile of the order ideal S with labelling parameter c. It is only useful for members of AnG(n,c)
#and AnSG(n,c) but makes sense for any order-ideal. Try:
#Prof({[2,0]},3,1);
Prof:=proc(S,n,c) local gu,i:
gu:=SigL(S,n,c):
[seq([gu[i][1] mod 2,gu[i][-1] mod 2], i=1..nops(gu))]:
end:


#ProfDB(n,c): Inputs a positive integer n. Creates a database, according to profiles of the members of AnSG(n,c):
#ProfDB(5,1);
ProfDB:=proc(n,c) local gu,gu1,mu,T,mu1,pu:
option remember:

 gu:=AnSG(n,c):
 mu:={seq(Prof(gu1,n,c),gu1 in gu)}:

for mu1 in mu do
 T[mu1]:={}:
od:

for gu1 in gu do
 pu:=Prof(gu1,n,c):
 T[pu]:=T[pu] union {gu1}:
od:

[mu,op(T)]:

end:

#AnSGP(n,c,P): The set of semi-good order-ideals in OIn(n) with spacing c with profile P. Try
#AnSGP(5,1,[[0,1]]);
AnSGP:=proc(n,c,P) local gu:
option remember:
gu:=ProfDB(n,c):

 if not member(P,gu[1]) then
  {}:
 else
  gu[2][P]:
 fi:
end:


#AnSGPdb(n,c,P):  a database for AnSGPdb(n,c,P), according to its a. Try:
#AnSGdb(7,1,[[0,0]]);
AnSGPdb:=proc(n,c,P) local gu,gu1,T,a:
option remember:
gu:=AnSGP(n,c,P) :

for a from 0 to n do
T[a]:={}:
od:

for gu1 in gu do
T[Dec(gu1,n)[1]]:= T[Dec(gu1,n)[1]] union  {gu1}:
od:

[seq(T[a],a=0..n)]:
end:

#AnSGPa(n,c,P,a): The set of semi-good order-ideals in OIn(n) with spacing c with profile P and whose main diagonal has a elements. Try
#AnSGPa(7,1,[[1,1]],1);
AnSGPa:=proc(n,c,P,a) 
option remember:
if not (0<=a and a<=n) then
 RETURN(FAIL):
else
AnSGPdb(n,c,P)[a+1]
fi:


end:


#KidsProfOld(n,c,P): all the children of the  profile P for AnSGP(n,c,P). Try:
#KidsProfOld(7,1,[[0,0]]);
KidsProfOld:=proc(n,c,P) local gu,gu1,T,a:

gu:=AnSGP(n,c,P) :

gu:={seq(Dec(gu1,n),gu1 in gu)}:

for a from 0 to n do
T[a]:={}:
od:

for gu1 in gu do
T[gu1[1]]:=T[gu1[1]] union {[Prof(gu1[2],a-1,c), Prof(gu1[3],n-a-1,c)]}:
od:

[seq(T[a],a=0..n)]:
end:

#KidsProf(n,c,P,a): all the children of the  profile P for AnSGPa(n,c,P,a). Try:
#KidsProf(7,1,[[0,0]],1);
KidsProf:=proc(n,c,P,a) local lu,lu1,gu,gu1,P1,P2:


lu:=AnSGPa(n,c,P,a) :

gu:={}:

for lu1 in lu do
gu1:=Dec(lu1,n):
if gu1[1]<>a then
 RETURN(FAIL):
fi:

P1:=Prof(gu1[2],a-1,c):
P2:=Prof(gu1[3],n-a-1,c):


gu:=gu union {[a,P1,P2]}:
od:
gu:
end:





#AddOne(L): inputs a list of lists of 0's and 1's reverses them
AddOne:=proc(L) local i,j:
[seq([seq(1-L[i][j],j=1..nops(L[i]))], i=1..nops(L))]:
end:


#JoP(a,P1,P2,c): inputs a non-negative integer a and two profiles P1 and P2, and labelling paramter c (0 or 1)
#Finds the profile of the parent. Try
#JoP(3,[[0,1]],[[0,1],[1,0]]);
JoP:=proc(a,P1,P2,c) local gu,P1a,P2a,i:



if a=0 then
 if P1<>[] then
   RETURN(FAIL):
 else
  RETURN(AddOne(P2)):
 fi:
fi:

if nops(P1)>a-1 then
   RETURN(FAIL):
fi:

if a mod 2=0 then
 P2a:=AddOne(P2):
else
 P2a:=P2:
fi:

	   
if c mod 2=1 then
 P1a:=AddOne(P1):
else
 P1a:=P1:
fi:

P1a:=[[1,a mod 2],op(P1a)]:



if P1a=[] then
 RETURN(P2a):
fi:

if P2a=[] then
 RETURN(P1a):
fi:

gu:=[]:

for i from 1 to min(nops(P1a),nops(P2a)) do

	 
 if P1a[i][2]=P2a[i][1] then
    RETURN(FAIL):
 else
  gu:=[op(gu),[P1a[i][1],P2a[i][2]] ]:
 fi:



od:

if nops(P1a)=nops(P2a) then
 RETURN(gu):

elif nops(P1a)>nops(P2a) then
 RETURN([op(gu),op(nops(P2a)+1..nops(P1a),P1a)]):
elif nops(P1a)<nops(P2a) then

 RETURN([op(gu),op(nops(P1a)+1..nops(P2a),P2a)]):
else
 print(`Something terrible happened`):
fi:

FAIL:

end:

#CheckJoP1(n,c,P,a): checks procedure JoP(a,P1,P2,c) for all the members of KidsProf(n,c,P,a).
#Try:
#CheckJoP1(7,1,[[1,1]]);
CheckJoP1:=proc(n,c,P,a) local pu,pu1,lu:

pu:=KidsProf(n,c,P,a):

lu:={seq(JoP(op(pu1),c),pu1 in pu)}:

evalb(lu={} or lu={P}):



end:

#CheckJoP(n,c): checks procedure JoP(a,P1,P2,c) for all members of AnSG(n,c):
#Try:
#CheckJoP(7,1);
CheckJoP:=proc(n,c) local pu,P,a:
pu:=ProfDB(n,c):

evalb({seq(seq(CheckJoP1(n,c,P,a),a=0..n), P in pu)}={true}):
end:



#SpP111(Z): Inputs a pair [a,b] and outputs the set {[[a,0],[1,b]],[[a,1],[0,b]]
SpP111:=proc(Z) local a,b:
a:=Z[1]: b:=Z[2]:
{ [[a,0],[1,b]], [[a,1],[0,b]]}:
end:

#SpP11(L): Given a list of pairs, finds all their splittings [L1,L2] with nops(L1)=nops(L2)=nops(L)  Try:
#SpP11([[1,1],[0,0]]):
SpP11:=proc(L) local L1,gu,mu,gu1,mu1:
if L=[] then
 RETURN({[[],[]]}):
fi:

L1:=[op(1..nops(L)-1,L)]:

gu:=SpP11(L1):
mu:=SpP111(L[nops(L)]):

{seq(  seq(    [ [op(gu1[1]),mu1[1]], [ op(gu1[2]), mu1[2]  ] ],      gu1 in gu),    mu1 in mu)}:

end:


#SpP1(L,r,s): Given a list of pairs, finds all their splittings [L1,L2] with nops(L1)=r and nops(L2)=s  Try:
#SpP1([[1,1],[0,0],[0,1]],1,3):
SpP1:=proc(L,r,s) local gu,L1,Zanav,gu1:

if nops(L)<>max(r,s) then
 RETURN(FAIL):

elif r=s then
 RETURN(SpP11(L)):

elif r>s then
  L1:=[op(1..s,L)]:
  Zanav:=[op(s+1..r,L)]:
  gu:=SpP11(L1):
 RETURN( { seq ( [ [op(gu1[1]),op(Zanav)], gu1[2]] ,   gu1 in gu)}):

else

  L1:=[op(1..r,L)]:
  Zanav:=[op(r+1..s,L)]:
  gu:=SpP11(L1):
 RETURN( { seq ( [gu1[1], [op(gu1[2]),op(Zanav)]] ,   gu1 in gu)}):
fi:

end:



#SpP(n,c,P,a): Inputs a profile in OI(n) and positive integers n and a, (0<=a<=n) and labelling parameter c
#(0 or 1). Outputs
#the set of pairs {[P1,P2]} where P1 is a profile of OI(a-1) and P2 is a profile of OI(n-a-1).
#Try
#SpP(8,1,[[1,1]],3);
SpP:=proc(n,c,P,a) local mu,mu1,mu1a, mu1b, gu,r,s,lu:

if nops(P)>n then
  RETURN(FAIL):
fi:

if a=0 then
 RETURN({[[],AddOne(P)]}):
fi:

if a=1 then
 if P=[] or P[1][1]=0 then
   RETURN({}):
 fi:
 lu:={[[],[[0,P[1][2]], op(2..nops(P),P)] ]}:
 if P[1][1]=[1,1] then
  lu:= lu union {[[],[]]}:
 fi:

fi:

mu:={}:

for r from 1 to a do
 for s from 0 to n-a-1 do
   if max(r,s)=nops(P) then
    mu:= mu union SpP1(P,r,s):
   fi:
 od:
od:

gu:={}:
for mu1 in mu do
  mu1a:=mu1[1]:
  mu1b:=mu1[2]:

 if mu1a[1]=[1, a mod 2] then
   mu1a:=[op(2..nops(mu1a),mu1a)]:
   if c=1 then
    mu1a:=AddOne(mu1a):
   fi:


 if a mod 2=0 then
   mu1b:=AddOne(mu1b):
 fi:

if JoP(a,mu1a,mu1b,c)<>FAIL then
 gu:=gu union {[mu1a,mu1b]}:
fi:

fi:
od:


gu:
end:

#CheckSpP1(n,c,P,a): checks procedure SpP(n,c,P,a) 
#Try:
#CheckSpP1(7,1,[[1,1]],3);
CheckSpP1:=proc(n,c,P,a) local pu,pu1,lu:

pu:=SpP(n,c,P,a) :


lu:={seq(JoP(a,op(pu1),c),pu1 in pu)} minus {FAIL}:


if lu<>{} and lu<>{P} then
 false:
else
true:
fi:

end:


#CheckSpP(n,P): Checks CheckSpP1(n,c,P,a) for all a and c. Try
#CheckSpP(8,[[1,1]]);
CheckSpP:=proc(n,P)  local a, c:
{seq(seq(CheckSpP1(n,c,P,a) ,a=0..n),c=0..1)}:
end:


#OIP(n,c,P): The set of semi-good order-ideals of P_{n+1,n+2} with labelling parameter c, and profile P.
#Try:
#OIP(3,1,[[1,1]]);
OIP:=proc(n,c,P) local gu,mu,mu1,ku,lu,lu1,i1,i,P1,fu,fu1: 
option remember:

if nops(P)>n then
 RETURN({}):
fi:

if n=-1 then
 RETURN({}):
fi:

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

if n=0 then
 if P=[] then
 RETURN({{}}):
 else
  RETURN({}):
fi:
fi:

#i=0
P1:=AddOne(P):
mu:=OIP(n-1,c,P1):
gu:={seq(Hag(mu1,[0,1]),mu1 in mu)}:

#i=n

if P[1]=[1, n mod 2] then
 P1:=[op(2..nops(P),P)]:
 if c=1 then
  P1:=AddOne(P1):
 fi:
 mu:=OIP(n-1,c,P1):

 ku:={seq([n-i1,i1-1],i1=1..n)}:

  gu:=gu union {seq(mu1 union ku,mu1 in mu)}:

fi:

#start 1<=i<=n-1
for i from 1 to n-1 do

ku:={seq([n-i1,i1-1],i1=1..i)}:

fu:=SpP(n,c,P,i):

for fu1 in fu do
 mu:=OIP(i-1,c,fu1[1]):
 mu:={seq(Hag(mu1,[n-i,0]),mu1 in mu)}:
 mu:={seq(mu1 union ku, mu1 in mu)}:

 lu:=OIP(n-i-1,c,fu1[2]):

lu:={seq(Hag(lu1,[0,i+1]),lu1 in lu)}:

 gu:=gu union {seq(seq(mu1 union lu1,mu1 in mu),lu1 in lu)}:
od:

od:
gu:

end:

#CheckOIP1(n,c,P): checks whether OIP(n,c,P)=AnSGP(n,c,P). Try:
#CheckOIP1(3,1,[[1,1]]);
CheckOIP1:=proc(n,c,P)
if OIP(n,c,P)=AnSGP(n,c,P) then
 RETURN(true):
else
print([P,OIP(n,c,P) minus AnSGP(n,c,P), AnSGP(n,c,P) minus OIP(n,c,P)]):
RETURN(false):
fi:

end:


#CheckOIP(n,c): checks whether OIP(n,c,P)=AnSGP(n,c,P) for all profiles. Try
#CheckOIP(3,1);
CheckOIP:=proc(n,c) local P, pu:
pu:=ProfDB(n,c)[1]:
{seq(CheckOIP1(n,c,P),P in pu)}:
end:


#GoodP1(r): all the good profiles of length r. Try
#GoodP1(4);
GoodP1:=proc(r) local gu,mu,mu1:
option remember:
if r=0 then
  RETURN({[]}):
fi:

if r=1 then
 RETURN({[[1,0]],[[1,1]]}):
fi:

mu:=GoodP1(r-1):

gu:={}:

for mu1 in mu do
 gu:=gu union {[op(mu1),[1-mu1[-1][2],0]],[op(mu1),[1-mu1[-1][2],1]]}:
od:

gu:
end:
#GoodP(r): all the good profiles of length <=r. Try
#GoodP(4);
GoodP:=proc(r) local r1:
{seq(op(GoodP1(r1)),r1=0..r)}:
end:

#OIG(n,c): AnG(n,c) done recursively, using the dynamical programming schme. Try
#OIG(5,1);
OIG:=proc(n,c) local pu,P:

pu:=GoodP(trunc((n+2)/2)):

{seq(op(OIP(n,c,P)), P in pu)}:
end:




#NuOIP(n,c,P): The number of semi-good order-ideals of P_{n+1,n+2} with labelling parameter c, and profile P.
#Try:
#NuOIP(3,1,[[1,1]]);
NuOIP:=proc(n,c,P) local gu,mu,lu,i,P1,fu,fu1: 
option remember:

if nops(P)>n then
 RETURN(0):
fi:

if n=-1 then
 RETURN(0):
fi:

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

if n=0 then
 if P=[] then
 RETURN(1):
 else
  RETURN(0):
fi:
fi:

#i=0
P1:=AddOne(P):
gu:=NuOIP(n-1,c,P1):


#i=n

if P[1]=[1, n mod 2] then
 P1:=[op(2..nops(P),P)]:
 if c=1 then
  P1:=AddOne(P1):
 fi:
 gu:=gu+NuOIP(n-1,c,P1):
fi:

#start 1<=i<=n-1
for i from 1 to n-1 do


fu:=SpP(n,c,P,i):

for fu1 in fu do
 mu:=NuOIP(i-1,c,fu1[1]):
 lu:=NuOIP(n-i-1,c,fu1[2]):

 gu:=gu+mu*lu:
od:

od:
gu:

end:

#NuOIG(n,c): The number of good order ideals of P_{n+1,n+2} with labelling parameter c. Try
#NuOIG(5,1);
NuOIG:=proc(n,c) local pu,P:

pu:=GoodP(trunc((n+2)/2)):

add(NuOIP(n,c,P), P in pu):

end:


#NuOIGr(n,r,c): The number of good order ideals of P_{n+1,n+2} with labelling parameter c and width <=r. Try
#NuOIGr(5,1,1);
NuOIGr:=proc(n,r,c) local pu,P:

pu:=GoodP(r):

add(NuOIP(n,c,P), P in pu):

end:

#NuOIGr1(n,r,c): The number of good order ideals of P_{n+1,n+2} with labelling parameter c and width=r. Try
#NuOIGr1(5,1,1);
NuOIGr1:=proc(n,r,c) local pu,P:

pu:=GoodP1(r):

add(NuOIP(n,c,P), P in pu):

end:


#OAseq(N): The first N terms of the two sequences related to the Odd-Armin sequence and  OEIS sequence A047749, and their
#interleavings. The output is  a list of these four sequences
#Try:
#OAseq(11);
OAseq:=proc(N) local gu0, gu1,mu0,mu1,i:


gu0:=[seq(NuOIG(i,0),i=1..N)]:
gu1:=[seq(NuOIG(i,1),i=1..N)]:
mu0:=[seq( op([gu1[2*i-1],gu0[2*i]]),i=1..trunc(N/2))]:
mu1:=[seq( op([gu0[2*i-1],gu1[2*i]]),i=1..trunc(N/2))]:
print(``):

[gu0,gu1,mu0,mu1]:

end:

#OAseqS(N): the story of two sequences related to the Odd-Armin sequence and  OEIS sequence A047749. Try:
#OAseqS(11);
OAseqS:=proc(N) local gu0, gu1,i:

print(`The sequence enumerating order-ideals in P_{n+1,n+2} that alternate colors, with coloring parameter c=0, up to the `, N, `-th term  is`):
gu0:=[seq(NuOIG(i,0),i=1..N)]:
print(``):
print(gu0):
print(``):


print(`The sequence enumerating order-ideals in P_{n+1,n+2} that alternate colors, with coloring parameter c=1, up to the `, N, `-th term  is`):
gu1:=[seq(NuOIG(i,1),i=1..N)]:
print(``):
print(gu1):
print(``):


print(`The sequence enumerating order-ideals in P_{n+1,n+2} that alternate colors, with coloring parameter n mod 2, up to the `, N, `-th term  is`):
print(`starting with n=1 is`):
print(`This sequence is our object of desire, the elusive sequence enumerating ODD (n+1,n+2)-core partitions `):
print(``):
lprint([seq( op([gu1[2*i-1],gu0[2*i]]),i=1..trunc(N/2))]):
print(``):

print(`The sequence enumerating order-ideals in P_{n+1,n+2} that alternate colors, with coloring parameter n+1 mod 2, up to the `, N, `-th term  is`):
print(`starting with n=1 is`):
print(`This sequence is (presumably) the twin-sequence, the good-looking A047749, defined by`):
print(``):
print(`If n=2m then C(3m,m)/(2m+1); if n=2m+1 then C(3m+1,m+1)/(2m+1) `):
print(``):
lprint([seq( op([gu0[2*i-1],gu1[2*i]]),i=1..trunc(N/2))]):
print(``):

print(`This took`, time(), `seconds `):

end:



#OAseqR(N,r): The first N terms of two sequences related to the restricted Odd-Armin sequence whose support is the r outmost diagonals
#and their
#interleavings. The output is  a list of these four sequences
#Try:
#OAseqR(11,2);
OAseqR:=proc(N,r) local gu0, gu1,mu0,mu1,i:


gu0:=[seq(NuOIGr(i,r,0),i=1..N)]:
gu1:=[seq(NuOIGr(i,r,1),i=1..N)]:
mu0:=[1,seq( op([gu1[2*i-1],gu0[2*i]]),i=1..trunc(N/2))]:
mu1:=[1,seq( op([gu0[2*i-1],gu1[2*i]]),i=1..trunc(N/2))]:
print(``):

[gu0,gu1,mu0,mu1]:

end:

#OAseqR1(N,r): The first N terms of two sequences related to the restricted Odd-Armin sequence whose width is EXACTLY r 
#and their
#interleavings. The output is  a list of these four sequences
#Try:
#OAseqR1(11,2);
OAseqR1:=proc(N,r) local gu0, gu1,mu0,mu1,i:


gu0:=[seq(NuOIGr1(i,r,0),i=1..N)]:
gu1:=[seq(NuOIGr1(i,r,1),i=1..N)]:
mu0:=[seq( op([gu1[2*i-1],gu0[2*i]]),i=1..trunc(N/2))]:
mu1:=[seq( op([gu0[2*i-1],gu1[2*i]]),i=1..trunc(N/2))]:
print(``):

[gu0,gu1,mu0,mu1]:

end:

#StoryO(N,R,t): tells the story of restriced odd (n+1,n+2)-core partitions up to width R, displaying N terms. Try:
#StoryO(30,2,t);
StoryO:=proc(N,R,t) local r,gu:
print(`Rational generating functions for the sequences enumerating (n+1,n+2) Core Partitions into Odd parts`):
print(`whose Corresponding Order-Ideals have width is  <=r for r from 1 to`, R, `and their first`, N, `terms. `):
print(``):
print(`By Shalosh  B. Ekhad `):
print(``):

for r from 1 to R do
print(`For the  width`, r, ` the generating function is `):
  gu:=Gk(r,t):
  print(``):
  print(gu):
  print(``):
  print(`and in Maple format`):
  print(``):
  lprint(gu):
print(``):
print(`The first`, N+1, `terms, staring with n=0, are `):

lprint([seq(coeff(taylor(gu,t=0,N+1),t,i),i=0..N)]):

od:

print(``):
print(`This ends this article, that took`, time(), `seconds to generate. `):

end:



#StoryAD(R,t,N): tells the story of (n+1,n+2)-core partitions where each part can repeat at most k times for k from 1 to R
#it out puts the sequences up to the N-th term
StoryAD:=proc(R,t,N) local k,gu,lu,i:
print(`The  generating functions for sequences enumerating (n+1,n+2) core partitions`):
print(`where each part can repeat at most k times for k from 1 to`, R):
print(``):
print(`By Shalosh  B. Ekhad `):
print(``):

for k from 1 to R do
  gu:=Fk(k,t):
print(``):
print(`The generating function enumerating (n+1,n+2)-core partitions where each part can repeat at most`, k, `times is`):
print(``):
print(gu):
print(``):
print(`and in Maple format`):
print(``):
lprint(gu):
print(``):
print(`For the sake of the OEIS, here are the first`, N+1, `terms, starting with n=0`):
print(``):
gu:=taylor(gu,t=0,N+1):
lprint([seq(coeff(gu,t,i),i=0..N)]):

od:

print(``):
print(`This ends this article, that took`, time(), `seconds to generate. `):

end:

#StoryAll(R,t,N): tells the story of enumerating order ideals of P_{n+1,n+2}  supported in the first k outermost diagonals
#for k from 1 to R
#it out puts the sequences up to the N-th term
StoryAll:=proc(R,t,N) local k,gu,i:
print(`The  generating functions for sequences enumerating Order-Ideals of P_{n+1,n+2}`):
print(`confined to the k outermost diagonals for k from 1 to`, R):
print(``):
print(`By Shalosh  B. Ekhad `):
print(``):

for k from 1 to R do
  gu:=Ak(k,t):
print(``):
print(`The generating function enumerating (n+1,n+2)-core partitions supported in the `, k, `outermost diagonals is`):
print(``):
print(gu):
print(``):
print(`and in Maple format`):
print(``):
lprint(gu):
print(``):
print(`For the sake of the OEIS, here are the first`, N+1, `terms, starting with n=0`):
print(``):
gu:=taylor(gu,t=0,N+1):
lprint([seq(coeff(gu,t,i),i=0..N)]):

od:
print(``):
print(`This ends this article, that took`, time(), `seconds to generate. `):

end:




##################start new OIr


#Wid(S,n): Given an order ideal, S, or  P_{n+1,n+2} (i.e. a member of OI(n)) outputs its width. Try
#Wid({[3,0]},4);
Wid:=proc(S,n) local s:

if S={} then
 RETURN(0):
fi:

n-min(seq(s[1]+s[2],s in S)):
end:














#GkSlow(r,t,N): Slow version, by guessing, using N
#data points, of Gk(r,t).(q.v.) For checking purposes only
#GkSlow(2,t,30);

GkSlow:=proc(r,t,N) local gu,lu:
  gu:=OAseqR(N,r):
 lu:=GuessRec(gu[3]):
 if lu<>FAIL then
   lu:=CtoR(lu,t):
  RETURN(lu):
 else
  RETURN(FAIL):
 fi:
end:




#OI(n): the set of order-ideals of n
OI:=proc(n) local gu,mu,mu1,muA,muB,ku,mu2,i1,i:
option remember:

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

#[n-1,0] unoccupied
mu:=OI(n-1):
gu:={seq(Hag(mu1,[0,1]),mu1 in mu)}:

#[n-1,0], ..., [n-1,i-1] occupied (i<=n)

for i from 2 to n do
 muA:=OI(i-2):
 muA:={seq(Hag(mu1,[n-i+1,0]),mu1 in muA)}:
 ku:={seq([n-i1,i1-1],i1=1..i-1)}:
 muA:={seq(mu1 union ku,mu1 in muA)}:
 muB:=OI(n-i):
  muB:={seq(Hag(mu1,[0,i]),mu1 in muB)}:
 gu:=gu union {seq(seq(mu1 union mu2, mu1 in muA),mu2 in muB)}:
od:

#[n-1,0], ..., [0,n-1] are all occupied 
 muA:=OI(n-1):
 ku:={seq([n-i1,i1-1],i1=1..n)}:
 muA:={seq(mu1 union ku,mu1 in muA)}:
 gu:=gu union muA:
gu:
end:


 
NuOI:=proc(n) local i:
option remember:
if n=0 then
 1:
else
2*NuOI(n-1)+add(NuOI(i-2)*NuOI(n-i),i=2..n):
fi:
end:




#OIr(n,r): the set of order-ideals of n with width <=r
OIr:=proc(n,r) local gu,mu,mu1,muA,muB,ku,mu2,i1,i:
option remember:

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

if r=0 then
RETURN({{}}):
fi:

#[n-1,0] unoccupied
mu:=OIr(n-1,r):
gu:={seq(Hag(mu1,[0,1]),mu1 in mu)}:

#[n-1,0], ..., [n-1,i-1] occupied (i<=n)

for i from 2 to n do
 muA:=OIr(i-2,r-1):
 muA:={seq(Hag(mu1,[n-i+1,0]),mu1 in muA)}:
 ku:={seq([n-i1,i1-1],i1=1..i-1)}:
 muA:={seq(mu1 union ku,mu1 in muA)}:
 muB:=OIr(n-i,r):
  muB:={seq(Hag(mu1,[0,i]),mu1 in muB)}:
 gu:=gu union {seq(seq(mu1 union mu2, mu1 in muA),mu2 in muB)}:
od:

#[n-1,0], ..., [0,n-1] are all occupied 
 muA:=OIr(n-1,r-1):
 ku:={seq([n-i1,i1-1],i1=1..n)}:
 muA:={seq(mu1 union ku,mu1 in muA)}:
 gu:=gu union muA:
gu:
end:

#NuOIr(n,r): the number of order-ideals of n with width <=r
NuOIr:=proc(n,r) local i:
option remember:

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

NuOIr(n-1,r)+NuOIr(n-1,r-1)+
add(NuOIr(i-2,r-1)*NuOIr(n-i,r),i=2..n):
end:

#OIgf(t): the generating function for NuOI(n)
OIgf:=proc(t) local f:
solve(f=(1+t*f)^2,f)[2]:
end:



#Ak(k,t): the generating function for NuOIr(n,k) for a fixed k. Try:
#Ak(k,t);
Ak:=proc(k,t) local gu: option remember:

if k=0 then
 RETURN(1/(1-t)):
else
 gu:=Ak(k-1,t):
 RETURN( normal((1+t*gu)/(1-t-t^2*gu))):
fi:
end:

##################From core.txt
with(combinat):

Ids:=proc(k) local ab,m,S,j: 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:

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

AmSk:=proc(m,S,k) local a,S1,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:

#FkSlow(k,t): inputs a pos. integer k and a variable name t, outputs the generating function
# whose Maclaurin coeff. of t^n is the number of (n+1,n+2)-core partitions, where
#no part is repeated more than k times. Try:
#FkSlow(3,t);
FkSlow:=proc(k,t) local ids,f,eq,var,ab,j,id,x,y,rhs,S1,S2,S,m:
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:

##################End from core.txt





#FAST version of FkSlow
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:


#########Odd parts, OIs confined to first k diagonals


#IDs(r,c): All the IDs of OIs confined to the first r diagonals,
#with diagonal spacing c. A set of pairs of the form [A, B],
#where A is a list of size r, with A[i]=1 iff the ith element along the
#bottom edge (starting from the bottom right corner) is filled,
#and B is a list of size <=r of pairs of parities (0=even, 1=odd).
#(If A[i]=1, then this fixes B[i][1]=1+(k-1)*c mod 2.)
IDs:=proc(k,c) local S1,S,s,f,p,a: option remember:
if k=0 then
	#return {[[1],[[1,0]]], [[1],[[1,1]]], [[0],[]], 
	#  [[0],[[0,0]]], [[0],[[1,0]]], [[0],[[1,1]]], [[0],[[0,1]]]}:
	return {[[],[]]}:
end:
S1:=IDs(k-1,c):
S:={}:

for s in S1 do
	f:=s[1]:
	p:=s[2]:
	a:=1+(k-1)*c mod 2:
	#kth spot filled:
	if f=[] or f[-1]=1  then 
		S:= S union {[[op(f),1],[op(p),[a,0]]], [[op(f),1],[op(p),[a,1]]]}:
	fi:
	#kth spot empty, no new profile row:
	S:= S union {[[op(f),0],p]}:
	#kth spot empty, new profile row:
	if nops(p)=k-1  then 
		for a in {[0,0],[0,1],[1,0],[1,1]} do
		   S:= S union {[[op(f),0],[op(p),a]]}:
		od:
	fi:
	
od:
S:
end:



#IsChild(x,y): Is id x the child of id y?
IsChild:=proc(x,y) local fx,px,fy,py: option remember:
fx:=x[1]: px:=x[2]:
fy:=y[1]: py:=y[2]:
if nops(px) < nops(py)-1 or nops(px)>nops(py) then return false: fi:
#check profile matching on top row
if nops(px)<nops(py) and fy[nops(py)]=1  then
	if py[-1][1]<>py[-1][2] then return false: fi:
elif nops(px)<nops(py)  then 
	return false: 
elif fy[nops(py)]=1 then
	if not px[-1]=[py[-1][1],1-py[-1][2]] then return false: fi:
else
	if not px[-1]=[1-py[-1][1],1-py[-1][2]] then return false: fi:
fi:
if nops(py)=1 then return true: fi:
if nops(py)>=2 and fy[nops(py)]=1 and fx[nops(py)-1]=0 then 
	return false:
fi:
return IsChild([fx, px[1..min(nops(px),nops(py)-1)]], [fy, py[1..nops(py)-1]]):
end:


#CheckP(p): Does the parity profile p correspond to an odd partition?
CheckP:=proc(p) local j:
if nops(p)=0 then return true: fi:
if p[1][1]<>1 then return false: fi:
for j from 2 to nops(p) do
	if p[j-1][2]=p[j][1] then return false: fi:
od:
true:
end:

#Gk1(r,c,t): Generating function for alternating-parity OIs 
#confined to first r diagonals, with spacing c between diagonals
Gk1:=proc(r,c,t) local S,eq,var,f,x,y,rhs,g,id:
S:=IDs(r,c):
var:={seq(f[id],id in S)}:
eq:={f[S[1]]=1+t*f[S[1]]}:
g:=f[S[1]]:
for y in S minus {S[1]} do
	rhs:=0:
	for x in S do
		if IsChild(x,y) then
			rhs:=rhs+t*f[x]:
		fi:
	od:
	eq:=eq union {f[y]=rhs}:
	if CheckP(y[2]) then g:=g+f[y]: fi:
od:
#return solve(eq,var):
normal(subs(solve(eq,var),g)):
end:

#Odd(f,t): odd part of f(t)
Odd:=proc(f,t)
normal((f-subs(t=-t,f))/2):
end:

#Even(f,t): even part of f(t)
Even:=proc(f,t)
normal((f+subs(t=-t,f))/2):
end:

#Gk(r,t): generating function for odd partitions whose OIs are confined to first r diagonals
Gk:=proc(r,t)
normal(Odd(Gk1(r,1,t),t)+Even(Gk1(r,0,t),t)):
end:



########