###############################################################################
# A Maple program written to verify the new partition conjectures             #
# in the article:                                                             #
#                                                                             #
#                                                                             #
#            Staircases to analytic sum-sides for many new integer            #
#                partition identities of Rogers-Ramanujan type                #
#                                                                             #
#                                                                             #
###############################################################################


print(`A Maple program written`):
print(`to verify the new partition conjectures`):
print(`in the article:`):
print(``):
print(`STAIRCASES TO ANALYTIC SUM-SIDES FOR MANY NEW INTEGER PARTITION`): 
print(`IDENTITIES OF THE ROGERS-RAMANUJAN TYPE`):
print(``):
print(`This program requires following Maple packages:`):  
print(`combinat, ListTools, PolynomialsTools.`):
print(``):
print(`Throughout, we adopt the numbering of identities as in the`):
print(`accompanying paper.`):
print(``):
print(`IMPORTANT: Throughout, we let 'q' to be the formal variable.`):
print(`Please make sure that this variable is free.`):
print(``):
print(`For help, type Help();`);


Help:=proc()

if args=NULL then

print(`Main procedure to verify the new identities in the paper is:`):
print(`checkAllNew(cutoff:=500,verbose:=false)`):
print(``):
print(`For a list of procedures type:`):
print(`Help(procedures);`):
print(``):
print(`For more information on a specific procedure, type:`):
print(`Help(<procedure name>);`):
print(`For example:`):
print(`Help(prodmake)`):

elif args=procedures then

print(`Procedures are:`):
print(``):
print(`prodmake(f,t:=30,var:=q)`):
print(`ser(f,terms,var:=q)`):
print(`pochhammer(ap,qpow,n,terms:=-1)`):
print(`inverseEulerTransform(L,upto,terms:=-1)`):
print(`analyticSum3Variables(signfunc,qexp,qpoch,cutoff:=500)`):
print(`analyticSum4Variables(signfunc,qexp,qpoch,cutoff:=500)`):
print(`checkId(idname,cutoff:=500,verbose:=false)`):  
print(`checkAllNew(cutoff:=500,verbose:=false)`):

elif args=prodmake then

print(`prodmake(f,t:=30,var:=q)`):
print(`Given a function 'f' in variable 'var',`): 
print(`find its Eulerian product form up to 't' factors.`):
print(`Return NULL if this can not be done.`):
print(`Try:`):
print(`prodmake(1+q+q^2+q^3+q^4+q^5+q^6,6,q);`):

elif args=ser then

print(`ser(f,terms,var:=q)`):
print(`Given a function 'f' in variable 'var',`):
print(`find its power series around var=0 up to 'terms' terms.`):
print(`Return f as it is if terms=-1.`):
print(`Try:`):
print(`ser(1/(1-q),30);`):

elif args=pochhammer then

print(`pochhammer(ap,qpow,n,terms:=-1)`):
print(`Calculates Pochhammer symbol (ap;qpow)_n expanded as a series up to 'terms'`):
print(`number of terms.`):
print(`Try:`):
print(`pochhammer(q,q^2,10)`):

elif args=inverseEulerTransform then

print(`inverseEulerTransform(L,upto,terms:=-1)`):
print(`First, copy the list 'L' so as to get 'upto' number of terms.`):
print(`suppose the new list is 'Lnew'.`):
print(`then, expand the following product as a series upto 'terms' number of terms.`):
print(`Prod_{i = 1,2,..,upto} (1-q^i)^(-Lnew[i])`):
print(`Try:`):
print(`inverseEulerTransform([0,1],25)`):

elif args=analyticSum3Variables then

print(`analyticSum3Variables(signfunc,qexp,qpoch,cutoff:=500)`):
print(`Makes a q-series as in the paper involving a summation over 3 variables.`):
print(`'signfunc': a function of 3 variables is the sign prefactor in every summand`):
print(`'qexp': a function of 3 variables that is the exponent of q in every summand`):
print(`'qpoch': a sequence of 3 terms of the form q^number appearing in the`): 
print(`         denominators in the summand`):
print(`'cutoff': limit up to which identities need to be verified. `):
print(`The following is Identity 1:`):
print(`analyticSum3Variables((i,j,k)->(-1)^k,`):
print(`                      (i,j,k)->(i+2*j+3*k)*(i+2*j+3*k-1)+i+6*j + 3*k^2+6*k,`):
print(`                      [q,q^4,q^6]);`):

elif args=analyticSum4Variables then

print(`analyticSum4Variables(signfunc,qexp,qpoch,cutoff:=500)`):
print(`Makes a q-series as in the paper involving a summation over 4 variables.`):
print(`'signfunc': a function of 4 variables is the sign prefactor in every summand`):
print(`'qexp': a function of 4 variables that is the exponent of q in every summand`):
print(`'qpoch': a sequence of 4 terms of the form q^number appearing in the`): 
print(`         denominators in the summand`):
print(`'cutoff': limit up to which identities need to be verified. `):
print(`The following is Identity 7:`):
print(`analyticSum4Variables((i,j,k,l)->(-1)^l,`):
print(`                      (i,j,k,l)-> (i+2*j+3*k+4*l)*(i+2*j+3*k+4*l-1)/2+2*l^2`):
print(`                                   + i + 3*j + 6*k + 6*l,`):
print(`                      [q,q^2,q^3,q^4]):`):

elif args=checkId then

print(`checkId(idname, cutoff:=500, verbose:=false)`):
print(`Verify Identity 'idname' upto q^cutoff.`):
print(`'idname' has to be entered as one of the following strings:`):
print(`"1", "2", "3", "4", "4a", "5", "5a", "6", "6a", "7", "7a", "8"`):
print(`For verbose output, set 'verbose' flag to true.`):
print(`Try:`):
print(`checkId("3",100,true);`):

elif args=checkAllNew then

print(`checkAllNew(cutoff:=500, verbose:=false)`):
print(`Check all of the Identities 1, 2, 3, 4, 4a, 5, 5a, 6, 6a, 7, 7a, 8, 8a.`):
print(`For verbose output, set 'verbose' flag to true.`):
print(`Try:`):
print(`checkAllNew(100,true):`):

else 

print(`There is no help for`,args):

fi:

end:



###############################################################################
#                                                                             #
# Required maple packages:                                                    #
#                     combinat, ListTools, PolynomialsTools                   #
#                                                                             #
###############################################################################


with(combinat):
with(ListTools):
with(PolynomialTools):


###############################################################################
#                                                                             #
#                             Helping procedures                              #
#                                                                             #
###############################################################################


###############################################################################
#
# prodmake(f,t:=30,var:=q)
# 
# Given a function `f` in variable `var`, 
# find its Eulerian product form up to `t` factors. 
# Return NULL if this can not be done.
#
# Try:
# prodmake(1+q+q^2+q^3+q^4+q^5+q^6,6,q);

prodmake:=proc(f,t:=30,var:=q) 

local s:=CoefficientList(ser(f,t+1),var),ans:=[],i,temp1,temp2,j,divseq,addup: 

# get a sequence of divisors of ndivseq. 
#
# Try:
# divseq(12);

divseq:=proc(ndivseq) option remember: local ddivseq:
return [seq(`if`(ndivseq mod ddivseq=0, ddivseq, 0),ddivseq=1..ndivseq)];
end:

# Add elements of a list `L`.
#
# Try:
# addup([1,2,3]);

addup:=proc(Laddup) local iaddup: 
return add(iaddup, iaddup in Laddup):
end:

s:=[op(s),seq(0,i=1..t-nops(s)+2)]:

if s[1]<>1 then: return NULL: fi:
s:=s[2..nops(s)]:

for i from 1 to t do:
temp1:= addup(zip(`*`,divseq(i),ans)); 
temp2:=add( s[i-j] * addup(zip(`*`,divseq(j),ans))  ,   j=1..i-1);
ans:=[op(ans), s[i] - (temp1 +temp2)/i]:
od:
return -ans:
end:


###############################################################################
#
# ser(f,terms,var:=q)
# 
# Given a function `f` in variable `var`, 
# find its power series around var=0 up to `terms` terms.
# Return f as it is if terms=-1.
#
# Try:
# ser(1/(1-q),30);

ser:=proc(f,terms,var:=q) local a,an:=0,i:
if (terms=-1) then: return f: fi:
a:=convert(series(f,var=0,terms+1),polynom):
return a:
end:


###############################################################################
#
# pochhammer(ap,qpow,n,terms:=-1)
# 
# Calculates Pochhammer symbol (ap;qpow)_n expanded as a series up to `terms` 
# number of terms.
#
# Try:
# pochhammer(q,q^2,10)

pochhammer:=proc(ap,qpow,n,terms:=-1) local iter,ans,seqx: option remember:
if(n=0) then: return 1: fi:
ans:=[seq(1-ap*(qpow^iter),iter=0..n-1)];
return ser(product(ans[seqx],seqx=1..nops(ans)),terms);
end:


###############################################################################
#
# inverseEulerTransform(L,upto,terms:=-1)
# 
# First, copy the list `L` so as to get `upto` number of terms.
# suppose the new list is `Lnew`.
# then, expand the following product as a series upto `terms` number of terms.
# Prod_{i = 1,2,..,upto} (1-q^i)^(-Lnew[i])
#
# Try:
# inverseEulerTransform([0,1],25)

inverseEulerTransform:=proc(L,upto,terms:=-1) local i, n:=nops(L),ans:=1:
for i from 1 to upto do:
ans:=ser(ans*ser( (1-q^i)^(-L[`if`(i mod n=0,n,i mod n)]) ,terms),terms);
od:
return ans:
end:


###############################################################################
#
# analyticSum3Variables(signfunc,qexp,qpoch,cutoff:=500)
# 
# Makes a q-series as in the paper involving a summation over 3 variables.
#
# `signfunc`: a function of 3 variables is the sign prefactor in every summand
# `qexp`: a function of 3 variables that is the exponent of q in every summand 
# `qpoch`: a sequence of 3 terms of the form q^number appearing in the 
#          denominators in the summand
# `cutoff`: limit up to which identities need to be verified. 
#
# The following is Identity 1:
# analyticSum3Variables((i,j,k)->(-1)^k,
#                       (i,j,k)->(i+2*j+3*k)*(i+2*j+3*k-1)+i+6*j + 3*k^2+6*k,
#                       [q,q^4,q^6]);

analyticSum3Variables:=proc(signfunc,qexp,qpoch,cutoff:=500) 

# `1.5` below is an empirical constant, 
# may need to be fiddled for higher cutoffs / more number of variables

local v1,v2,v3,ans:=0,qe, bound:=ceil(1.5*sqrt(cutoff)):

for v1 from 0 to bound do:
for v2 from 0 to bound do:
for v3 from 0 to bound do:
qe:=qexp(v1,v2,v3):

if(qe>cutoff) then: next: fi:

ans:=ans + ser(signfunc(v1,v2,v3)*q^(qe)
/pochhammer(qpoch[1],qpoch[1],v1)
/pochhammer(qpoch[2],qpoch[2],v2)
/pochhammer(qpoch[3],qpoch[3],v3),cutoff):
od:
od:
od:
return ans:
end:


###############################################################################
#
# analyticSum4Variables(signfunc,qexp,qpoch,cutoff:=500)
# 
# Makes a q-series as in the paper involving a summation over 4 variables.
# For a discription of arguments, see analyticSum3Variables.
#
# The following is Identity 7:
# analyticSum4Variables((i,j,k,l)->(-1)^l,
#                       (i,j,k,l)-> (i+2*j+3*k+4*l)*(i+2*j+3*k+4*l-1)/2+2*l^2
#                                    + i + 3*j + 6*k + 6*l,
#                       [q,q^2,q^3,q^4]):

analyticSum4Variables:=proc(signfunc,qexp,qpoch,cutoff:=500) 

# `1.5` below is an empirical constant, 
# may need to be fiddled for higher cutoffs / more number of variables

local v1,v2,v3,v4,ans:=0,qe, bound:=ceil(1.5*sqrt(cutoff)):

for v1 from 0 to bound do:
for v2 from 0 to bound do:
for v3 from 0 to bound do:
for v4 from 0 to bound do:

qe:=qexp(v1,v2,v3,v4):

if(qe>cutoff) then: next: fi:

ans:=ans + ser(signfunc(v1,v2,v3,v4)*q^(qe)
/pochhammer(qpoch[1],qpoch[1],v1)
/pochhammer(qpoch[2],qpoch[2],v2)
/pochhammer(qpoch[3],qpoch[3],v3)
/pochhammer(qpoch[4],qpoch[4],v4),cutoff):
od:
od:
od:
od:
return ans:
end:


###############################################################################
#
# checkId(idname, cutoff:=500, verbose:=false)
# 
# Verify Identity `idname` upto q^cutoff.
# `idname` has to be entered as one of the following strings:
#     "1", "2", "3", "4", "4a", "5", "5a", "6", "6a", "7", "7a", "8"
# For verbose output, set `verbose` flag to true.
#
# Try:
# checkId("3",100,true);

checkId:=proc(idname,cutoff:=500,verbose:=false) 

local Idsum, Idprod,starttime,endtime,match:

starttime := time();

if idname="1" then

Idsum:=analyticSum3Variables((i,j,k)->(-1)^k,
(i,j,k)->(i+2*j+3*k)*(i+2*j+3*k-1)+i+6*j + 3*k^2+6*k,
[q,q^4,q^6],
cutoff):
Idprod:=inverseEulerTransform([1,0,0,1,0,1,0,1,0,0,1,0],cutoff,cutoff):

elif idname="2" then

Idsum:=analyticSum3Variables((i,j,k)->(-1)^k,
(i,j,k)->(i+2*j+3*k)*(i+2*j+3*k-1)+2*i+2*j + 3*k^2+6*k,
[q,q^4,q^6],
cutoff):
Idprod:=inverseEulerTransform([0,1,1,1,0,-1,0,1,1,1,0,0],cutoff,cutoff):

elif idname="3" then

Idsum:=analyticSum3Variables((i,j,k)->(-1)^k,
(i,j,k)->(i+2*j+3*k)*(i+2*j+3*k-1)+4*i+6*j + 3*k^2+12*k,
[q,q^4,q^6],
cutoff):
Idprod:=inverseEulerTransform([0,0,0,1,1,1,1,1,0,0,0,0],cutoff,cutoff):

elif idname="4" then

Idsum:=analyticSum3Variables((i,j,k)->(-1)^k,
(i,j,k)-> i^2+4*i*j+6*i*k+4*j^2+12*j*k+12*k^2+j,
[q,q^4,q^6],
cutoff):
Idprod:=inverseEulerTransform([1,0,0,1,1,0,0,0,1,0,1,0],cutoff,cutoff):

elif idname="4a" then

Idsum:=analyticSum3Variables((i,j,k)->(-1)^k,
(i,j,k)-> i^2+4*i*j+6*i*k+4*j^2+12*j*k+12*k^2+i-3*j,
[q,q^4,q^6],
cutoff):
Idprod:=inverseEulerTransform([1,0,0,0,1,0,1,1,1,0,0,0],cutoff,cutoff):

elif idname="5" then

Idsum:=analyticSum3Variables((i,j,k)->(-1)^k,
(i,j,k)-> i^2+4*i*j+6*i*k+4*j^2+12*j*k+12*k^2-2*j-3*k,
[q,q^4,q^6],
cutoff):
Idprod:=inverseEulerTransform([1,1,-1,0,1,1,0,0,1,1,0,0],cutoff,cutoff):

elif idname="5a" then

Idsum:=analyticSum3Variables((i,j,k)->(-1)^k,
(i,j,k)-> i^2+4*i*j+6*i*k+4*j^2+12*j*k+12*k^2+i+2*j+3*k,
[q,q^4,q^6],
cutoff):
Idprod:=inverseEulerTransform([0,1,1,0,0,1,1,0,-1,1,1,0],cutoff,cutoff):

elif idname="6" then

Idsum:=analyticSum3Variables((i,j,k)->(-1)^k,
(i,j,k)-> i^2+4*i*j+6*i*k+4*j^2+12*j*k+12*k^2-j,
[q,q^4,q^6],
cutoff):
Idprod:=inverseEulerTransform([1,0,1,0,0,0,1,1,0,0,1,0],cutoff,cutoff):

elif idname="6a" then

Idsum:=analyticSum3Variables((i,j,k)->(-1)^k,
(i,j,k)-> i^2+4*i*j+6*i*k+4*j^2+12*j*k+12*k^2+2*i+3*j+6*k,
[q,q^4,q^6],
cutoff):
Idprod:=inverseEulerTransform([0,0,1,1,1,0,1,0,0,0,1,0],cutoff,cutoff):

elif idname="7" then

Idsum:=analyticSum4Variables((i,j,k,l)->(-1)^l,
(i,j,k,l)-> (i+2*j+3*k+4*l)*(i+2*j+3*k+4*l-1)/2 + 2*l^2
+ i + 3*j + 6*k + 6*l,
[q,q^2,q^3,q^4],
cutoff):
Idprod:=
inverseEulerTransform([1,0,1,1,0,1,0,1,1,0,1,0],cutoff,cutoff):

elif idname="7a" then

Idsum:=analyticSum4Variables((i,j,k,l)->(-1)^l,
(i,j,k,l)-> (i+2*j+3*k+4*l)*(i+2*j+3*k+4*l-1)/2 + 2*l^2
+ 3*i + 5*j + 6*k + 10*l,
[q,q^2,q^3,q^4],
cutoff):
Idprod:=
inverseEulerTransform([0,0,1,1,1,1,1,1,1,0,0,0],cutoff,cutoff):

elif idname="8" then

Idsum:=analyticSum4Variables((i,j,k,l)->(-1)^l,
(i,j,k,l)-> (i+2*j+3*k+4*l)*(i+2*j+3*k+4*l-1)/2 + 2*l^2
+ 2*i + 3*j + 5*k + 6*l,
[q,q^2,q^3,q^4],
cutoff):
Idprod:=
inverseEulerTransform([0,1,1,1,1,0,0,1,1,0,1,0],cutoff,cutoff):

elif idname="8a" then

Idsum:=analyticSum4Variables((i,j,k,l)->(-1)^l,
(i,j,k,l)-> (i+2*j+3*k+4*l)*(i+2*j+3*k+4*l-1)/2 + 2*l^2
+ 2*i + 3*j + 4*k + 6*l,
[q,q^2,q^3,q^4],
cutoff)
+ 
analyticSum4Variables((i,j,k,l)->(-1)^l,
(i,j,k,l)-> (i+2*j+3*k+4*l+1)*(i+2*j+3*k+4*l-1+1)/2 + 2*l^2
+ 2*i + 3*j + 4*k + 6*l+1,
[q,q^2,q^3,q^4],
cutoff)
+ 
analyticSum4Variables((i,j,k,l)->(-1)^l,
(i,j,k,l)-> (i+2*j+3*k+4*l+2)*(i+2*j+3*k+4*l-1+2)/2 + 2*l^2
+ 2*i + 3*j + 4*k + 6*l+2,
[q,q^2,q^3,q^4],
cutoff)
- 
analyticSum4Variables((i,j,k,l)->(-1)^l,
(i,j,k,l)-> (i+2*j+3*k+4*l+3)*(i+2*j+3*k+4*l-1+3)/2 + 2*l^2
+ 2*i + 3*j + 4*k + 6*l+4,
[q,q^2,q^3,q^4],
cutoff):

Idprod:=
inverseEulerTransform([1,0,1,1,0,0,1,1,1,1,0,0],cutoff,cutoff):


else

print(`Identity name must be entered as one of the following strings:`): 
print(`"1", "2", "3", "4", "4a", "5", "5a", "6", "6a", "7", "7a", "8","8a"`):
return():

fi:

match:=Idsum-Idprod:

endtime := time();

if match=0 then:
print("Identity ", idname, " is verified up to partitions of ", cutoff);
print("Verification took ", endtime-starttime, " seconds.");
else:
print("Identity ", idname, " verification FAILED.");
print("Verification took ", endtime-starttime, " seconds.");
fi:

if verbose then:
print("The difference Sum-Product is: ",match);
print("Sum is: ",Idsum);
print("Product is: ",Idprod);
end:

end:


###############################################################################
#
# checkAllNew(cutoff:=500, verbose:=false)
# 
# Check all of the Identities 1, 2, 3, 4, 4a, 5, 5a, 6, 6a, 7, 7a, 8, 8a.
# For verbose output, set `verbose` flag to true.
#
# Try:
# checkAllNew(100,true):

checkAllNew:=proc(cutoff:=500,verbose:=false)
checkId("1" ,cutoff,verbose);
checkId("2" ,cutoff,verbose);
checkId("3" ,cutoff,verbose);
checkId("4" ,cutoff,verbose);
checkId("4a",cutoff,verbose);
checkId("5" ,cutoff,verbose);
checkId("5a",cutoff,verbose);
checkId("6" ,cutoff,verbose);
checkId("6a",cutoff,verbose);
checkId("7" ,cutoff,verbose);
checkId("7a",cutoff,verbose);
checkId("8" ,cutoff,verbose);
checkId("8a",cutoff,verbose);
end: