#######################################################################
##stCorePlus.txt: Save this file as  stCorePlus.txt                   #
## To use it, stay in the                                             #
##same directory, get into Maple (by typing: maple <Enter> )          #
##and then type:  read stCorePlus.txt<Enter>                          #
##Then follow the instructions given there                            #
##                                                                    #
##Written by Anthony Zaleski and Doron Zeilberger, Rutgers University #
#zeilberg at math dot rutgers dot edu                                 #
#######################################################################
 
#Created: Nov. 2017 [Adding to  the previous package stCore.txt 
#This version (adding Odd and Distinct): Nov. 3. 2017
with(linalg):

print(` Created:  Nov. 2017`):
print(` This is stCorePlus.txt `):
print(`It is the Maple package that accompanies the article `):
print(``):
print(`On the Intriguing Problem of Counting (n+1,n+2)-core partitions into Odd Parts`):
print(`by Anthony Zaleski and Doron Zeilberger`):
print(``):
print(`It is an exetension of the Maple package stCore written in 2015 by Doron Zeilberger`):
print(`The paper and the package are`):
print(`available from Zeilberger's website`):
print(``):
print(`Please report bugs to zeilberg at math dot rutgers dot edu`):
print(``):
 print(`The most current version of this  Maple packgage is available from`):
 print(` are  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/tokhniot/stCorePlus.txt  .`):

 print(`The most current version of the article is available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/mamarim/mamarimhtml/stcore.html  .`):

print(`---------------------------------------`):
 print(`For a list of the new procedures regarding odd and/or distint, (added Nov. 2017) type ezraP();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

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 MAIN procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):
print(`---------------------------------------`):

with(combinat):

ezraP:=proc()

if args=NULL then
 print(` The new procedures (added for the new vesrion) , regarding odd and/or distinct partitons, are: IsAD, IsOdd, NuPnn1,`):
 print(`NuPnn1Odd, ParsCoreAD, ParsCoreO, ParsCoreSm,  `):
 print(` ParsCoreSmO, Pnn1, Pnn1Odd `):

else
ezra(args):
fi:

end:

ezraC:=proc()

if args=NULL then
 print(` The Checking procedures are: IsLrk, Pnrk`):

 print(``):

else
ezra(args):
fi:

end:

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: AveAndMoms, Ca, Conj, CheckA, CheckAf, CodePa , Fij, HL, FN,GFcore, GFpa,`):
 print(` GP1, GuessC, GuessC1, GuessCfd , MomData, NuSLrk, PAst,ParsCore, ParsCoreF,  Pars, ParsCore1, PAst1,SLrk, SLrkOdd, stP,  Wt, VecToP, Yafe, Yeled `):

 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are:   Fab,    SipurMoms, SipurMomsFD, stCoreMoms, stCoreMomsFD, VSmoms  `):
 print(` `):


elif nops([args])=1 and op(1,[args])=AveAndMoms then
print(`AveAndMoms(f,x,N): Given a probability generating function`):
print(`f(x) (a polynomial, rational function and possibly beyond) `):
print(` returns a list whose first entry is the average  `):
print(` (alias expectation, alias mean) `):
print(` followed by the variance, the third moment (about the mean) and `):
print(` so on, until the N-th moment (about the mean). `):
print(` If f(1) is not 1, than it first normalizes it by dividing `):
print(`by f(1) (if it is not 0) . `):
print(` For example, try: `):
print(` AveAndMoms(((1+x)/2)^100,x,4); `):


elif nops([args])=1 and op(1,[args])=Ca then
print(`Ca(L,r,q): The Carlitz determinant. Try:`):
print(` Generating function according to area of paths included in L. Try:`):
print(`Ca([4,2,1],2,q);`):

elif nops([args])=1 and op(1,[args])=CheckA then
print(`CheckA(s,t): checks Drew Armstong's conjecture for s,t. Try`):
print(`CheckA(3,5); `):


elif nops([args])=1 and op(1,[args])=CheckAf then
print(`CheckAf(s,t): checks Drew Armstong's conjecture for s,t, the fast way. s and t must be rel. prime and s<t. Try`):
print(`CheckAf(3,5); `):


elif nops([args])=1 and op(1,[args])=CodePa then
print(` CodePa(P,a,b): The code of the path (b,a)-path P. Try: `):
print(` CodePa([[0,0],[1,0],[2,0],[2,1]],1,2); `):
print(` CodePa(VecToP([1,0,1,0,1,1,0,0,1,0,0,0,0]),5,8); `):

elif nops([args])=1 and op(1,[args])=Conj then
print(`Conj(lam): the conjugate of partiton lam. Try:`):
print(`Conj([5,3,2]);`):


elif nops([args])=1 and op(1,[args])=GFcore then
print(`GFcore(S,q): the generating function, in q,  according to size of the set of S-core partitions where S is`):
print(`a set of positive integers mutually relatively prime. Try:`):
print(`GFcore({3,5},q)`):

elif nops([args])=1 and op(1,[args])=Fab then
print(`Fab(a,b,q): the generating function, in q,  of all (a,b)-Dyck paths, according the statistic given by Anderson's bijection. Try:`):
print(`In other words, it is a fast way of doing GFpa(a,b,q);  Try: `):
print(`Fab(3,5,q);`):

elif nops([args])=1 and op(1,[args])=Fij then
print(`Fij(a,b,i,j,q,z): the weight-enumerator of walks from [0,0] to [i,j] staying above the line a*y-b*x=0. Try: `):
print(` Fij(3,5,3,5,q,z); `):

elif nops([args])=1 and op(1,[args])=FN then
print(`FN(L): the frequency notation of the partition L`):

elif nops([args])=1 and op(1,[args])=GFpa then
print(`GFpa(s,t,q): the generating function, in q,  of all (s,t)-Dyck paths, according toe statistic given by Anderson's bijection. Try:`):
print(`GFpa(3,5,q);`):

elif nops([args])=1 and op(1,[args])=GP1 then
print(`GP1(L,x,d): inputs a set of values [a,b], with cardinality >=d+2, a variable x, and a positive integer d, finds`):
print(`a polynomial in x, let's call it P(x), of degree d, such that for every [a,b] in the L, P(a)=b. For example try:`):
print(`GP1({seq([(2*i+1),(2*i+1)^3],i=1..6)},x,3);`):



elif nops([args])=1 and op(1,[args])=GuessC then
print(`GuessC(s,t,i,d): inputs variables s and t, and positive integers i and d, guesses a`):
print(`a polynomial expression of degree<=d in s,t, for the i-th moment (and when i>=2, it is about the mean)`):
print(`of the r.v. "size of (s,t)-core patition" for (s,t) relatively prime `):
print(` Try: `):
print(` GuessC(s,t,1,2); `):

elif nops([args])=1 and op(1,[args])=GuessC1 then
print(`GuessC1(s0,t,i,d): inputs a positive integer s0 and a variable t, and a positive integer d, tries to guess`):
print(`a polynomial expression of degree<=d in t, for the i-th moment (and when i>=2, it is about the mean)`):
print(` of the r.v. "size of (s0,t)-core patition" for t relatively prime to s0. `):
print(` Try: `):
print(` GuessC1(3,t,1,2); `):

elif nops([args])=1 and op(1,[args])=GuessCfd then
print(`GuessCfd(s,c,i,d): inputs a variable s and and a positive integer c , and a positive integer d, tries to guess`):
print(`a polynomial expression of degree<=d in t, for the i-th moment (and when i>=2, it is about the mean) `):
print(` of the r.v. "size of (s0,t)-core patition" for t relatively prime to s0. `):
print(` Try: `):
print(` GuessCfd(s,1,1,3); `):

elif nops([args])=1 and op(1,[args])=IsAD then
print(`IsAD(L,k): Is the partiton L such that the largest number of repeated parts is k? Try:`):
print(`IsAD([4,3,2],1);`):

elif nops([args])=1 and op(1,[args])=IsLrk then
print(`IsLrk(Pa,L,r,k): Does the partition Pa obey the restrictions of (L,r,k)? Try:`):
print(`IsLrk([3,1,1], [[[1,3],{3,4}]],3,1);`):

elif nops([args])=1 and op(1,[args])=IsOdd then
print(`IsOdd(L): Is the partition L odd?`):

elif nops([args])=1 and op(1,[args])=HL then
print(`HL(lam): the list-of lists of hook-lengths of the partition lam. Try:`):
print(`HL([3,2,2]);`):


elif nops([args])=1 and op(1,[args])=MomData then
print(`MomData(N,i): collects data for the i-th moment (for i>=2, about the mean) for all pairs (s,t) relatively prime`):
print(` 2<=s<t<=N. (Divided by (s-1)*(t-1). Try: `):
print(` MomData(4,1); `):

elif nops([args])=1 and op(1,[args])=NuPnn1 then
print(` NuPnn1(n): the number of (n+1,n+2), core partitions done the new way. Try:`):
print(`NuPnn1(5);`):

elif nops([args])=1 and op(1,[args])=NuPnn1Odd then
print(` NuPnn1Odd(n): the number of ODD (n+1,n+2), core partitions done the new way. Try:`):
print(`NuPnn1Odd(5);`):

elif nops([args])=1 and op(1,[args])=NuSLrk then
print(`NuSLrk(L,r,k): nops(SLrk(L,r,k)) done faster. Try: `):
print(` NuSLrk([[[1,6],{7,8}]] ,3,2); `):

elif nops([args])=1 and op(1,[args])=NuSLrkOdd then
print(`NuSLrkOdd(L,r,k): inputs  `):
print(`NuSLrkOdd(L,r,k): nops(SLrkOdd(L,r,k)) done faster. Try: `):
print(` NuSLrkOdd([[[1,6],{7,8}]] ,3,1); `):


elif nops([args])=1 and op(1,[args])=Pars then
print(`Pars(n): the set of partitions of n. Try:`):
print(`Pars(10);`):

elif nops([args])=1 and op(1,[args])=ParsCore then
print(`ParsCore(s,t): inputs  positive integer s,t, that must be relatively prime`):
print(`outputs set of partitions whose hook-lenght are never s and never t . Try:`):
print(`ParsCore(3,5);`):

elif nops([args])=1 and op(1,[args])=ParsCoreAD then
print(`ParsCoreAD(s,t,k): The set of (s,t)-core partitions where each part repeats at most k times. Try:`):
print(`ParsCoreAD(5,6,2);`):

elif nops([args])=1 and op(1,[args])=ParsCoreF then
print(`ParsCoreF(s,t): Faster version of ParsCore(s,t) using the Anderson bijection. Inputs  positive integer s,t, that must be relatively prime`):
print(`outputs set of partitions whose hook-lenght are never s and never t . Try:`):
print(`ParsCoreF(3,5);`):

elif nops([args])=1 and op(1,[args])=ParsCoreO then
print(`ParsCoreO(s,t): The set of (s,t)-core partitions into odd parts. Try:`):
print(`ParsCoreO(5,6);`):

elif nops([args])=1 and op(1,[args])=ParsCoreSm then
print(`ParsCoreSm(s,k): The set of (s+1,s+2)-core partitions such that nops(L)-1+L[1]<=n+(n+1)*(k-1). Try:`):
print(`ParsCoreSm(5,2);`):

elif nops([args])=1 and op(1,[args])=ParsCoreSmO then
print(`ParsCoreSmO(s,k): The set of ODD (s+1,s+2)-core partitions such that nops(L)-1+L[1]<=n+(n+1)*(k-1). Try:`):
print(`ParsCoreSmO(5,2);`):

elif nops([args])=1 and op(1,[args])=ParsCore1 then
print(`ParsCore1(n,S): inputs a positive integer n, and a set of positive integers. Outputs the`):
print(`set of partitions of n whose set of hook-lenghts is disjoint from S. Try:`):
print(`ParsCore1(10,{3});`):

elif nops([args])=1 and op(1,[args])=PAst1 then
print(`PAst1(s,t,i,j): the set of paths from [0,0] to [i,j] staying in the region y/x>s/t. Try:`):
print(`PAst1(3,5,3,5);`):

elif nops([args])=1 and op(1,[args])=PAst then
print(`PAst(s,t): the set of paths from [0,0] to [s,t] staying in the region y/x>s/t. Try:`):
print(`PAst(3,5);`):

elif nops([args])=1 and op(1,[args])=Pnn1 then
print(`Pnn1(n):  Pnn1(n): the set of (n+1,n+2), core partitions done the new way. Try:`):
print(`Pnn1(5);`):

elif nops([args])=1 and op(1,[args])=Pnn1Odd then
print(`Pnn1Odd(n):  Pnn1Odd(n): the set of (n+1,n+2), core partitions into odd parts, done the new way.`):
print(`Should be the same as ParsCoreO(n+1,n_2). Try:`):
print(`Pnn1Odd(5);`):

elif nops([args])=1 and op(1,[args])=Pnrk then
print(`Pnrk(n,r,k): the set of partitions with largest part n, number of parts r, and smallest part k.`):
print(`Try:`):
print(`Pnrk(5,3,2);`):

elif nops([args])=1 and op(1,[args])=SipurMoms then
print(`SipurMoms(s,t,M): Verbose version of stCoreMoms(s,t,M) (q.v.) `):
print(` Try: `):
print(` SipurMoms(s,t,3); `):

elif nops([args])=1 and op(1,[args])=SipurMomsFD then
print(`SipurMomsFD(s,c,M): Verbose version of stCoreMomsFD(s,c,M) (q.v.)`):
print(` Try: `):
print(` SipurMomsFD(s,1,3); `):


elif nops([args])=1 and op(1,[args])=SLrk then
print(`SLrk(L,r,k): inputs  `):
print(` L: a list of pairs [[Interval1,ForbiddenHookLength1], [Interval2,ForbiddenHookLength2], ...] `):
print(` r, and k positive integers `):
print(` outputs: `):
print(` The set of partitions whose largest part is LastInterval[2], number of parts is r, and smallest part is k `):
print(` that has the property that all hook-lenght below Interval1 avoid  ForbiddenHookLength1. `):
print(` all hook-lenght below Interval2 avoid  ForbiddenHookLength2. etc. Try: `):
print(` SLrk([[[1,6],{7,8}]] ,3,2); `):

elif nops([args])=1 and op(1,[args])=SLrkOdd then
print(`SLrkOdd(L,r,k): inputs  `):
print(` L: a list of pairs [[Interval1,ForbiddenHookLength1], [Interval2,ForbiddenHookLength2], ...] `):
print(` r, and k positive integers `):
print(` outputs: `):
print(` The set of partitions into ODD parts, whose largest part is LastInterval[2], number of parts is r, and smallest part is k `):
print(` that has the property that all hook-lengths below Interval1 avoid  ForbiddenHookLength1. `):
print(` all hook-lenght below Interval2 avoid  ForbiddenHookLength2. etc. Try: `):
print(` SLrkOdd([[[1,6],{7,8}]] ,3,2); `):

elif nops([args])=1 and op(1,[args])=SLrkBF then
print(`SLrkBF(L,r,k): Like SLr(L,r,k) (q.v.), but done directly, for checking purposes only. Try: `):
print(` SLrkBF([[[1,6],{7,8}]] ,3,2); `):

elif nops([args])=1 and op(1,[args])=stCoreMoms then
print(`stCoreMoms(a,b,M): inputs variables a,b, and a positive integer M, outputs a list of the conjectured (BUT ABSOLUTELY CERTAIN)`):
print(` regarding the random variable "size of an (a,b)-core partition" `):
print(`expressions for the (i) average (ii) variance  followed by  the third though M's moments about the mean, followed by`):
print(`the limits, expressions in a, for fix a, followed by the limits, where b-a=constant. Warning: M=4 takes a while`):
print(`and M>4 take a long time. Try:`):
print(` stCoreMoms(a,b,3); `):

elif nops([args])=1 and op(1,[args])=stCoreMomsFD then
print(`stCoreMomsFD(a,c,M): inputs a variable a, and a positive integer c, and a positive integer M, `):
print(`Regarding the random variable "size of an (a,a+c)-core partition" (where c is rel. prime to a)`):
print(` outputs a list of the conjectured (BUT ABSOLUTELY CERTAIN) `):
print(` expressions in a, for the (i) average (ii) variance  and the third though M's moments about the mean, followed by `):
print(` the  list consisting of the average and variance (again), and the limits of the standardized moments from  `):
print(`from the third to the M-th. Try `):
print(`stCoreMomsFD(a,1,4);`):

elif nops([args])=1 and op(1,[args])=stP then
print(`stP(s,t): The largest partition living in an s by t rectangle. Try:`):
print(`stP(5,7); `):


elif nops([args])=1 and op(1,[args])=VecToP then
print(`VecToP(v): inputs a list in 0,1 and makes a path that starts at [0,0] with is a vertical step.`):
print(`Try: `):
print(`VecToP([1,1,,0,0,1]); `):


elif nops([args])=1 and op(1,[args])=VSmoms then
print(`VSmoms(N): inputs a positive integer and outputs a`):
print(`list of length N whose first and second entries `):
print(` are the average (1/12) and variance (1/360) of the `):
print(` VS-test distribution, and the the i-th item `):
print(`is the standardized i-th moment for i from 3 to N. Try: `):
print(`VSmoms(N); `):

elif nops([args])=1 and op(1,[args])=Wt then
print(`Wt(a,b,i,j,q,z): The weight of the vertical step [i,j]->[i,j+1] in an (a,b)-generalized Dyck path. Try: `):
print(` Wt(3,5,0,2,q,z); `):

elif nops([args])=1 and op(1,[args])=Yafe then
 print(`Yafe(L): prints nicely a list of lists`):

elif nops([args])=1 and op(1,[args])=Yeled then
print(` Yeled(L,k): Inputs  `):
print(` (i) a list, L, of conditions about forbidden hook-lengths, where L hs length r, say`):
print(` here each component is of the form [[StartingCol,EndingCol], SetOfForbiddenHookLengthForThese Columns] `):
print(` (ii) a positive integer k, <=NumberOfColumns `):
print(`Outputs: `):
print(`The list of restrictions implied for the partitions obtained from removing the bottom row, if that row has k  boxes .`):
print(`Try: `):
print(`Yeled([[[1,6],{7,8}]],3);`):

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

end:


###Start from HISTABRUT
#AveAndMoms(f,x,N): Given a probability generating function
#f(x) (a polynomial, rational function and possibly beyond)
#returns a list whose first entry is the average 
#(alias expectation, alias mean)
#followed by the variance, the third moment (about the mean) and
#so on, until the N-th moment (about the mean).
#If f(1) is not 1, than it first normalizes it by dividing
#by f(1) (if it is not 0) .
#For example, try:
#AveAndMoms(((1+x)/2)^100,x,4);

AveAndMoms:=proc(f,x,N) local mu,F,memu1,gu,i:
mu:=simplify(subs(x=1,f)):

if mu=0 then
print(f, `is neither a prob. generating function nor can it be made so`):
RETURN(FAIL):
fi:

F:=f/mu:


memu1:=simplify(subs(x=1,diff(F,x))):

gu:=[memu1]:

F:=F/x^memu1:

F:=x*diff(F,x):
for i from 2 to N do
 F:=x*diff(F,x):
 gu:=[op(gu),simplify(subs(x=1,F))]:
od:

gu:

end:


#AlphaSeq(f,x,N): Given a probability generating function
#f(x) (a polynomial, rational function and possibly beyond)
#returns a list, of length N, whose 
#(i) First entry is the average
#(ii): Second entry is the variance
#for i=3...N, the i-th entry is the so-called alpha-coefficients
#that is the i-th moment about the mean divided by the
#variance to the power i/2 (For example, i=3 is the skewness
#and i=4 is the Kurtosis)
#If f(1) is not 1, than it first normalizes it by dividing
#by f(1) (if it is not 0) .
#For example, try:
#AlphaSeq(((1+x)/2)^100,x,4);
AlphaSeq:=proc(f,x,N) local gu,i:
 gu:=AveAndMoms(f,x,N):
 if gu=FAIL then
  RETURN(gu):
 fi:

if gu[2]=0 then
print(`The variance is 0`):
  RETURN(FAIL):
fi:

[gu[1],gu[2],seq(gu[i]/gu[2]^(i/2),i=3..N)]:

end:
###end from HISTABRUT
IsPar:= proc(lam)local i:        #this procedure tests that lam valid

 if not type (lam,list) then     #is lam a list?
   RETURN(false):
 fi:

 for i from 1 to nops(lam) do    #are the list elements of lam integers?
   if not (type(lam[i], integer) and lam[i]>0) then
    RETURN(false):
        fi:
 od:

 for i from 1 to nops(lam) -1 do #is lam in the correct (non-increasing) order?
   if lam[i] - lam[i+1] < 0 then
     RETURN(false):
        fi:
 od:

 true:                           #if we pass the tests, then we've got it

 end: #IsPar


#--------------------- Chop -----------------------
Chop:=proc(lam) local i:

for i from 1 to nops(lam) while lam[i]>0 do
od:
i:=i-1:

[op(1..i,lam)]:


end:  #Chop

#--------------------- Conj -----------------------
#Conj(lam): the conjugate of partiton lam. Try:
#Conj([5,3,2]);
Conj:=proc(lam)
local i,lam1,Lam1,n :

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

if IsPar(lam) then

        n:=nops(lam):

        lam1:=Chop([ seq( lam[i]-1, i=1..n) ]):

        Lam1:=Conj(lam1):

        [n, op(Lam1)]:

 else
   ERROR(`lam must be an integer list in non-increasing order`):
fi:

end:  #Conj


#stP(s,t): The largest partition living in an s by t rectangle. Try
#stP(5,7);
stP:=proc(s,t) local i,g,gu:

gu:=[]:

for i from 1 to s-1 do
 g:=trunc(t/s*(s-i)):
 if g=0 then
    RETURN(gu):
  else
 gu:=[op(gu),g]:
 fi:
od:
gu:
end:


qbin:=proc(a,b,q) local i:

if a<b or b<0 then
 0:
else
normal(mul(1-q^(a-i),i=0..b-1)/mul(1-q^i,i=1..b)):
fi:

end:

#Ca2(L,q): The Carlitz determinant for r=2, rediscovered by Amdeberhan and Sergel.
# Generating function according to area of paths included in L. Try:
#Ca2([4,2,1],q);
Ca2:=proc(L,q) local i,j,n:
n:=nops(L):

expand(det([seq([seq(q^(binomial(j-i+1,2))*qbin(L[j]+1,j-i+1,q),j=1..n)],i=1..n)])):

end:


#Ca(L,r,q): The Carlitz determinant  
# Try:
#Ca([4,2,1],r,q);
Ca:=proc(L,r,q) local i,j,n:
n:=nops(L):

expand(q^convert(L,`+`)*det([seq([seq(q^(binomial(i-j,2))*qbin(L[j]+r-1,r-i+j-1,q),i=1..n)],j=1..n)])):

end:


#CheckA(s,t): checks Drew Armstong's conjecture for s,t
CheckA:=proc(s,t) local gu,q:
if gcd(s,t)<>1 then
 print(s,t, `should be relatively prime`):
 RETURN(FAIL):
fi:

gu:=GFcore({s,t},q):

[subs(q=1,diff(gu,q))/subs(q=1,gu),(s+t+1)*(s-1)*(t-1)/24]:

end:

#CheckAf(s,t): checks Drew Armstong's conjecture for s,t
CheckAf:=proc(s,t) local gu,q:

if gcd(s,t)<>1 then
 print(s,t, `should be relatively prime`):
 RETURN(FAIL):
fi:


if s>t then
 print(`s should be less than t`):
 RETURN(FAIL):
fi:


gu:=GFpa(s,t,q):

[subs(q=1,diff(gu,q))/subs(q=1,gu),(s+t+1)*(s-1)*(t-1)/24]:

end:



#HL(lam): the list-of lists of hook-lengths of the partition lam. Try:
#HL([3,2,2]);
HL:=proc(lam) local lam1,i,j:
option remember:
lam1:=Conj(lam):
[seq([seq(lam[i]-j+lam1[j]-i+1,j=1..lam[i])],i=1..nops(lam))]:
end:

#HL1(lam): the set of hook-lengths of the partition lam. Try:
#HL1([3,2,2]);
HL1:=proc(lam) local lam1,i,j:
option remember:
lam1:=Conj(lam):
{seq(seq(lam[i]-j+lam1[j]-i+1,j=1..lam[i]),i=1..nops(lam))}:
end:



Rev:=proc(L) local i: [seq(L[nops(L)-i+1],i=1..nops(L))]:end:

#Pars(n): the set of partitions of n. Try:
#Pars(10);
Pars:=proc(n) local gu,gu1:
option remember:
gu:=partition(n):
{seq(Rev(gu1),gu1 in gu)}:
end:


#ParsCore1(n,S): inputs a positive integer n, and a set of positive integers. Outputs the
#set of partitions of n whose set of hook-lenghts is disjoint from S. Try:
#ParsCore1(10,{3});
ParsCore1:=proc(n1,S) local gu,mu,mu1:
mu:=Pars(n1):

gu:={}:


for mu1 in mu do
 if HL1(mu1) intersect S={} then
   gu:=gu union {mu1}:
 fi:
od:

gu:

end:


#GFcore(S,q): the generating function, in q,  according to size of the set of S-core partitions where S is
#a set of positive integers mutually relatively prime. Try:
#GFcore({3,5},q)
GFcore:=proc(S,q) local ma,i,j,n1:
option remember:
if nops(S)<2 then
 print(`S must have at least two elements`):
 RETURN(FAIL):
fi:

if igcd(S[1],S[2])<>1 then
 print(S,`should be relatively prime`):
 RETURN(FAIL):
fi:

ma:=(S[1]^2-1)*(S[2]^2-1)/24:

if nops(S)>2 then
for i from 1 to nops(S) do
 for j from i+1 to nops(S) do
   if igcd(i,j)<>1 then
     print(i,j, `should be relatively prime `):
     RETURN(FAIL):
   fi:

   ma:=min(ma,(S[i]^2-1)*(S[j]^2-1)/24):
 od:
od:
fi:


add(nops(ParsCore1(n1,S))*q^n1,n1=0..ma):
end:



#GP1(L,x,d): inputs a set of values [a,b], with cardinality >=d+2, a variable x, and a positive integer d, finds
#a polynomial in x, let's call it P(x), of degree d, such that for every [a,b] in the L, P(a)=b. For example try:
#GP1({seq([(2*i+1),(2*i+1)^3],i=1..6)},x,3);
GP1:=proc(L,x,d) local P,a,i,var,eq:
if nops(L)<d+2 then
 print(L, `should have at least`, d+2, `members `):
 RETURN(FAIL):
fi:

P:=add(a[i]*x^i,i=0..d):
var:={seq(a[i],i=0..d)}:

eq:={seq(subs(x=L[i][1],P)-L[i][2],i=1..nops(L))}:

var:=solve(eq,var):

if var=NULL then
  RETURN(FAIL):
else
subs(var,P):
fi:

end:




#GuessC1(s0,t,i,d): inputs a positive integer s0 and a variable t, and a positive integer d, tries to guess
#a polynomial expression of degree<=d in t, for the i-th moment (and when i>=2, it is about the mean)
#of the r.v. "size of (s0,t)-core patition" for t relatively prime to s0.
#Try:
#GuessC1(3,t,1,2);
GuessC1:=proc(s0,t,i,d) local M,t0,q,ku:

M:=[]:



for t0 from 1 while nops(M)<=d+2 do
 if igcd(s0,t0)=1 then
   ku:=[t0, AveAndMoms(Fab(s0,t0,q),q,i)[i]]:
  M:=[op(M), ku ]:
 fi:


od:



GP1(M,t,d):
  
end:


#GuessC(s,t,i,d): inputs variables s and t, and positive integers i and d, guesses a
#a polynomial expression of degree<=d in s,t, for the i-th moment (and when i>=2, it is about the mean)
#of the r.v. "size of (s,t)-core patition" for (s,t) relatively prime 
#Try:
#GuessC(s,t,1,2);
GuessC:=proc(s,t,i,d) local s0,M,ku:

M:=[]:


s0:=2:

while nops(M)<d+2 do

 ku:=GuessC1(s0,t,i,d):


   if ku=FAIL then
      RETURN(FAIL):
   else
  M:=[op(M),[s0,ku]]:
  fi:

   s0:=nextprime(s0):

od:

factor(GP1(M,s,d)):
end:





#MomData(N,i): collects data for the i-th moment (for i>=2, about the mean) for all pairs (s,t) relatively prime
#2<=s<t<=N (divided by (s-1)*(t-1). Try:
#MomData(4,1);
MomData:=proc(N,i) local S,lu,a,b,q,ku:

S:=[]:

for a from 2 to N do
 for b from a+1 to N do
   if gcd(a,b)=1 then
     lu:=GFcore({a,b},q):
      ku:=AveAndMoms(lu,q,i)[i]:
     S:=[op(S),[[a,b],ku/(a-1)/(b-1)]]:
   fi:
 od:
od:
S:
end:

   



#PAst1(s,t,i,j): the set of paths from [0,0] to [i,j] staying in the region y/x>s/t. Try:
#PAst1(3,5,3,5);
PAst1:=proc(s,t,i,j) local lu1,lu2,lu1a,lu2a:
option remember:

if i<0 or j<0 then
  RETURN({}):
fi:

if i=0 and j=0 then
  RETURN({[[0,0]]}):
fi:

if j*s-i*t<0 then
    RETURN({}):
fi:

lu1:=PAst1(s,t,i-1,j):
lu2:=PAst1(s,t,i,j-1):


{seq([op(lu1a),[i,j]],lu1a in lu1), seq([op(lu2a),[i,j]],lu2a in lu2)}:

end:



#IsVer(P1,P2): Is the step from P1 to P2 vertical?
IsVer:=proc(P1,P2)
evalb(P1[1]=P2[1]):
end:

#IsHor(P1,P2): Is the step from P1 to P2 horizontal?
IsHor:=proc(P1,P2)
evalb(P1[2]=P2[2]):
end:


#CodePaNG(P,a,b): The code of the path (b,a)-path P. Try:
#CodePaNG([[0,0],[1,0],[2,0],[2,1]],1,2);
CodePaNG:=proc(L,a,b) local i,j,gu,r,lu,k:
gu:=[]:

for r from  1 to nops(L)-1 do

if IsVer(L[r],L[r+1]) then

i:=L[r][1]:
j:=L[r][2]:

lu:=a*j-b*i-b:
if lu>0 then
 gu:=[op(gu),lu]:
fi:

fi:

if IsHor(L[r],L[r+1]) then
 if (r>1 and IsHor(L[r-1],L[r])) then
  i:=L[r][1]:
  j:=L[r][2]-1:

lu:=a*j-b*i-b:
if lu>0 then
 gu:=[op(gu),lu]:
fi:
fi:
fi:

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


end:




#VecToP(v): inputs a list in 0,1 and makes a path that starts at [0,0] with is a vertical step.
#Try:
#VecToP([1,1,,0,0,1]);
VecToP:=proc(v) local i,gu,pt:
pt:=[0,0]:
gu:=[pt]:

for i from 1 to nops(v) do
 if v[i]=1 then
  pt:=[pt[1],pt[2]+1]:
 gu:=[op(gu),pt]:
elif v[i]=0 then
   pt:=[pt[1]+1,pt[2]]:
 gu:=[op(gu),pt]:
else
 RETURN(FAIL):
 fi:
od:
gu:
end:




#PAst(s,t):
#PAst1(3,5,3,5);
PAst:=proc(s,t):
PAst1(s,t,s,t):
end:



ParsCore:=proc(s,t) local ma,n1:
ma:=(s^2-1)*(t^2-1)/24:
{seq(op(ParsCore1(n1,{s,t})),n1=0..ma)}:
end:


ParsCoreF:=proc(s,t) local gu,gu1:
option remember:
gu:=PAst(s,t):
{seq(CodePa(gu1,s,t),gu1 in gu)}:
end:


#GFpa(s,t,q): the generating function, in q,  of all (s,t)-Dyck paths, according toe statistic given by Anderson's bijection. Try:
#GFpa(3,5,q)
GFpa:=proc(s,t,q) local gu,gu1:
option remember:

if igcd(s,t)<>1 then
 RETURN(FAIL):
fi:


if s>t then 
RETURN(FAIL):
fi:

gu:=ParsCoreF(s,t):

add(q^convert(gu1,`+`),gu1 in gu):

end:





#CodePa(P,a,b): The code of the path (b,a)-path P. Try:
#CodePa([[0,0],[1,0],[2,0],[2,1]],1,2);
CodePa:=proc(L,a,b) local i,j,gu,r,i1,k:
gu:={}:

if igcd(a,b)<>1 then
 RETURN(FAIL):
fi:

if a>b then 
RETURN(FAIL):
fi:


for r from  1 to nops(L)-1 do

if IsVer(L[r],L[r+1]) then

i:=L[r][1]:
j:=L[r][2]:


gu:={op(gu),seq(max(a*j-b*i1-b,0),i1=i..trunc(b/a*j))}:


fi:
od:


gu:=gu minus {0}:


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

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

gu:=[seq(gu[i1]-(k-i1),i1=1..k)]:

for i1 from 1 to nops(gu) while gu[i1]<>0 do od:
[op(1..i1-1,gu)]:

end:







#Wt(a,b,i,j,q,z): The weight of the vertical step [i,j]->[i,j+1] in an (a,b)-generalized Dyck path. Try:
#Wt(3,5,0,2,q,z);
Wt:=proc(a,b,i,j,q,z)  local i1,gu:
gu:=1:

for i1 from i to trunc((a*j-b)/b)+1 do

if a*j-b*i1-b>0 then

  gu:=gu*q^(a*j-b*i1-b)*z:

fi:

od:

gu:

end:


#Fij(a,b,i,j,q,z): the weight-enumerator of walks from [0,0] to [i,j] staying above the line a*y-b*x=0. Try:
Fij:=proc(a,b,i,j,q,z) option remember:

if i<0 or j<0 then
  RETURN(0):
fi:

if i=0 and j=0 then
  RETURN(1):
fi:

if j*a-i*b<0 then
    RETURN(0):
fi:

expand(Fij(a,b,i-1,j,q,z)+Fij(a,b,i,j-1,q,z)*Wt(a,b,i,j-1,q,z)):

end:



#Fab(a,b,q): GFpa(a,b,q) the fast way.
Fab:=proc(a,b,q) local i,gu,z,mu:
option remember:


gu:=Fij(a,b,a,b,q,z):
mu:=0:

for i from 0 to degree(gu,z) do
 mu:=expand(mu+coeff(gu,z,i)*q^(-binomial(i,2))):
od:

sort(mu):

end:



#stCoreMoms(a,b,M): inputs variables a,b, and a positive integer M, 
#Regarding the random variable "size of an (a,b)-core partition"
#outputs a list of the conjectured (BUT ABSOLUTELY CERTAIN)
#expressions for the (i) average (ii) variance  and the third though M's moments about the mean, followed by
#the limits, expressions in a, for fix a, followed by the limits, where b-a=constant. Warning: M=4 takes awhile
#and M>5 take a long time. Try:
#stCoreMoms(a,b,3);
stCoreMoms:=proc(a,b,M) local i,gu,mu,lu:
option remember:
gu:=[]:

for i from 1 to M do
 gu:=[op(gu),GuessC(a,b,i,2*i)]:
od:

mu:=[op(1..2,gu),seq( lcoeff(gu[i],a)/(lcoeff(gu[2],a))^(i/2),i=3..nops(gu))]:

lu:=[op(1..2,gu),seq( simplify(lcoeff(subs(b=a+1,gu[i]),a)/(lcoeff(subs(b=a+1,gu[2]),a))^(i/2)),i=3..nops(gu))]:

[gu,mu,lu]:

end:






#GuessCfd(s,c,i,d): inputs a variable s and and a positive integer c , and a positive integer d, tries to guess
#a polynomial expression of degree<=d in t, for the i-th moment (and when i>=2, it is about the mean)
#of the r.v. "size of (s0,t)-core patition" for t relatively prime to s0.
#Try:
#GuessCfd(s,1,1,3);
GuessCfd:=proc(s,c,i,d) local M,s0,q,ku:

M:=[]:



for s0 from 1 while nops(M)<=d+2 do
 if igcd(s0,c)=1 then
   ku:=[s0, AveAndMoms(Fab(s0,s0+c,q),q,i)[i]]:
   M:=[op(M), ku ]:
 fi:


od:



factor(GP1(M,s,d)):
  
end:


#stCoreMomsFD(a,c,M): inputs a variables a, a positive integer c,, and a positive integer M, 
#Regarding the random variable "size of an (a,a+c)-core partition" (for a relatively prime to a)
#outputs a list of the conjectured (BUT ABSOLUTELY CERTAIN)
#expressions for the (i) average (ii) variance  and the third though M's moments about the mean, followed by
#the limits, expressions in a, for fix a, followed by the limits, where b-a=constant. Warning: M=4 takes awhile
#and M>5 take a long time. Try:
#stCoreMomsFD(a,1,3);
stCoreMomsFD:=proc(a,c,M) local i,gu,mu:
option remember:
gu:=[]:

for i from 1 to M do
 gu:=[op(gu),GuessCfd(a,c,i,3*i)]:
od:

mu:=[op(1..2,gu),seq( simplify(lcoeff(gu[i],a)/(lcoeff(gu[2],a))^(i/2)),i=3..nops(gu))]:

[gu,mu]:

end:






#SipurMoms(s,t,M): Verbose version of stCoreMoms(s,t,M) (q.v.)
#Try:
#SipurMoms(s,t,3);
SipurMoms:=proc(s,t,M) local gu,mu,lu,i,t0:

t0:=time():

if not type(M,integer) and M>=2 then
 print(M, `should have been at least 2`):
 RETURN(FAIL):
fi:

gu:=stCoreMoms(s,t,M):

print(``):
print(`--------------------------------------------------------------------------------`):
print(``):
print(`Explicit Polynomial Expressions for the Average, Variance, and all Moments up to the`, M, `-th `):
print(` for the Random Variable "Size of an (s,t)-core Partition" `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Let s and t be relatively prime positive integers, and consider partitions whose hook-lengths`):
print(`are NEVER s and NEVER t. It has been shown that there are only finitely many of them, namely `):
print(``):
print((s+t-1)!/(s!*t!), `. `):
print(``):
print(`Drew Armstrong conjectured that the average is given by `):
print(``):
print( (s-1)*(t-1)*(s+t+1)/24, `. ` ):
print(``):
print(`This has been first proved (summer 2015) by P. Johsnon and shortly after reproved by Victor Y. Wang. `):
print(``):
print(`In the present article, we will conjecture explicit polynomial expressions for`):
print(`the variance, and all moments (about the mean) up to the`, M, `-th  moment `):
print(`We will also state expressions, in t, for the limiting standardized moments for fixed s`):
print(`and the magic numbers when s and t both tend to infinity`):
print(`Note that while these are, officially, "only"  conjectures, they are ABSOLUTELY CERTAIN, and it may not be worth`):
print(`the effort to prove them "rigorously", unless the proofs give considerable insight. `):

mu:=gu[2]:
lu:=gu[3]:
gu:=gu[1]:

print(`First we rederive, without effort, the fact that the AVERAGE size of an (s,t)-core partition is indeed the one conjectured `):
print(`by Drew Armstrong`):
print(``):
print(`namely `, gu[1]):
print(``):
print(`and in Maple notation: `):
lprint(gu[1]):
print(``):
print(`Going beyond the average (aka mean, aka expectation), the VARIANCE (aka the second moment [about the mean]) is given by the polynomial `):
print(``):
print(gu[2]):
print(`and in Maple notation: `):
lprint(gu[2]):
print(``):

for i from 3 to M do
print(`The `, i, `-th moment [about the mean] is given by the polynomial expression `):
print(``):
print(gu[i]):
print(`and in Maple notation: `):
lprint(gu[i]):
od:


for i from 3 to M do
print(``):
print(`For Fixed t, as s goes to infinity,  the limiting stadardized`, i, `-th moment [about the mean], as an expression in t, is given by  `):
print(``):
print(mu[i]):
print(`and in Maple notation: `):
lprint(mu[i]):
od:
print(``):
print(`Finally, when s and t both become large, the limist of the standardized moments from the third to the`, M, `-th are as follows. `):
print(``):
for i from 3 to M do
print(``):
print(`The limit of the`, i, `-th stadardized moment, as both s and t go to infinity is the number`):
print(``):
print(lu[i], `. `):
print(``):
print(`and in Maple notation. `):
print(``):
lprint(lu[i]):
od:


print(``):
print(`This ends this EMPIRICAL, yet ABSOLUTELY CERTAIN, fascinating article, that took`,  time()-t0, `seconds to produce. `):
print(``):
print(`--------------------------------------------------------------------------------`):
print(``):
end:


#SipurMomsFD(s,c,M): Verbose version of stCoreMomsFD(s,c,M) (q.v.)
#Try:
#SipurMomsFD(s,1,3);
SipurMomsFD:=proc(s,c,M) local gu,mu,i,t0:

t0:=time():

if not type(M,integer) and M>=2 then
 print(M, `should have been at least 2`):
 RETURN(FAIL):
fi:

gu:=stCoreMomsFD(s,c,M):

print(``):
print(`--------------------------------------------------------------------------------`):
print(``):
print(`Explicit Polynomial Expressions for the Average, Variance, and all Moments up to the`, M, `-th `):
print(` for the Random Variable "Size of an ( `, s,s+c, ` ) -core Partition" `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
if c>1 then
print(`Let s be a positive integer, relatively prime to`, c, `and consider partitions whose hook-lengths`):
else
print(`Let s be a positive integer, and consider partitions whose hook-lengths`):
fi:
print(`are NEVER s and NEVER`,  s+c,  `. It has been shown that there are only finitely many of them, namely `):
print(``):
print((s+s+c-1)!/(s!*(s+c)!), `. `):
print(``):
print(`It follows from Drew Armstrong conjecture, first proved (for the special case (s,s+1)),`):
print(`by Richard Stanley and Fabrizio Zanello, that the average is given by `):
print(``):
print( (s-1)*(s+c-1)*(s+s+c+1)/24, `. ` ):
print(``):
print(`In the present article, we will conjecture explicit polynomial expressions for`):
print(`the variance, and all moments (about the mean) up to the`, M, `-th  moment. `):
print(`We will also state the limits of the standardized moments from the third through the`, M, `-th. `):
print(`as s goes to infinity . `):
print(`Note that while these are officially "only"  conjectures, they are ABSOLUTELY CERTAIN, and it may not be worth`):
print(`the effort to prove them "rigorously", unless the proofs give considerable insight. `):

mu:=gu[2]:
gu:=gu[1]:

print(`First we rederive, without effort, the fact that the AVERAGE size of an (s,t)-core partition is indeed the one conjectured `):
print(`by Drew Armstrong`):
print(``):
print(`namely `, gu[1]):
print(``):
print(`and in Maple notation: `):
lprint(gu[1]):
print(``):
print(`Going beyond the average (aka mean, aka expectation), the VARIANCE (aka the second moment [about the mean]) is given by the polynomial `):
print(``):
print(gu[2]):
print(`and in Maple notation: `):
lprint(gu[2]):
print(``):

for i from 3 to M do
print(``):
print(`The `, i, `-th moment [about the mean] is given by the polynomial expression `):
print(``):
print(gu[i]):
print(`and in Maple notation: `):
lprint(gu[i]):
od:

print(``):
print(`Finally, when s  becomes large, the limits of the standardized moments from the third to the`, M, `-th are as follows. `):
print(``):
for i from 3 to M do
print(``):
print(`The limit of the`, i, `-th stadardized moment, as  s goes to infinity is the number`):
print(``):
print(mu[i], `. `):
print(``):
print(`and in Maple notation. `):
print(``):
lprint(mu[i]):
od:


print(``):
print(`This ends this EMPIRICAL, yet ABSOLUTELY CERTAIN, fascinating article, that took`,  time()-t0, `seconds to produce. `):
print(``):
print(`--------------------------------------------------------------------------------`):
print(``):
end:



#VSmomsOld(N): inputs a positive integer and outputs a
#list of length N whose first and second entries
#are the average (1/12) and variance (1/360) of the
#VS-test distribution, and the the i-th item
#is the standardized i-th moment for i from 3 to N. Try:
#VSmomsOld(9);
VSmomsOld:=proc(N) local gu,t,i,m,k:
gu:=exp(add(t^i/2^i/i*Zeta(2*i)/Pi^(2*i),i=1..N)):
gu:=taylor(gu,t=0,N+1):
gu:=[seq(coeff(gu,t,i)*i!,i=1..N)]:
m:=gu[1]:
gu:=[seq(add(binomial(k,i)*(-m)^(k-i)*gu[i],i=1..k)+(-m)^k,k=1..N)]:
gu:=[gu[1],gu[2],seq(gu[i]/gu[2]^(i/2),i=3..N)]:
simplify(gu):
end:

#VSmoms(N): inputs a positive integer and outputs a
#list of length N whose first and second entries
#are the average (1/12) and variance (1/360) of the
#VS-test distribution, and the the i-th item
#is the standardized i-th moment for i from 3 to N. Try:
#VSmoms(9);
VSmoms:=proc(N) local gu,t,i,m,k:
gu:=taylor(t/sin(t),t=0,2*N+2):
gu:=add(coeff(gu,t,i)*t^i,i=0..2*N):
gu:=simplify(subs(t=sqrt(t/2),gu)):
gu:=[seq(coeff(gu,t,i)*i!,i=1..N)]:
m:=gu[1]:
gu:=[seq(add(binomial(k,i)*(-m)^(k-i)*gu[i],i=1..k)+(-m)^k,k=1..N)]:
gu:=[gu[1],gu[2],seq(gu[i]/gu[2]^(i/2),i=3..N)]:
simplify(gu):
end:


###Nes  Stuff, Nov. 3, 2017

#IsOdd(L): Is the partition L odd?
IsOdd:=proc(L) local i:
for i from 1 to nops(L) do
 if L[i] mod 2=0 then
  RETURN(false):
 fi:
od:
true:
end:

#ParsCoreO(s,t): The set of (s,t)-core partitions into odd parts. Try:
#ParsCoreO(5,6);
ParsCoreO:=proc(s,t) local gu,mu,mu1:
option remember:
 mu:=ParsCoreF(s,t):

gu:={}:

for mu1 in mu do
 if IsOdd(mu1) then
  gu:=gu union {mu1}:
 fi:
od:
gu:
end:

#FN(L): the frequency notation of the partition L
FN:=proc(L) local i,L1,a:
option remember:
if L=[] then
 RETURN([]):
fi:

a:=L[1]:

for i from 1 to nops(L) while L[i]=a do od:
L1:=[op(i..nops(L),L)]:

[[a,i-1], op(FN(L1))]:
end:



#IsAD(L,k): Is the partiton L such that the largest number of repeated parts is k? Try
#IsAD([4,3,2],1);
IsAD:=proc(L,k) local i,L1:
L1:=FN(L):

for i from 1 to nops(L1) do
 if L1[i][2]>k then
  RETURN(false):
 fi:
od:

true:
end:

#ParsCoreAD(s,t,k): The set of (s,t)-core partitions where each part repeats at most k times. Try:
#ParsCoreAD(5,6,2);
ParsCoreAD:=proc(s,t,k) local gu,mu,mu1:
option remember:
 mu:=ParsCoreF(s,t):

gu:={}:

for mu1 in mu do
 if IsAD(mu1,k) then
  gu:=gu union {mu1}:
 fi:
od:
gu:
end:



#ParsCoreSm(s,k): The set of (s+1,s+2)-core partitions such that nops(L)-1+L[1]<=n+(n+1)*(k-1). Try:
#ParsCoreSm(5,2);
ParsCoreSm:=proc(s,k) local mu,mu1,gu:
mu:=ParsCoreF(s+1,s+2) minus {[]}:

gu:={[]}:

for mu1 in mu do

  if mu1[1]+nops(mu1)-1<=s+(s+1)*(k-1) then
   gu:=gu union {mu1}:
  fi:

od:

gu:


end:

#ParsCoreSmO(s,k): The set of (s+1,s+2)-core partitions such that nops(L)-1+L[1]<=n+(n+1)*(k-1) and all the parts are odd. Try:
#ParsCoreSmO(5,2);
ParsCoreSmO:=proc(s,k) local mu,mu1,gu:
mu:=ParsCoreF(s+1,s+2) minus {[]}:

gu:={[]}:

for mu1 in mu do

  if IsOdd(mu1) and mu1[1]+nops(mu1)-1<=s+(s+1)*(k-1) then
   gu:=gu union {mu1}:
  fi:

od:

gu:


end:


#Yafe(L): prints nicely a list of lists
Yafe:=proc(L) local i:

for i from 1 to nops(L) do
lprint(op(L[i])):
od:
end:


#MinusOne(A): The set A minus one
MinusOne:=proc(A) local a:
{seq(a-1,a in A)}:
end:






#Yeled(L,k): Inputs 
#(i) a list, L, of conditions about forbidden hook-lengths, where L hs length r, say
#here each component is of the form [[StartingCol,EndingCol], SetOfForbiddenHookLengthForThese Columns]
#(ii) a positive integer k, <=NumberOfColumns
#Outputs
#The list of restrictions implied for the partitions obtained from removing the bottom row, if that row has k  boxes .
#Try:
#Yeled([[[1,6],{7,8}]],3);

Yeled:=proc(L,k) local L1,i,ik,n:

n:=L[nops(L)][1][2]:

if k=n then
 RETURN([seq([L[i][1],MinusOne(L[i][2])],i=1..nops(L))]):
fi:

L1:=[]:

for i from 1 to nops(L) while k>=L[i][1][2] do
  L1:=[op(L1),[L[i][1], MinusOne(L[i][2])]]:
od:



ik:=i:


if ik<=nops(L) then
 if L[ik][1][1]<=k then
  L1:=[op(L1),[[L[ik][1][1],k],MinusOne(L[ik][2])]]:
 fi:

 if k+1<=L[ik][1][2] then
 L1:=[op(L1),[[k+1,L[ik][1][2]],L[ik][2]]]:
fi:
fi:

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

L1:
end:

#SLrk(L,r,k): inputs 
#L: a list of pairs [[Interval1,ForbiddenHookLength1], [Interval2,ForbiddenHookLength2], ...]
#r, and k positive integers
#outputs:
#The set of partitions whose largest part is LastInterval[2], number of parts is r, and smallest part is k
#that has the property that all hook-lenght below Interval1 avoid  ForbiddenHookLength1.
#all hook-lenght below Interval2 avoid  ForbiddenHookLength2. etc. Try:
#SLrk([[[1,6],{7,8}]] ,3,2);
SLrk:=proc(L,r,k) local L1,A1,n,i,j,k1,lu,lu1,i1:
option rememeber:

 n:=L[nops(L)][1][2]:
if r=1 then
 
   if n<>k then
      RETURN({}):
   fi:

 for i from 1 to nops(L) do
   L1:=L[i][1]:
   A1:=L[i][2]:
   if {seq(n-j+1,j=L1[1]..L1[2])} intersect A1<>{} then
    RETURN({}):
   fi:
od:
RETURN({[n]}):

fi:


for i from 1 to nops(L) while k>=L[i][1][1] do
 if {seq(k-i1+1,i1=L[i][1][1]..min(k,L[i][1][2]))} intersect L[i][2]<>{} then
   RETURN({}):
 fi:
od:

L1:=Yeled(L,k):

lu:= {seq(op(SLrk(L1,r-1,k1)),k1=k..n)}:

{seq([op(lu1),k] , lu1 in lu)}:

end:


#Pnrk(n,r,k): the set of partitions with largest part n, number of parts r, and smallest part k.
#Try:
#Pnrk(5,3,2);
Pnrk:=proc(n,r,k) local gu,lu,lu1,k1:
option remember:
 if r=1 then
   if n=k then
     RETURN({[n]}):
   else
    RETURN({}):
   fi:
 fi:

gu:={}:

for k1 from k to n do
  lu:=Pnrk(n,r-1,k1):
  gu:=gu union {seq([op(lu1),k],lu1 in lu)}:
od:
gu:
end:




#IsLrk(Pa,L,r,k): Does the partition Pa obey the restrictions of (L,r,k)? Try:
#IsLrk([3,1,1], [[[1,3],{3,4}]],3,1);
IsLrk:=proc(Pa,L,r,k) local n,i,j,Pa1,a:
n:=Pa[1]:


if L[nops(L)][1][2]<>n then
 RETURN(false):
fi:

if  nops(Pa)<>r then
  RETURN(false):
fi:

 if Pa[nops(Pa)]<>k then
   RETURN(false):
 fi:



Pa1:=Conj(Pa):
  for  a from 1 to nops(L) do
    if {seq(seq(Pa[i]-j+Pa1[j]-i+1,j=L[a][1][1]..min(L[a][1][2],Pa[i])),i=1..nops(Pa))} intersect L[a][2]<>{} then
     RETURN(false):
    fi:
  od:

true:
end:


#SLrkBF(L,r,k): A brute force version of SLrk(L,r,k) (q.v.). For checking puroposes only. Try:
#SLrkBF([[[1,6],{7,8}]] ,3,2);
SLrkBF:=proc(L,r,k) local mu,mu1,n,gu:
n:=L[nops(L)][1][2]:
mu:=Pnrk(n,r,k):

gu:={}:

for mu1 in mu do
  if IsLrk(mu1,L,r,k) then
   gu:=gu union {mu1}:
  fi:
od:

gu:

end:




#Pnn1(n):  Pnn1(n): the set of (n+1,n+2), core partitions done the new way. Try:
#Pnn1(5);
Pnn1:=proc(n) local n1,r1,k1,gu,L:
gu:={[]}:

for n1 from 1 to binomial(n+1,2) do

  L:= [[[1,n1],{n+1,n+2}]]:

   for r1 from 1 to binomial(n+1,2) do
    
    for k1 from 1 to n1 do
      gu:=gu union  SLrk(L,r1,k1):
    od:

  od:

od:

gu:

end:


  


#NuSLrk(L,r,k): inputs 
#L: a list of pairs [[Interval1,ForbiddenHookLength1], [Interval2,ForbiddenHookLength2], ...]
#r, and k positive integers
#outputs:
#The number of elements  in the set of partitions whose largest part is LastInterval[2], number of parts is r, and smallest part is k
#that has the property that all hook-lenght below Interval1 avoid  ForbiddenHookLength1.
#all hook-lenght below Interval2 avoid  ForbiddenHookLength2. etc. Try:
#NuSLrk([[[1,6],{7,8}]] ,3,2);
NuSLrk:=proc(L,r,k) local L1,A1,n,i,j,k1,i1:
option rememeber:

 n:=L[nops(L)][1][2]:
if r=1 then
 
   if n<>k then
      RETURN(0):
   fi:

 for i from 1 to nops(L) do
   L1:=L[i][1]:
   A1:=L[i][2]:
   if {seq(n-j+1,j=L1[1]..L1[2])} intersect A1<>{} then
    RETURN(0):
   fi:
od:
RETURN(1):

fi:


for i from 1 to nops(L) while k>=L[i][1][1] do
 if {seq(k-i1+1,i1=L[i][1][1]..min(k,L[i][1][2]))} intersect L[i][2]<>{} then
   RETURN(0):
 fi:
od:

L1:=Yeled(L,k):

add(NuSLrk(L1,r-1,k1),k1=k..n):

end:

#NuPnn1(n):  the number of (n+1,n+2), core partitions done the new way. Try:
#NuPnn1(5);
NuPnn1:=proc(n) local n1,r1,k1,gu,L:
gu:=1:

for n1 from 1 to binomial(n+1,2) do

  L:= [[[1,n1],{n+1,n+2}]]:

   for r1 from 1 to binomial(n+1,2) do
    
    for k1 from 1 to n1 do
      gu:=gu + NuSLrk(L,r1,k1):
    od:

  od:

od:

gu:

end:


#########################Start Odd
#SLrkOdd(L,r,k): inputs 
#L: a list of pairs [[Interval1,ForbiddenHookLength1], [Interval2,ForbiddenHookLength2], ...]
#r, and k an odd positive integers
#outputs:
#The set of partitions whose largest part is LastInterval[2], and is odd , number of parts is r, and smallest part is k
#that has the property that all hook-lenght below Interval1 avoid  ForbiddenHookLength1.
#all hook-lenght below Interval2 avoid  ForbiddenHookLength2. etc. Try:
#SLrkOdd([[[1,6],{7,8}]] ,3,1);
SLrkOdd:=proc(L,r,k) local L1,A1,n,i,j,k1,lu,lu1,i1:
option rememeber:

if k mod 2 =0 then
  RETURN({}):
fi:

 n:=L[nops(L)][1][2]:
 if n mod 2=0 then
  RETURN({}):
 fi:
 
 if r=1 then
 
   if n<>k then
      RETURN({}):
   fi:

 for i from 1 to nops(L) do
   L1:=L[i][1]:
   A1:=L[i][2]:
   if {seq(n-j+1,j=L1[1]..L1[2])} intersect A1<>{} then
    RETURN({}):
   fi:
od:
RETURN({[n]}):

fi:


for i from 1 to nops(L) while k>=L[i][1][1] do
 if {seq(k-i1+1,i1=L[i][1][1]..min(k,L[i][1][2]))} intersect L[i][2]<>{} then
   RETURN({}):
 fi:
od:

L1:=Yeled(L,k):


lu:= {seq(op(SLrkOdd(L1,r-1,k+2*k1)),k1=0..(n-k)/2)}:

{seq([op(lu1),k] , lu1 in lu)}:

end:

#Pnn1Odd(n):  Pnn1Odd(n): the set of (n+1,n+2), core partitions into odd parts, done the new way. 
#Should be the same as ParsCoreO(n+1,n+2).
#Try:
#Pnn1Odd(5);
Pnn1Odd:=proc(n) local n1,r1,k1,gu,L:
gu:={[]}:

for n1 from 1 to binomial(n+1,2) do

  L:= [[[1,n1],{n+1,n+2}]]:

   for r1 from 1 to binomial(n+1,2) do
    
    for k1 from 1 by 2 to n1 do
      gu:=gu union  SLrkOdd(L,r1,k1):
    od:

  od:

od:

gu:

end:


#NuSLrkOdd(L,r,k): inputs 
#L: a list of pairs [[Interval1,ForbiddenHookLength1], [Interval2,ForbiddenHookLength2], ...]
#r, and k an odd positive integers
#outputs:
#The Number of of ODD partitions whose largest part is LastInterval[2], and is odd , number of parts is r, and smallest part is k
#that has the property that all hook-lenght below Interval1 avoid  ForbiddenHookLength1.
#all hook-lenght below Interval2 avoid  ForbiddenHookLength2. etc. Try:
#NuSLrkOdd([[[1,6],{7,8}]] ,3,1);
NuSLrkOdd:=proc(L,r,k) local L1,A1,n,i,j,k1,i1:
option rememeber:

if k mod 2 =0 then
  RETURN(0):
fi:

 n:=L[nops(L)][1][2]:
 if n mod 2=0 then
  RETURN(0):
 fi:
 
 if r=1 then
 
   if n<>k then
      RETURN(0):
   fi:

 for i from 1 to nops(L) do
   L1:=L[i][1]:
   A1:=L[i][2]:
   if {seq(n-j+1,j=L1[1]..L1[2])} intersect A1<>{} then
    RETURN(0):
   fi:
od:
RETURN(1):

fi:


for i from 1 to nops(L) while k>=L[i][1][1] do
 if {seq(k-i1+1,i1=L[i][1][1]..min(k,L[i][1][2]))} intersect L[i][2]<>{} then
   RETURN(0):
 fi:
od:

L1:=Yeled(L,k):


add(NuSLrkOdd(L1,r-1,k+2*k1),k1=0..(n-k)/2 ):


end:

#NuPnn1Odd(n):  the number of ODD (n+1,n+2), core partitions done the new way. Try:
#NuPnn1Odd(5);
NuPnn1Odd:=proc(n) local n1,r1,k1,gu,L:
gu:=1:

for n1 from 1 to binomial(n+1,2) do

  L:= [[[1,n1],{n+1,n+2}]]:

   for r1 from 1 to binomial(n+1,2) do
    
    for k1 from 1 to n1 do
      gu:=gu + NuSLrkOdd(L,r1,k1):
    od:

  od:

od:

gu:

end: