######################################################################
##HANS: Save this file as  HANS                                      #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read HANS<Enter>                                   #
##Then follow the instructions given there                           #
#                                                                    #
##Written by Andrew V. Sills, Georgia Southern University, and       #
##Doron Zeilberger, Rutgers University ,                             #
#Please report bugs to:                                              #
#zeilberg at math dot rutgers dot edu                                #
###################################################################### 
#Created: Aug-Oct., 2011
#First posted, Sept. 29, 2011
# Version of March 13, 2012 (adding procedure E01s)


with(numtheory):
print(`Created:  Aug.-Oct., 2011`):
print(`This version: March 13, 2012 (adding procedure E01s)`):
print(` This is HANS `):
print(`It accompanies the article `):
print(`"Rademacher's Infinite Partial Fraction Conjecture is False"`):
print(`by Andrew V. Sills and Doron Zeilberger`):
print(`and also available from Sills' and Zeilberger's websites.`):
print(``):
print(`Please report bugs to zeilberg at math dot rutgers dot edu`):
print(``):
 print(`The most current version of this  package and paper`):
 print(` are  available from`):
 print(`http://www.math.rutgers.edu/~zeilberg/  .`):
 print(`For a list of the procedures type ezra();, for help with`):
 print(`a specific procedure, type ezra(procedure_name);   .`):
 print(``):

with(combinat):

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: Ana1, AsyB`):
 print(` C01d, C1hklN, ChkFormulaList, Chkgf, L32, pmngf, RecEqn`):
 print(` SeriesExpansionCoeffs, shk, whk, Zinn `):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: Ana,  ChkFormula, ChklInfinity, ChklN, `):
 print(` ChklNdirect, ChklNfloat,ChklSeq, ChklSeqF,  CloseEncounters `):
 print(`ClosestN, CNi, CNid, `):
 print(` C01, C01f, C01Seq, C01Seqf, C01Table, C01TableF, C01TableMinus `):
 print(` HansTopDownAutoPaper, PaV, PNz, SipurAsy01, SipurAsy011 `):


elif nops([args])=1 and op(1,[args])=Ana then
print(`Ana(L,n): inputs a sequence of real numbers L, and outputs`):
print(`conjenctured asymptotics in terms of a list of length`):
print(`P, where P is the "period" and the ith-term of the`):
print(`output list gives an estimate for the asymptotic behavior of`):
print(`L[P*n+i], for i from 1 to P-1, and the P-th term an`):
print(`for L[P*n]. `):
print(`It also returns the list [max,min] where`):
print(`max is the integer mod P where the local maxima occur`):
print(`and min is the integer mod P where the local minima occur`):
print(`For example, try:`):
print(`read HANSC01: Ana(evalf(C01r[1]),n);`):

elif nops([args])=1 and op(1,[args])=Ana1 then
print(`Ana1(L): analyzes the pattern of differences between peaks and valleys`):
print(`of the sequence L, for example, once you have the data file HANS800, and uploaded it by typing`):
print(`read HANS800: `):
print(` try `):
print(`Ana1(HANS800[1]);`):


elif nops([args])=1 and op(1,[args])=AsyB then
print(`AsyB:=proc(L): Given a list of positive numbers, tries`):
print(`to find an asympototic expression for them as C*mu^n*n^theta`):
print(`it returns the approx. conjenctured asymptotics `):
print(`in the form  C*mu^n*n^theta. `):
print(`AsyB([seq(5*2^i*i^2,i=1..10)],n);`):

elif nops([args])=1 and op(1,[args])=C01Table then
print(`C01Table(K,L): A table of C_{01l}(N) for l from 1 to L`):
print(`and N from 1 to K, it is given as a list of lists`):
print(`such that T[l][N] is C[0,1,l](N). For example, try:`):
print(`C01Table(50,10);`):

elif nops([args])=1 and op(1,[args])=C01TableMinus then
print(`C01TableMinus(K,L): A table of C_{01l}(N) for l from 0 to -L`):
print(`and N from 1 to K, it is given as a list of lists`):
print(`such that T[l+1][N] is C[0,1,-l](N). For example, try:`):
print(`C01TableMinus(50,10);`):

elif nops([args])=1 and op(1,[args])=C01TableF then
print(`C01TableF(K,L,Dig): A table, in floating-point`):
print(`with Dig digits`):
print(`of C_{01l}(N) for l from 1 to L`):
print(`and N from 1 to K, it is given as a list of lists`):
print(`such that T[l][N] is C[0,1,l](N). For example, try:`):
print(`C01TableF(50,10,40);`):

elif nops([args])=1 and op(1,[args])=ChklN then
print(`ChklN(h,k,l,N): the coeff. of 1/(x-exp(2*Pi*h/k))^l in `):
print(`1/(1-x)/(1-x^2)/.../(1-x^N).  computed by the recurrence`):
print(`For example, try:`):
print(`ChklN(0,1,1,5);`):

elif nops([args])=1 and op(1,[args])=ChkFormula then
print(`ChkFormula(h,k,r,n,R): returns a formula, in symbolic n for`):
print(`the Rademacher coefficient C[h,k,n-R](k*n+r) ` ):
print(`For example, try:`):
print(`ChkFormula(1,2,1,n,2)`):

elif nops([args])=1 and op(1,[args])=ChkFormulaList then
print(`ChkFormulaList(h,k,r,n,R): returns a list of formulas for`):
print(`the Rademacher coefficient C[h,k,n-j](k*n+r) for j from 0 to R-1`):
print(`For example try:`):
print(`ChkFormulaList(1,2,1,n,3)`):

elif nops([args])=1 and op(1,[args])=Chkgf then
print(`Chkgf(N,x,h,k): returns the Taylor series, up to`);
print(` O( (x-exp(2*Pi*I*h/k)^N ) of the function`):
print(` x^(floor(N/k))/( (1-x)*(1-x^2)* . . . * (1-x^N) )`):
print(` the first floor(N/k) coefficients of which are `):
print(` Rademacher's coefficients C[h,k,n-j](k*n+r), j=0..n-1 `):
print(` For example, try:`):
print(`Chkgf(6,x,1,4)`):

elif nops([args])=1 and op(1,[args])=ChklInfinity then
print(`ChklInfinity(h,k,l): Rademacher's definition of C[h,k,l](infinity)`):
print(`Eq. (130.6) in p. 302`):
print(`of his book "Topics in Analytic Number Theory"`):
print(`(Springer, 1973, Grundlehren series Band 169)`):
print(`For example, try:`):
print(`ChklInfinity(0,1,1);`):

elif nops([args])=1 and op(1,[args])=ChklNfloat then
print(`ChklNfloat(h,k,l,N,Dig): the coeff. of 1/(x-exp(2*Pi*h/k))^l in `):
print(`1/(1-x)/(1-x^2)/.../(1-x^N).  computed by the recurrence`):
print(`but using floating point with Dig digits.`):
print(`For example, try:`):
print(`ChklNfloat(0,1,1,5,30);`):

elif nops([args])=1 and op(1,[args])=C1hklN then
print(`C1hklN(h,k,l,N): the coeff. of 1/(x-exp(2*Pi*h/k))^(trunc(N/k)-l) in`):
print(`1/(1-x)/(1-x^2)/.../(1-x^N).  computed by the recurrence`):
print(`For example, try:`):
print(`C1hklN(0,1,1,5);`):

elif nops([args])=1 and op(1,[args])=ChklNdirect then
print(`ChklNdirect(h,k,l,N): the coeff. of 1/(x-exp(2*Pi*h/k))^l in`):
print(`1/(1-x)/(1-x^2)/.../(1-x^N).  computed directly`):
print(`For example, try:`):
print(`ChklNdirect(0,1,1,5);`):


elif nops([args])=1 and op(1,[args])=ChklSeq then
print(`ChklSeq(h,k,l,N0):The first N0 terms of the Rademacher sequence`):
print(`C_{hkl}(0,1,1,60). Try: `):
print(`ChklSeq(0,1,1,100);`):

elif nops([args])=1 and op(1,[args])=ChklSeqF then
print(`ChklSeqF(h,k,l,N0,Dig):The first N0 terms of the Rademacher sequence`):
print(`C_{hkl}(N). In floating-point with Dig Digits. Try: `):
print(`ChklSeqF(0,1,1,100,20);`):

elif nops([args])=1 and op(1,[args])=CloseEncounters then
print(`CloseEncounters(K1,L1,N0): The list of the closest`):
print(`approach to Rademacher's C_{hkl}(infinity)`):
print(`together with the deviations for all 0<=h<=k<=K1`):
print(` 1<=l<=L1 for N<=N0. For example try:` ):
print(`CloseEncounters(3,4,50);`):

elif nops([args])=1 and op(1,[args])=ClosestN then
print(`ClosestN(h,k,l,K): The integer N between 1 and K that comes`):
print(`closest to ChklInfinity(h,k,l), together with the difference`):
print(`For example, try:`):
print(`ClosestN(0,1,1,100);`):

elif nops([args])=1 and op(1,[args])=CNi then
print(`CNi(N,i): the coeff. of z^i in z^(N-1)/((z+1)^2-1))/((z+1)^3-1)...`):
print(`by the recurrence. For example, try:`):
print(`CNi(5,4);`):

elif nops([args])=1 and op(1,[args])=CNid then
print(`CNid(N,i): the coeff. of z^i in z^(N-1)/((z+1)^2-1))/((z+1)^3-1)...`):
print(`computed directly. For example, try:`):
print(`CNid(5,4);`):

elif nops([args])=1 and op(1,[args])=C01 then
print(`C01(j,N): The Rademacher coeff. C_{01j}(N)`):
print(`computed by the efficient recurrence. For example, try:`):
print(`C01(1,8);`):

elif nops([args])=1 and op(1,[args])=C01d then
print(`C01d(j,N): The Rademacher coeff. C_{01j}(N)`):
print(`computed directly. For example, try:`):
print(`C01(1,8);`):

elif nops([args])=1 and op(1,[args])=C01f then
print(`C01f(j,N,Dig): The Rademacher coeff. C_{01j}(N) in FLOATING point`):
print(`using Dig Digits`):
print(`could give wrong (and often nonsense) answers if Dig is not`):
print(`large enough. Try:`):
print(`C01f(1,8,20);`):

elif nops([args])=1 and op(1,[args])=C01Seq then
print(`C01Seq(j,N0):The first N0 terms of the Rademacher sequence`):
print(`C_{01j}(N). Try: `):
print(`C01Seq(1,100);`):

elif nops([args])=1 and op(1,[args])=C01Seqf then
print(`C01Seqf(j,N0,d):The first N0 terms of the Rademacher sequence`):
print(`C_{01j}(N), in floating-point to d digits. Try: `):
print(`C01Seqf(1,100,20);`):

elif nops([args])=1 and op(1,[args])=E01s then
print(`E01s(N): the first N terms of the sequence {C_{011}(i)*(2*i)!}`):
print(`Try: `):
print(`E01s(10);`):

elif nops([args])=1 and op(1,[args])=HansTopDownAutoPaper then
print(`HansTopDownAutoPaper(h,k,r,n,C,R): `):
print(`automatically generates a paper which states and proves`):
print(`formulas for Rademacher coefficients C[h,k,n-j](k*n+r) `):
print(`for j from 0 to R-1`):
print(`For example, try:`):
print(`HansTopDownAutoPaper(3,4,0,n,C,3): `):

elif nops([args])=1 and op(1,[args])=L32 then
print(`L32(y1): L_{3/2} in Rademacher's formula on p.302`):
print(`of his book "Topics in Analytic Number Theory"`):
print(`(Springer, 1973, Grundlehren series Band 169)`):
print(`For example, try`):
print(`L32(-4);`):

elif nops([args])=1 and op(1,[args])=PaV then
print(`PaV(L) the locations and values of the local max and min of `):
print(`a sequence. For example, try:`):
print(`PaV([3,2,1,5,7,9,5,4]);`):

elif nops([args])=1 and op(1,[args])=pmngf then
print(`pmngf(m,x): the generating function for the number of partitions`):
print(`of n into at most m parts, factored into cyclotomic polynomials`):
print(`For example, try`):
print(`pmngf(4,x)`):

elif nops([args])=1 and op(1,[args])=PNz then
print(`PNz(N,z): The rational generating function`):
print(` Sum(CNi(N,i)*z^i,i=0..infinity)`):
print(`For example, try:`):
print(`PNz(2,z);`):


elif nops([args])=1 and op(1,[args])=QuasiZinn then 
print(`QuasiZinn(L,p): given a list of numbers L, and a pos. integer p`):
print(`tries to guess, for each n mod p, the appx. asymptotics`):
print(`in the form C*mu^n*n^theta, the output is a list of triples`):
print(`[C,theta,mu] for each residue class`):
print(`For example, try:`):
print(`read HANS500: QuasiZinn(lu[1],32);`):


elif nops([args])=1 and op(1,[args])=RecEqn then
print(`ReqEqn(h,k,r,n,j,C): returns the recurrence equation that results`):
print(`from comparing the coefficients of`, (x-exp(2*Pi*I*h/k))^j ):
print(`on either side of the generating function recurrence equation`):
print(`(1-x^(k*n+r))*(1-x^(k*n+r-1))*...*(1-x^(k*n+r-k+1))* Chkgf(k*n+r,x,h,k) = (x-exp(2*Pi*I*h/k)) * Chkgf(k*n+r-k,x,h,k)`):


elif nops([args])=1 and op(1,[args])=SeriesExpansionCoeffs then
print(`SeriesExpansionCoeffs(h,k,r,n,s): `):
print(`the coefficient of (x- exp(2*Pi*I*h/k))^s in the Taylor expansion`):
print(`of (1-x^(k*n+r))*(1-x^(k*n+r-1))*...* (1-x^(k*n+r-k+1))`):
print(`about x = exp(2*Pi*I*h/k) `):
print(`For example, try`):
print(`SeriesExpansionCoeffs(1,2,0,n,3)`):


elif nops([args])=1 and op(1,[args])=SipurAsy01 then 
print(`SipurAsy01(n): The guessed asymptotics of C01l(32n+r) for l from 1 to 15 and r from 0 to 31.`):
print(`requires file http://www.math.rutgers.edu/~zeilberg/tokhniot/HANSC01f  .`):
print(`Type:`):
print(`SipurAsy01(n):`):

elif nops([args])=1 and op(1,[args])=SipurAsy011 then 
print(`SipurAsy011(n): The guessed asymptotics of C011(32n+r) for r from 0 to 31.`):
print(`requires file http://www.math.rutgers.edu/~zeilberg/tokhniot/oHANS12  .`):
print(`Type:`):
print(`SipurAsy011(n):`):


elif nops([args])=1 and op(1,[args])=shk then 
print(`shk(h,k): The Dedekind sum Eq. (5.2.6) in George Andrews's`):
print(`The Theory of Partitions" (Chap. 5.2, p. 72)`):
print(`For example, try:`):
print(`shk(11,5);`):

elif nops([args])=1 and op(1,[args])=whk then 
print(`whk(h,k): The 24-th root of unity of Rademacher`):

elif nops([args])=1 and op(1,[args])=Zinn then 
print(`Zinn(resh): Zinn-Justin's method to estimate`):
print(`the mu, and theta such that`):
print(`resh[i] is appx. Const*mu^i*i^theta. `):
print(`For example, try:`):
print(`Zinn([seq(5*i*2^i,i=1..30)]);`):

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



ez:=proc(): print(`C01(j,N), C01d(j,N), CNid(N,i), CNi(N,i) `): end:

#C01d(j,N): The Rademacher coeff. C_{01j}(N)
#computed directly
C01d:=proc(j,N) local gu,y,i:
gu:=normal((y-1)^N/mul(1-y^i,i=1..N)):
gu:=subs(y=1+y,gu):
gu:=taylor(gu,y=0,N-j+1):
coeff(gu,y,N-j):
end:


#CNid(N,i): the coeff. of z^i in 1/*(1+(1+z))*(1+(1+z)+(1+z+z^2)*...))
#computed directly
CNid:=proc(N,j) local gu,y,i:
gu:=normal((y-1)^N/mul(1-y^i,i=1..N)):
gu:=subs(y=1+y,gu):
gu:=taylor(gu,y=0,j+1):
coeff(gu,y,j):
end:

#CNi(N,i): the coeff. of z^i in z^(N-1)/((z+1)^2-1))/((z+1)^3-1)...
#by the recurrence. For example, try:
#CNi(5,4);
CNi:=proc(N,i) local r:
option remember:

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

if N<1 then
 RETURN(FAIL):
fi:

if N=1 then
 if i=0 then
   RETURN(-1):
  else
 RETURN(0):
 fi:
fi:

(-CNi(N-1,i)-add(binomial(N,r+1)*CNi(N,i-r),r=1..N-1))/N:

end:

#C01(j,N): The Rademacher coeff. C_{01j}(N)
#computed by the recurrence
C01:=proc(j,N):
CNi(N,N-j):
end:


#CNif(N,i,Dig): the coeff. of z^i in 1/*(1+(1+z))*(1+(1+z)+(1+z+z^2)*...))
#in Floating-point
CNif:=proc(N,i,Dig) local r:
option remember:
Digits:=Dig:
if i<0 then
 RETURN(0):
fi:
if i=0 then
 RETURN(evalf((-1)^N/N!)):
fi:

if N=0 then
 RETURN(0):
fi:

(-CNif(N-1,i,Dig)-add(binomial(N,r+1)*CNif(N,i-r,Dig),r=1..i))/N:

end:


#C01f(j,N,Dig): The Rademacher coeff. C_{01j}(N)
#computed in floating-point
C01f:=proc(j,N,Dig):
CNif(N,N-j,Dig):
end:

#C01Seq(j,N0):The first N0 terms of the Rademacher sequence
#C_{01j}(N). Try: 
#C01Seq(1,100);
C01Seq:=proc(j,N0) local N:
[seq(C01(j,N),N=1..N0)]:
end:

#C01Seqf(j,N0,d):The first N0 terms of the Rademacher sequence
#in floading point with d digits
#C_{01j}(N). Try: 
#C01Seqf(1,100,20);
C01Seqf:=proc(j,N0,d) local N:
[seq(C01f(j,N,d),N=1..N0)]:
end:


#PNzold(N,z): The rational generating function Sum(CNi(N,i)*z^i,i=0..infinity)
#For example, try:
#PNzold(2,z);
PNzold:=proc(N,z) local r:
option remember:
if N=0 then
 RETURN(1):
else
-PNzold(N-1,z)/(add(binomial(N,r+1)*z^r,r=0..N)):
fi:
end:

#PNz(N,z): The rational generating function Sum(CNi(N,i)*z^i,i=0..infinity)
#For example, try:
#PNz(2,z);
PNz:=proc(N,z) option remember:
if N=0 then
 RETURN(1):
else
-PNz(N-1,z)/normal (((1+z)^N-1)/z) :
fi:
end:





PaV1:=proc(L) local i:

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

if L[2]>=L[1] then
for i from 3 to nops(L) while L[i]>=L[i-1] do od:
 RETURN([i-1,L[i-1]]):
else
for i from 3 to nops(L) while L[i]<L[i-1] do od:
 RETURN([i-1,L[i-1]]):
fi:
end:

#PaVa(L) the locations and values of the local max and min of 
#a sequence. For example, try:
#PaVa([3,2,1,5,7,9,5,4]);
PaVa:=proc(L) local i,mu,L1,gu1:

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

mu:=PaV1(L):

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

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


gu1:=PaVa(L1):

[mu,seq([gu1[i][1]+mu[1]-1,gu1[i][2]],i=1..nops(gu1))]:

end:

#PaV(L) the locations and values of the local max and min of 
#a sequence. For example, try:
#PaV([3,2,1,5,7,9,5,4]);
PaV:=proc(L) local gu:
gu:=PaVa(L):
[op(1..nops(gu)-1,gu)]:
end:


sn:=proc(resh,n1):
-1/log(op(n1+1,resh)*op(n1-1,resh)/op(n1,resh)^2):
end:
 
#Zinn(resh): Zinn-Justin's method to estimate
#the C,mu, and theta such that
#resh[i] is appx. Const*mu^i*i^theta
#For example, try:
#Zinn([seq(5*i*2^i,i=1..30)]);
Zinn:=proc(resh)
local s1,s2,theta,mu,n1,i:
if nops({seq(sign(resh[i]),i=1..nops(resh))})<>1 then
 RETURN(FAIL):
fi:

n1:=nops(resh)-1:
s1:=sn(resh,n1):
s2:=sn(resh,n1-1):
theta:=evalf(2*(s1+s2)/(s1-s2)^2):
mu:=evalf(sqrt(op(n1+1,resh)/op(n1-1,resh))*exp(-(s1+s2)/((s1-s2)*s1))):
[theta,mu]:
end:
 
 
 
#QuasiZinn(L,p): given a list of numbers L, and a pos. integer p
#tries to guess, for each n mod p, the appx. asymptotics
#in the form C*mu^n*n^theta, the output is a list of triples
#[C,theta,mu] for each residue class.
#For example, try:
#read HANS500: QuasiZinn(lu[1],32);

QuasiZinn:=proc(L,p) local lu,mu,i,N,T:
N:=nops(L):
if N-4*p-1<1 then
 print(`List too small`):
 RETURN(FAIL):
fi:

for i from 0 to p-1 do
mu:=evalf([L[N-3*p-i],L[N-2*p-i],L[N-p-i], L[N-i]]):
lu:=Zinn(mu):
T[(N-i) mod p]:=[L[N-i]/lu[2]^((N-i)/p)/(N-i)^lu[1], lu[1],lu[2]^(1/p)]:
od:

[seq(T[i],i=1..p-1),T[0]]:

end:




#ChklNdirect(h,k,l,N): the coeff. of 1/(x-exp(2*Pi*h/k))^l in
#1/(1-x)/(1-x^2)/.../(1-x^N). For example, try:
#ChklNdirect(0,1,1,5);
ChklNdirect:=proc(h,k,l,N) local gu,i,w,h1,x,gu1,r,co:
co:=0:
w:=exp(2*Pi*I/k):
gu:=1:
for i from 1 to N do
if i mod k<>0 then
gu:=gu/(x^i-1):
else
gu1:=1:
co:=co+1:
for r from 0 to i-1 do
 
 if r/i<>h/k then
    gu1:=gu1*(x-exp(2*Pi*I*r/i)):
 fi:
od:
gu:=gu/gu1:
fi:
od:
gu:=subs(x=x+w^h,gu):

(-1)^N*coeff(taylor(gu,x=0,co-l+3),x,co-l):

end:


#C1hklN(h,k,l,N): the coeff. of 1/(x-exp(2*Pi*h/k))^(trunc(N/k)-l) in
#1/(1-x)/(1-x^2)/.../(1-x^N). Using the recurrence
#For example, try:
#C1hklN(0,1,1,5);
C1hklN:=proc(h,k,l,N) local om,a,lu,r:
option remember:

om:=exp(2*Pi*I*h/k):

if N=0 then
 if l=0 then
   RETURN(1):
 else
  RETURN(0):
 fi:
fi:

if l<0 then
 RETURN(0):
fi:

if N mod k=0 then
 lu:=om/N*C1hklN(h,k,l,N-1)-
add(binomial(N,a+1)/N/om^a*C1hklN(h,k,l-a,N),a=1..l):
 lu:=evalc(lu):
RETURN(lu):
else
r:=N mod k:
lu:=C1hklN(h,k,l,N-1) - add(binomial(N,a)*om^(r-a)*C1hklN(h,k,l-a,N),
a=1..l):
lu:=lu/(om^r-1):
lu:=evalc(lu):
RETURN(lu):
fi:
end:





#ChklN(h,k,l,N): the coeff. of 1/(x-exp(2*Pi*h/k))^l in
#1/(1-x)/(1-x^2)/.../(1-x^N). Using the recurrence
#For example, try:
#ChklN(1,2,1,5);
ChklN:=proc(h,k,l,N):
(-1)^N*C1hklN(h,k,trunc(N/k)-l,N):
end:




#C1hklNfloat(h,k,l,N,Dig): the coeff. of 1/(x-exp(2*Pi*h/k))^(trunc(N/k)-l) in
#1/(1-x)/(1-x^2)/.../(1-x^N). Using the recurrence with floating point
#with Dig Digits
#For example, try:
#C1hklNfloat(0,1,1,5,20);
C1hklNfloat:=proc(h,k,l,N,Dig) local om,a,lu,r:
option remember:
Digits:=Dig:
om:=evalf(exp(2*Pi*I*h/k)):

if N=0 then
 if l=0 then
   RETURN(1.):
 else
  RETURN(0):
 fi:
fi:

if l<0 then
 RETURN(0):
fi:

if N mod k=0 then
 lu:=om/N*C1hklNfloat(h,k,l,N-1,Dig)-
add(binomial(N,a+1)/N/om^a*C1hklNfloat(h,k,l-a,N,Dig),a=1..l):
 lu:=evalf(lu):
RETURN(lu):
else
r:=N mod k:
lu:=C1hklNfloat(h,k,l,N-1,Dig) - 
add(binomial(N,a)*om^(r-a)*C1hklNfloat(h,k,l-a,N,Dig),
a=1..l):
lu:=lu/(om^r-1):
lu:=evalf(lu):
RETURN(lu):
fi:
end:

#ChklNfloat(h,k,l,N,Dig): the coeff. of 1/(x-exp(2*Pi*h/k))^l in
#1/(1-x)/(1-x^2)/.../(1-x^N). Using the recurrence
#For example, try:
#ChklNfloat(1,2,1,5,30);
ChklNfloat:=proc(h,k,l,N,Dig):
(-1)^N*C1hklNfloat(h,k,trunc(N/k)-l,N,Dig):
end:

#ChklSeq(h,k,l,N0):The first N0 terms of the Rademacher sequence
#C_{hkl}(0,1,1,60). Try: 
#ChklSeq(0,1,1,100);
ChklSeq:=proc(h,k,l,N0) local N:
[seq(ChklN(h,k,l,N),N=1..N0)]:
end:

#ChklSeqF(h,k,l,N0,Dig):The first N0 terms of the Rademacher sequence
#C_{hkl}(0,1,1,60) in floating-point. Try: 
#ChklSeqF(0,1,1,100,20);
ChklSeqF:=proc(h,k,l,N0,Dig) local N:
[seq(ChklNfloat(h,k,l,N,Dig),N=1..N0)]:
end:




#shk(h,k): The Dedekind sum Eq. (5.2.6) in George Andrews's
#The Theory of Partitions" (Chap. 5.2, p. 72)
#For example, try:
#shk(11,5);
shk:=proc(h,k) local m:
add((m/k-trunc(m/k)-1/2)*(h*m/k-trunc(h*m/k)-1/2),m=1..k-1):
end:

#exp(Pi*I*shk(h,k)): The Dedekind sum Eq. (5.2.6) in George Andrews's
whkF:=proc(h,k):
exp(Pi*I*shk(h,k)):
end:

#whk(h,k): The 24-th root of unity of Rademacher
whk:=proc(h,k) local h1:
for h1 from 1 to k-1 while h*h1+1 mod k<>0 do od:
 
if h mod 2=1 then
legendre(-k,h)*
exp(
-Pi*I*
(
1/4*(2-h*k-h)
+1/12*(k-1/k)*(2*h-h1+h^2*h1)
)
):
elif k mod 2=1 then
legendre(-h,k)*
exp(
-Pi*I*
(
1/4*(k-1)
+1/12*(k-1/k)*(2*h-h1+h^2*h1)
)
):
fi:

end:


#L32(y1): L_{3/2} in Rademacher's formula on p.302
#of his book "Topics in Analytic Number Theory"
#(Springer, 1973, Grundlehren series Band 169)
#For example, try
#L32(-4);
L32:=proc(y1) local y:
if evalf(y1)>=0 then
 print(`argument of L32 must he negative`):
RETURN(FAIL):
fi:

y:=sqrt(-y1):
-(2*cos(2*y)-sin(2*y)/y)/2/sqrt(Pi)/y^2:
end:


#ChklInfinity(h,k,l): Rademacher's definition of C[h,k,l](infinity)
#Eq. (130.6) in p. 302
#of his book "Topics in Analytic Number Theory"
#(Springer, 1973, Grundlehren series Band 169)
#For example, try
#For example, try
#ChklInfinity(0,1,1);
ChklInfinity:=proc(h,k,l) local gu,a,mu,alpha:
gu:=-2*Pi*(Pi/12)^(3/2)*whkF(h,k)*exp(2*Pi*I*h*l/k)/k^(5/2):

alpha:=1/24:
mu:=add(binomial(l-1,a)*(-1)^(l-1-a)*L32(-Pi^2/6/k^2*(alpha+1+a)),
a=0..l-1):


simplify(evalc(gu*mu)):

end:


#ClosestN(h,k,l,K): The integer N between 1 and K that comes
#closest to ChklInfinity(h,k,l), together with the difference
#For example, try:
#ClosestN(0,1,1,100);
ClosestN:=proc(h,k,l,K) local lu,i,halev,aluf,si,si1:
lu:=ChklInfinity(h,k,l):
aluf:=1:
si:=evalf(abs(ChklNfloat(h,k,l,1)-lu)):
if lu<>0 then
si1:=evalf(abs(ChklNfloat(h,k,l,1)/lu)):
else
 si1:=FAIL:
fi:

for i from 2 to K do
halev:=evalf(abs(ChklNfloat(h,k,l,i)-lu)):
if halev<si then
  si:=halev:
   aluf:=i:
if lu<>0 then
si1:=evalf(abs(ChklNfloat(h,k,l,i)/lu)):
else
 si1:=FAIL:
fi:

fi:
od:

aluf,si,si1:
end:


#CloseEncounters(K1,L1,N0): The list of the closest
#approach to Rademacher's C_{hkl}(infinity)
#together with the deviations for all 0<=h<=k<=K1
#1<=l<=L1 for N<=N0. For example try:
#CloseEncounters(3,4,50);
CloseEncounters:=proc(K1,L1,N0) local N,i,x,h,k,l,lu,t0:

t0:=time():
print(`Close Encounters of the Rademacher kind`):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`The great analytic number theorist, Hans Rademacher`):
print(`in his posthumously published book`):
print(`"Topics in Analytic Number Theory", Springer, 1973 `):
print(`(pp. 299-302) conjectured that C_{hkl}(N),`):
print(`the coefficient of (x-exp(2*Pi*I*h/k)^(-l) `):
print(`in the partial fraction decomposition of`):
print(1/Prod(1-x^i,i=1..N)):
print(`converges to a certain number defined in his book and`):
print(` implemented in procedure ChklInfinity(h,k,l) in`):
print(`The Maple package HANS, available, free of charge from`):
print(`http://www.math.rutgers.edu/~zeilberg/tokhniot/HANS  .`):
print(`It turns out (empirically), that these sequences diverge very`):
print(`badly, and as N goes to infinity, sooner or later they oscillate`):
print(`widely with peaks that go to infinity (exponentially fast!) and`):
print(`valleys that go to negative infinity (also exponentially fast!)`):
print(``):
print(`But it seems that one can tweak Hans's conjecture`):
print(`if instead of going all the way to infinity, you only go so far.`):
print(`In this article we will find, for h and k, 1<=h<k<=`, K1 ):
print(`with gcd(h,k)=1, and l<=`, L1 , `, `):
print(`that integer N, between 1 and`, N0):
print(`such that C_{hkl}(N) gets closest to Rademacher's conjectured`):
print(`value for the "limit", what he called, with wishful thinking`):
print(`C_{hkl}(infinity) `):

for l from 1 to L1 do
lu:=ClosestN(0,1,l,N0):
print(`For `, C_[0,1,l](N), `N=`, lu[1] , `gets closest to`):
print(C_[0,1,l](infinity)=evalf(ChklInfinity(0,1,l))):
print(`the difference being`, lu[2]):
print(`and the ratio being`, lu[3]):
od:

for k from 2 to K1 do
for h from 1 to k-1 do
 if gcd(h,k)=1 then
  for l from 1 to L1 do
  lu:=ClosestN(h,k,l,N0):
  print(`The closest`, C_[h,k,l](N), ` gets to`):
  print(C_[h,k,l](infinity)=evalf(ChklInfinity(h,k,l))):
  print(`is when `, N=lu[1], ` , `):
  print(`the abs. value of the difference being`, lu[2]):
  print(`and the absolute value of the ratio being`, lu[3]):
  od:
 fi:
od:
od:

print(`The whole web-book took`, time()-t0, `seconds to generate.`):
end:

SymbolicDegree:=proc(term);
## for a term of the form -x^p or x^p, returns p.
op(-1,term/TermCoeff(term))
end proc:

TermCoeff:=proc(term);
## for a term of the form a*x^p or x^p, returns a or 1 respectively.
if type(term,`*`) then op(1,term) else 1 fi
end proc:

pmngf:=proc(m,x) local i;
(-1)^m/mul( numtheory:-cyclotomic(i,x)^(floor(m/i) )     ,i=1..m)
end proc:


Chkgf:=proc(N,x,h,k);
evalc(
taylor( (x-RootOf1(h,k))^(floor(N/k)) * pmngf(N,x),
   x=RootOf1(h,k), N)
)
end proc:


SeriesExpansionCoeffs:=proc(h,k,r,n,s) local len,ex,i,x,signs,powers;
  ex:=expand(mul( (1-x^(k*n-k+i+r)), i=1..k ));
  ex:=simplify(ex,power,symbolic) - 1;
  signs:= [ seq(TermCoeff(op(i,ex) ), i=1..nops(ex)) ];
  powers:=[ seq(SymbolicDegree(op(i,ex)), i=1..nops(ex))];
  if k=1 then signs:=[-1]; powers:=[n] fi:
  if k=1 then len:=1 else len:=nops(ex) fi:
  add( op(i,signs)*binomial( op(i,powers),s )
         * evalc( RootOf1(h,k)^(subs(n=0,op(i,powers))-s mod k) )
     ,i=1..len )
end proc:

RecEqn:= proc(h,k,r,n,j,C) local eq,rr;
  eq:=
   add(  expand(SeriesExpansionCoeffs(h,k,r,n,rr)) *C[h,k,n-j+rr,k*n+r],
     rr=1..j) = C[h,k,n-j,k*n-k+r];
  if j=1 then eq:= simplify(eq) fi:
  eq
end proc:

## Finds a formula for the Rademacher coefficient C[h,k,n-R](k*n+r)
ChkFormula:=proc(h,k,r,n,R)  option remember:
op(R+1, ChkFormulaList(h,k,r,n,R+1) )
end proc:

## Returns a list a formulas for C[h,k,n-j](k*n+r) for j=0..R-1
ChkFormulaList:=proc(h,k,r,n,R) local eq,init,adjust,j,rcf,y,C, x;
option remember:
if gcd(h,k)>1 then RETURN(`ERROR: h and k must be relatively prime`) fi:
if r>= k then RETURN(`ERROR: r must be in the range 0,1,2,...,k-1`) fi:
rcf:=[ ];
for j from 1 to R do
  eq:= RecEqn(h,k,r,n,j,C);
  adjust:={seq( C[h,k,n-y+1,k*n+r]=op(y,rcf),y=1..j-1)};
  adjust:={op(adjust), C[h,k,n-j,k*n-k+r]=f(n-1), C[h,k,n-j+1,k*n+r]=f(n)};
  eq:= subs(adjust,eq);
#  init:= f(j) = coeff(Chkgf(k*j+r,x,h,k),x-RootOf1(h,k),j-1);
  init:= f(j) = simplify(ChklN(h,k,j,k*j+r)):
  rcf:=[op(rcf), rsolve({eq,init}, f(n)) ]
od;
rcf
end proc:

RootOf1:=proc(h,k); cos(2*Pi*h/k) + I*sin(2*Pi*h/k) end proc:


## Automatically generates a paper which states and proves
## formulas for Rademacher coefficients C[h,k, n-j](k*n+r)
## for j from 0 to R-1
HansTopDownAutoPaper:=proc(h,k,r,n,C,R) local t0,i,results,j,
G, N, Theorem,       l, x;
t0:=time();
print(` TOP DOWN FORMULAS FOR SOME RADEMACHER COEFFICIENTS `):
print(` by Tapuah I. Sefer, Georgia Southern University `):
print(` Dedicated to my mentor, Professor Shalosh B. Ekhad `):
print(` `):
print(`Abstract: In this paper, I will prove formulas for several`):
print(`Rademacher coefficients.`):
print(` `):
print(`Definition: Let`,C[h,k,l](N) ,` be the numerator of the term `):
print(`whose denominator is, ` (x-exp(2*Pi*h*I/k))^l, `in the partial fraction`):
print(`decomposition of`, Product( 1/(1-x^j), j=1..N)):
results:=ChkFormulaList(h,k,r,n,R);
print(`Definition: Let`, G[N]=(x-RootOf1(h,k))^(floor(N/k))/Product(1-x^j, j=1..N));
print(`Thus G satsifies the following recurrence (R):`);
print( mul(1-x^(j+r), j=1..k) * G[k*n+r] = (x-RootOf1(h,k)) * G[k*n-k+r] ):
print(` Expand both sides of (R) as a Taylor series about x=`, RootOf1(h,k) ):

for i from 1 to R do
  print(`-------------------------------------------------------------------`):
  print(`Theorem`, i):
  print(C[h,k,n+1-i](k*n+r)= op(i,results)):
  print(`Proof :`):
  print(`Comparing coefficients of`, (x-RootOf1(h,k))^i, `on both sides of (R)` ):
  print(`we obtain the recurrence (R1) :`):
  print(RecEqn(h,k,r,n,i,C) ):
  print(`Solve this recurrence (R1) with the initial condition`,
    C[h,k,i-i+1](k*i+r) = coeff(Chkgf(k*i+r,x,h,k),x-RootOf1(h,k),i-1)):
  if i>1 then print(`taking into account the results of Theorems`,
     seq(j, j=1..i-1) ) fi:
  print(`and the result follows.`)
od:
print(`This paper was produced in`, time()-t0, `seconds of CPU time.`):
end proc:












################C01Table

#OSC01(v,N): applies one step in the recurrence with vectors for
#CNi
OSC01:=proc(v,N) local r,a,kha,v1:
v1:=[]:

for r from 0 to nops(v)-1 do
 kha:=-v[r+1]/N-add(binomial(N,a+1)/N*v1[r-a+1],a=1..r):
 v1:=[op(v1),kha]:
od:

v1:

end:




#C01Table(K,L): A table of C_{01l}(N) for l from 1 to L
#and N from 1 to K, it is given as a list of lists
#such that T[l][N] is C[0,1,l](N). For example, try:
#C01Table(50,10);
C01Table:=proc(K,L) local v,T,N,j:
v:=[1,0$K]:
T[1]:=[]:
for j from 2 to L do
 T[j]:=[]:
od:

for N from 1 to K do
v:=OSC01(v,N ):
T[1]:=[op(T[1]),v[N]]:
 for j from 2 to L do
  if N-j<=-1 then
 T[j]:=[op(T[j]),0]:
  else
 T[j]:=[op(T[j]),v[N-j+1]]:
 fi:
 od:
od:
[seq(T[j],j=1..L)]:
end:


#OSC01f(v,N,Dig): applies one step in the recurrence with vectors for
#CNi
OSC01f:=proc(v,N,Dig) local r,a,kha,v1:
v1:=[]:
Digits:=Dig:
for r from 0 to nops(v)-1 do
 kha:=-v[r+1]/N-add(binomial(N,a+1)/N*v1[r-a+1],a=1..r):
 kha:=evalf(kha):
 v1:=[op(v1),kha]:
od:

v1:

end:




#C01TableF(K,L,Dig): A table of C_{01l}(N) for l from 1 to L
#and N from 1 to K, it is given as a list of lists
#such that T[l][N] is C[0,1,l](N)
C01TableF:=proc(K,L,Dig) local v,T,N,j:
v:=[1.,0$K]:
T[1]:=[]:
for j from 2 to L do
 T[j]:=[]:
od:

for N from 1 to K do
v:=OSC01f(v,N,Dig ):
T[1]:=[op(T[1]),v[N]]:
 for j from 2 to L do
  if N-j<=-1 then
 T[j]:=[op(T[j]),0]:
  else
 T[j]:=[op(T[j]),v[N-j+1]]:
 fi:
 od:
od:
[seq(T[j],j=1..L)]:
end:




################End C01Table



#AsyB:=proc(L,n): Given a list of positive numbers, tries
#to find an expression for them as C*mu^n*n^theta
#it returns the approx. conjenctured asymptotics 
#in the form  C*mu^n*n^theta. 
#For example, try:
#AsyB([seq(5*2^i*i^2,i=1..10)],n):
AsyB:=proc(L,n) local L1,s,eq,var,C1,mu1,t1,n1,i:

if nops(L)<6 then
 print(`List too short`):
 RETURN(FAIL):
fi:
if nops({seq(sign(L[i]),i=nops(L)-5..nops(L))})<>1 then
 RETURN(FAIL):
fi:

if L[1]<0 then
L1:=[seq(-L[i],i=1..nops(L))]:
s:=-1:
else
L1:=L:
s:=1:
fi:

eq:=evalf({seq(log(L1[n1])-C1-n1*mu1-t1*log(n1),n1=nops(L1)-2..nops(L1))}):
var:={C1,mu1,t1}:
var:=solve(eq,var):
C1:=subs(var,C1):
C1:=exp(C1):
if abs(coeff(C1,I,1))>10^(-8) then
  RETURN(FAIL):
fi:
C1:=evalf(coeff(C1,I,0),10):
mu1:=evalf(subs(var,mu1),10):
t1:=evalf(subs(var,t1),10):
s*C1*(exp(mu1))^n*n^t1:
end:



#Ana1(L): analyzes the pattern of differences between peaks and valleys
#of the sequence L, for example, try 
#Ana1(HANS700[1]);
Ana1:=proc(L) local lu,i:
lu:=PaV(L):
lu:=[seq(lu[i][1],i=1..nops(lu))]:
[seq(lu[i+1]-lu[i],i=1..nops(lu)-1)]:
end:

#Ana(L,n): inputs a sequence of real numbers L, and outputs
#conjenctured asymptotics in terms of a list of length
#P, where P is the "period" and the ith-term of the
#output list gives an estimate for the asymptotic behavior of
#L[P*n+i], for i from 1 to P-1, and the P-th term an
#for L[P*n]. 
#It also returns the list [max,min] where
#max is the integer mod P where the local maxima occur
#and min is the integer mod P where the local minima occur
#For example, try:
#read HANS800: Ana(HANS800[1],n);
Ana:=proc(L,n) local lu,i,P,L1,n1,vu,gu,lu1,ru1,ru2,ru:
L1:=evalf(L):
lu:=PaV(L1):
lu1:=[seq(lu[i][1],i=1..nops(lu))]:
lu1:=[seq(lu1[i+1]-lu1[i],i=1..nops(lu)-1)]:

if nops(lu1)<6 then
print(`Not long enough`):
RETURN(FAIL):
fi:

P:={op(nops(lu1)-4..nops(lu1),lu1)}:

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

ru1:=lu[nops(lu)-1][1]:
ru2:=lu[nops(lu)-2][1]:


P:=2*P[1]:

if lu[nops(lu)-1][2]>0 then
 ru:=[ru1 mod P, ru2 mod P]:
else
  ru:=[ru2 mod P, ru1 mod P]:

fi:

gu:=[]:
for i from 1 to P-1 do
lu:=[seq(L1[P*n1+i],n1=1..trunc((nops(L1)-i)/P))]:

vu:=AsyB(lu,n):
gu:=[op(gu),vu]:

od:



lu:=[seq(L1[P*n1],n1=1..trunc(nops(L1)/P))]:
vu:=AsyB(lu,n):
gu:=[op(gu),vu]:


gu,ru:


end:



#C01TableMinus(K,L): A table of C_{01l}(N) for l from 0 to -L
#and N from 1 to K, it is given as a list of lists
#such that T[l+1][N] is C[0,1,-l](N). For example, try:
#C01TableMinus(50,10);
C01TableMinus:=proc(K,L) local v,T,N,j:
v:=[1,0$(K+L)]:
for j from 0 to L do
 T[-j]:=[]:
od:

for N from 1 to K do
v:=OSC01(v,N ):
T[0]:=[op(T[0]),v[N+1]]:
 for j from 1 to L do
  T[-j]:=[op(T[-j]),v[N+j]]:
 od:
od:
[seq(T[-j],j=0..L)]:
end:


#SipurAsy01(n,N): The guessed asymptotics of C01l(32n+r) for l from 1 to 15 and r from 0 to 31.
#Type:
#SipurAsy01(n,N):
SipurAsy01:=proc(n) local gu,l,r,ma,mi,P:

read HANSC01f:

print(`Conjectured Asymptotic Formulas for the Rademacher Numbers C[0,1,l](N) for l from 1 to 15`):
print(``):
print(`By Shalosh B. Ekhad`):
print(``):
for l from 1 to 15 do
print(`-------------------------------------------------------`):
print(`Doing l=`, l):
gu:=Ana(C01f[l],n):

if gu=FAIL then
 next:
fi:

P:=nops(gu[1]):

ma:=gu[2][1]:
mi:=gu[2][2]:
gu:=gu[1]:


print(`The local maxima of`, C[0,1,l](N) , `seem to be  when N is congruent to`, ma, `modulo `, P):
if ma<>0 then
print(`The asymptotic behavior of `, C[0,1,l](P*n+ma), `is approximately`):
print(gu[ma]):
else
print(`The asymptotic behavior of `, C[0,1,l](P*n+ma), `is approximately`):
print(gu[P]):
print(`and in Maple input format`):
lprint(gu[P]):
fi:

print(`The local minima of`, C[0,1,l](N) , `seem to be  when N is congruent to`, mi, `modulo `, P):
if mi<>0 then
print(`The asymptotic behavior of `, C[0,1,l](P*n+mi), `is approximately`):
print(gu[mi]):
else
print(`The asymptotic behavior of `, C[0,1,l](P*n+mi), `is approximately`):
print(gu[P]):
print(`and in Maple input format`):
lprint(gu[P]):
fi:


od :

print(`We will now list conjectured approximate asymptotic formulas for all residue classes mod`, P):


for l from 1 to 10 do
print(`-------------------------------------------------------`):
print(`Right now l=`, l):
gu:=Ana(C01f[l],n):

if gu=FAIL then
 next:
fi:

P:=nops(gu[1]):


gu:=gu[1]:

if gu[P]<>FAIL then
print(`The asymptotic behavior of `, C[0,1,l](P*n), `is approximately`):
print(gu[P]):
print(`and in Maple input format`):
lprint(gu[P]):
fi:

for r from 1 to P-1 do
if gu[r]<>FAIL then
print(`The asymptotic behavior of `, C[0,1,l](P*n+r), `is approximately`):
print(gu[r]):
print(`and in Maple input format`):
lprint(gu[r]):
fi:
od :

od:

end:




#SipurAsy011(n): The guessed asymptotics of C011(32n+r) for r from 0 to 31.
#Type:
#SipurAsy011(n):
SipurAsy011:=proc(n) local gu,r,ma,mi,P,lu:

read oHANS12:

lu:=[seq(E011[i]/(2*i)!,i=1..nops(E011))]:
print(`Conjectured Asymptotic Formulas for the Rademacher Numbers C[0,1,1](N)`):
print(``):
print(`By Shalosh B. Ekhad`):
print(``):
gu:=Ana(lu,n):

if gu=FAIL then
 next:
fi:

P:=nops(gu[1]):

ma:=gu[2][1]:
mi:=gu[2][2]:
gu:=gu[1]:


print(`The local maxima of`, C[0,1,1](N) , `seem to be  when N is congruent to`, ma, `modulo `, P):
if ma<>0 then
print(`The asymptotic behavior of `, C[0,1,1](P*n+ma), `is approximately`):
print(gu[ma]):
else
print(`The asymptotic behavior of `, C[0,1,1](P*n+ma), `is approximately`):
print(gu[P]):
print(`and in Maple input format`):
lprint(gu[P]):
fi:

print(`The local minima of`, C[0,1,1](N) , `seem to be  when N is congruent to`, mi, `modulo `, P):
if mi<>0 then
print(`The asymptotic behavior of `, C[0,1,1](P*n+mi), `is approximately`):
print(gu[mi]):
else
print(`The asymptotic behavior of `, C[0,1,1](P*n+mi), `is approximately`):
print(gu[P]):
print(`and in Maple input format`):
lprint(gu[P]):
fi:




print(`We will now list conjectured approximate asymptotic formulas for all residue classes mod`, P):

print(`-------------------------------------------------------`):
gu:=Ana(lu,n):

if gu=FAIL then
 next:
fi:

P:=nops(gu[1]):


gu:=gu[1]:

if gu[P]<>FAIL then
print(`The asymptotic behavior of `, C[0,1,1](P*n), `is approximately`):
print(gu[P]):
print(`and in Maple input format`):
lprint(gu[P]):
fi:

for r from 1 to P-1 do
if gu[r]<>FAIL then
print(`The asymptotic behavior of `, C[0,1,1](P*n+r), `is approximately`):
print(gu[r]):
print(`and in Maple input format`):
lprint(gu[r]):
fi:
od :



end:



###start  addition of March 13, 2012


#E01(j,N): The Rademacher coeff. C_{01j}(N)
#computed by the recurrence
E01:=proc(j,N):
ENi(N,N-j)/(2*N)!:
end:

#ENi(N,i): the coeff. of z^i in z^(N-1)/((z+1)^2-1))/((z+1)^3-1)...
#by the recurrence. For example, try:
#ENi(5,4);
ENi:=proc(N,i) local r:
option remember:

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

if N<1 then
 RETURN(FAIL):
fi:

if N=1 then
 if i=0 then
   RETURN(-2):
  else
 RETURN(0):
 fi:
fi:

-ENi(N-1,i)*(4*N-2)-add(binomial(N,r+1)/N*ENi(N,i-r),r=1..N-1):

end:




#OS(L,N): going from the sequence [ENi(N-1,0),..., ENi(N-1,R)]
#to [ENi(N,0),..., ENi(N,M)]
#by the recurrence. For example, try:
#OS([-2,0,0],2);
OS:=proc(L,N) local gu,r,a:

gu:=[]:

for r from 0 to nops(L)-1 do
gu:=[op(gu),-(4*N-2)*L[r+1]-add(binomial(N,a+1)/N*gu[r-a+1],a=1..r)]:
od:

end:


#E01s(N): the first N terms of the sequence ENi(i,i-1)
E01s:=proc(N) local gu,mu,i:
mu:=[-2,0$N]:
gu:=[]:

for i from 2 to N+1 do
gu:=[op(gu),mu[i-1]]:
mu:=OS(mu,i):
od:

gu:
end:

###end addition of March 13, 2012