######################################################################
##EvenChange.txt: Save this file as  EvenChange.txt                  #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read EvenChange.txt<Enter>                         #
##Then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#DoronZeil at gmail dot com                                          #
######################################################################
 
#Created: Jan. 2019
 
print(`Created:  Jan. 2019`):
print(` This is EvenChange.txt `):
print(`It  accompanies the article `):
print(` In How many ways can I carry a total of n coins in my two pockets, and have the same amount in both pockets?`):
print(`by Shalosh B. Ekhad and Doron Zeilberger`):
print(` available from Doron Zeilberger's website`):
print(``):
print(`Dedicated to the memory of Gert Almkvist (1934-2018)`):
print(``):

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

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

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

with(combinat):

ezra1:=proc()

if args=NULL then
 print(` The supporting procedures are: ChangeInfoVV, fPABzt, IsPos, SeqPABd`):
 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are: ChangeInfo, ChangeInfoV, EvalQPS, GertPaper, GfPAB, GfPABinfo, GfPABinfoV, GFtoQPS, SeqPAB `):
 print(` `):


elif nops([args])=1 and op(1,[args])=ChangeInfo then
print(`ChangeInfo(A,t,n,N0): inputs a set of positive integes A, a continuous variable t, a discrete variable n and a large positive integer N0 outputs`):
print(`a list of`):
print(`(i) the rational function in t whose coefficient of t^n is the number of ways of having n coins in your pockets, where the denominations are`):
print(`from the set A and the total monetary amounts in the two pockets is the same`):
pritn(`(ii) The sum of quasi-polynomial in n of that quantity`):
print(`The actual value for n=N0. Try:`):
print(`ChangeInfo({1,5,10,25},t,n,10^100);`):

elif nops([args])=1 and op(1,[args])=ChangeInfoV then
print(`ChangeInfoV(A,t,n,N0): a verbose form of ChangeInfo(A,t,n,N0) (q.v.). Try: `):
print(`ChangeInfoV({1,5,10,25},t,n,10^100);`):

elif nops([args])=1 and op(1,[args])=ChangeInfoVV then
print(`ChangeInfoVV(A,t,n,N0): a very verbose form of ChangeInfo(A,t,n,N0) (q.v.). Try: `):
print(`ChangeInfoVV({1,5,10,25},t,n,10^100);`):

elif nops([args])=1 and op(1,[args])=EvalQPS then
print(`EvalQPS(L,n,n0): evaluates the sum of the `):
print(`quasi-polynomials in the list L of n at n0`):
print(`For example try:`):
print(`EvalQPS([[n^2],[2*n^3,n^3],[1,2,6]],n,5);`):

elif nops([args])=1 and op(1,[args])=fPABzt then
print(`fPABzt(P,A,B,z,t):P/mul(1-t*z^a,a in A)/mul(1-t/z^b,b in B):`):
print(`fPABzt((1+z)^2/2/z,{1,3},{1,3},z,t);`):

elif nops([args])=1 and op(1,[args])=GertPaper then
print(`GertPaper(N,N0):  Inputs a positive integer N and a very large positive integer N0 and outputs an article of interest to`):
print(`Gert Almkvist (1934-2018). Try:`):
print(`GertPaper(6,10^100);`):

elif nops([args])=1 and op(1,[args])=GfPAB then
print(`GfPAB(P,z,t,A,B): inputs a Laurent polynomial P of the variable z and outputs the constant term of`):
print(`P(z)/(mul(1-t*z^a, a in A)*mul(1-t/z^b, b in P). Try:`):
print(`GfPAB((1+z)^2/(2*z),z,t,{1,3},{1,3});`):

elif nops([args])=1 and op(1,[args])=GfPABinfo then
print(`GfPABinfo(P,z,t,A,B,n,N0): inputs (i) a Laurent polynomial in P, in the variable z`):
print(` (ii) a variable z `):
print(` (iii) sets of positive integers A and B `):
print(` (iv) a (continuous) variable t `):
print(` (v) a (discrete) variable n `):
print(`(vi) a large positive integer N0`):
print(` Outputs `):
print(` (i) the rational function in t that is the constant term, in z, of P/Product(1-z^a*t, a in A)/Product(1-z^(-b)*t, b in B) `):
print(`(ii) the first 31 terms in the Maclaurin expansion in t starting at n=0`):
print(`(iii) The sum of quasi-polynomial in n of that quantity`):
print(`(iv) The actual value for n=N0. Try:`):
print(`To get OEIS sequence A1971, type: `):
print(`GfPABinfo((1+z)^2/2/z,z,t,{1,3},{1,3},n,10^100);`):
print(`To get OEIS sequence A1975, type: `):
print(`GfPABinfo((1+z)^2/2/z,z,t,{1,3,5},{1,3,5},n,10^100);`):
print(`To get OEIS sequence A1979, type: `):
print(`GfPABinfo((1+z)^2/2/z,z,t,{1,3,5,7},{1,3,5,7},n,10^100);`):
print(`The following is not yet (Jan. 22, 2019) in the OEIS`):
print(`GfPABinfo((1+z)^2/2/z,z,t,{1,3,5,7,9},{1,3,5,7,9},n,10^100);`):

elif nops([args])=1 and op(1,[args])=GfPABinfoV then
print(`GfPABinfoV(P,z,t,A,B,n,N0): verbose form of GfPABinfo(P,z,t,A,B,n,N0) (q.v.). Try: `):
print(`To get OEIS sequence A1971, type: `):
print(`GfPABinfoV((1+z)^2/2/z,z,t,{1,3},{1,3},n,10^100);`):
print(`To get OEIS sequence A1975, type: `):
print(`GfPABinfoV((1+z)^2/2/z,z,t,{1,3,5},{1,3,5},n,10^100);`):
print(`To get OEIS sequence A1979, type: `):
print(`GfPABinfoV((1+z)^2/2/z,z,t,{1,3,5,7},{1,3,5,7},n,10^100);`):
print(`The following is not yet (Jan. 22, 2019) in the OEIS`):
print(`GfPABinfoV((1+z)^2/2/z,z,t,{1,3,5,7,9},{1,3,5,7,9},n,10^100);`):

elif nops([args])=1 and op(1,[args])=GFtoQPS then
print(`GFtoQPS(Gu,q,n,MaxK): An explicit expression, in terms of n,  as a sum of quasi-polynomials  in n for the coefficient of q^n in the `):
print(`rational function Gu of q that is whose denominator has roots that are roots of unity.`):
print(`MaxK is a guessing parameter. `):
print(` For example, Try: `):
print(`GFtoQPS(1/(1-t)/(1-t^2),t,n,400);`):

elif nops([args])=1 and op(1,[args])=IsPos then
print(`IsPos(f,z,t): inputs a factored polynomial f with one term of the form 1-z^a*t or t-z^a outputs true if the former. Try:`):
print(`IsPos((1-z*t^3)*(1+t),z,t);`):

elif nops([args])=1 and op(1,[args])=SeqPAB then
print(`The first N+1 terms of the constant term of the Taylor coefficients of  fPABzt(P,A,B,z,t)  done via the generating functin. Try:`):
print(`SeqPAB((1+z)^2/2/z,z,{1,3},{1,3},20);`):

elif nops([args])=1 and op(1,[args])=SeqPABd then
print(` The first N+1 terms of the constant term of the Taylor coefficients of  fPABzt(P,A,B,z,t)  done Directly. Try:`):
print(`SeqPABd((1+z)^2/2/z,z,{1,3},{1,3},20);`):


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



fnzt:=proc(n,z,t) local i:
(1+z)^2/2/z/mul(1-t*z^(n-2*i),i=0..n):
end:

#fPABzt(P,A,B,z,t):P/mul(1-t*z^a,a in A)/mul(1-t/z^b,b in B):
#fPABzt((1+z)^2/2/z,{1,3},{1,3},z,t);
fPABzt:=proc(P,A,B,z,t) local a,b:
P/mul(1-t*z^a,a in A)/mul(1-t/z^b,b in B):
end:

#The first N+1 terms of the Taylor coefficients of fnzt(n,z,t) done directly for n odd.
#Try: SeqGd(3,20);
SeqGd:=proc(n,N) local z,t,gu,i:

if not type(n,integer) and n>0 and n mod 2=1 then
 RETURN(FAIL):
fi:

gu:=fnzt(n,z,t):
gu:=taylor(gu,t=0, N+1):
[seq(coeff(expand(coeff(gu,t,i)),z,0),i=0..N)]:
end:

#The first N+1 terms of the Taylor coefficients of fnzt(n,z,t) done via
#fPABzt((1+z)^2/2/z,{1,3,..,n},{-1,-3, ..., -n},z,t); `):
SeqG:=proc(n,N) local z,t,gu,i:
if not type(n,integer) and n>0 and n mod 2=1 then
 RETURN(FAIL):
fi:


gu:=GfPAB((1+z)^2/2/z,z,t,{seq(2*i-1,i=1..(n+1)/2)}, {seq(2*i-1,i=1..(n+1)/2)}):
gu:=taylor(gu,t=0, N+1):
[seq(coeff(gu,t,i),i=0..N)]:
end:


#IsPos(f,z,t): inputs a factored polynomial f with one term of the form 1-z^a*t or t-z^a outputs true if the former.
#Try: IsPos((1-z*t^3)*(1+t),z,t);
IsPos:=proc(f,z,t) local i,f1:

for i from 1 to nops(f) do
f1:=op(i,f):
 
 if expand(diff(f1,z))<>0 then
   if not type(coeff(f1,t,1),integer) then
    RETURN(true):
   fi:
 fi:
od:
false:
end:



#GfPAB(P,z,t,A,B): inputs a Laurent polynomial P of the variable z and outputs the constant term of
#P(z)/(mul(1-t*z^a, a in A)*mul(1-t/z^b, b in P). Try:
#GfPAB((1+z)^2/(2*z),z,t,{1,3},{1,3});
GfPAB:=proc(P,z,t,A,B) local lu,a,b,gu,i:

lu:=P/mul(1-t*z^a,a in A)/mul(1-t/z^b,b in B):

lu:=convert(lu,parfrac,z):

gu:=0:

for i from 1 to nops(lu) do
 if IsPos(denom(op(i,lu)),z,t)  then
  gu:=gu+op(i,lu):
 fi:
od:

factor(normal(subs(z=0,gu))):

end:

 





#The first N+1 terms of the constant term of the Taylor coefficients of  fPABzt(P,A,B,z,t)  done directly 
#Try: SeqPABd((1+z)^2/2/z,z,{1,3},{1,3},20);
SeqPABd:=proc(P,z,A,B,N) local t,gu,i:

gu:=fPABzt(P,A,B,z,t):

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

#The first N+1 terms of the constant term of the Taylor coefficients of  fPABzt(P,A,B,z,t)  done via
#Try: SeqPABd((1+z)^2/2/z,z,{1,3},{1,3},20);
SeqPAB:=proc(P,z,A,B,N) local t,gu,i:

gu:=GfPAB(P,z,t,A,B):
gu:=taylor(gu,t=0, N+1):
[seq(coeff(expand(coeff(gu,t,i)),z,0),i=0..N)]:
end:


###Start from PARTITIONS

#EvalQP(Q,n,n0): evaluates the quasi-polynomial Q of n at n0
#For example try:
#EvalQP([2*n^3,n^3],n,5);
EvalQP:=proc(Q,n,n0) local i:

i:=n0 mod nops(Q):

if i=0 then i:=nops(Q): fi:

subs(n=n0,Q[i]):
end:



#EvalQPS(L,n,n0): evaluates the sum of the 
#quasi-polynomials in the list L of n at n0
#For example try:
#EvalQPS([[n^2],[2*n^3,n^3],[1,2,6]],n,5);
EvalQPS:=proc(L,n,n0) local i:
add(EvalQP(L[i],n,n0),i=1..nops(L)):
end:

#KhaberLists(S): given a set of lists of the same size, adds them
#up, returns FAIL otherwise. For example, type:
#KhaberLists({[1,2],[2,4]});
KhaberLists:=proc(S) local s,kv,k,i1:

kv:={seq(nops(s), s in S)}:

if nops(kv)<>1 then
 RETURN(FAIL):
fi:
k:=kv[1]:

[seq(sort(add(s[i1],s in S)),i1=1..k)]:

end:

#GuessPol1(L,d,n): guesses a polynomial of degree d in n for
# the list L, such that P(i)=L[i+1] for i=1..nops(L)
#For example, try: 
#GuessPol1([seq(i,i=1..10)],1,n);
GuessPol1:=proc(L,d,n) local P,i,a,eq,var:
if d>nops(L)-2 then
 ERROR(`the list is too small`):
fi:

P:=add(a[i]*n^i,i=0..d):
var:={seq(a[i],i=0..d)}:
eq:={seq(subs(n=i,P)-L[i],i=1..nops(L))}:

var:=solve(eq,var):

if var=NULL then
 RETURN(FAIL):
fi:

subs(var,P):

end:


#GuessQP11(L,d,r,n): Inputs a list L and outputs a list of polynomials
#of degree d
#of length r, let's call it M, such that L[i+m*r]=M[i](i+m*r) for
#i=1..r, For example try:
#GuessQP11:=proc([1,-1,1,-1,1,-1,1,-1,1,-1],0,2,n);
GuessQP11:=proc(L,d,r,n) local M,i,L1,gu,m:
M:=[]:
for i from 1 to r do
L1:=[seq(L[i+(m-1)*r],m=1..trunc((nops(L)-i)/r)+1)]:
gu:=GuessPol1(L1,d,m):
if gu=FAIL then
  RETURN(FAIL):
fi:
gu:=expand(subs(m=(n-i)/r+1,gu)):
M:=[op(M),gu]:
end:
M:
end:





#GuessQP2(L,d,n): Inputs a list L and outputs a list of polynomials
#of degree d, the list being
#length r, for the smallest possible
#r, let's call it M, such that L[i+m*r]=M[i](i+m*r) for
#i=1..r, 
#For example try:
#GuessQP2:=proc([1,-1,1,-1,1,-1,1,-1,1,-1],0,n);
GuessQP2:=proc(L,d,n) local r,gu:

for r from 1 to trunc(nops(L)/(d+2)) do
gu:=GuessQP11(L,d,r,n):
if gu<>FAIL then
 RETURN(gu):
fi:
od:
FAIL:
end:


#FindExp(P): Given an expression of the form
#P=C*pol(q)^i finds i, for example, try:
#FindExp(6*(1+q+q^2)^3);
FindExp:=proc(P) local P1:

if type(P,`^`) then
 RETURN(op(2,P)):
fi:

if type(P,`*`) then
  P1:=op(2,P):
  if type(P1,`^`) then
    RETURN(op(2,P1)):
  else
    RETURN(1):
  fi:
fi:
1:
end:

#pmn(m,n): An explicit expression  for p_m(n): the
#number of integer partitions on n with largest part <=m
#(equivalently, the number of parts<=m) as a sum
#of quasi-polynomials. 
#For example, try:
#pmn(3,n);
pmn:=proc(m,n) local gu,q,i,lu1,d,H,i1,gu1,T,makh,Gu,mu,K:
option remember:
K:=2*m:
if m=1 then
 RETURN([[1]]):
fi:

Gu:=mul(1/(1-q^i),i=1..m):
gu:=convert(Gu,parfrac):

H:={}:

for i from 1 to nops(gu) do
 gu1:=op(i,gu):
 d:=FindExp(denom(gu1))-1:
 lu1:=taylor(gu1,q=0,K+1):
 lu1:=[seq(coeff(lu1,q,i1),i1=1..K)]:
 lu1:= GuessQP2(lu1,d,n):
   if lu1=FAIL then
     RETURN(FAIL):
   else
    H:=H union {lu1}:
   fi:
od:

makh:={seq(nops(H[i1]),i1=1..nops(H))}:
makh:=sort(convert(makh,list)):

for i1 from 1 to nops(makh) do
 T[makh[i1]]:= {}:
od:

for i1 from 1 to nops(H) do
T[nops(H[i1])]:=T[nops(H[i1])] union {H[i1]}:
od:

H:=[seq(T[makh[i1]],i1=1..nops(makh))]:

H:=[seq(KhaberLists(H[i1]),i1=1..nops(H))]:

gu:=taylor(Gu,q=0,2*K+1+2*m):
gu:=[seq(coeff(gu,q,i),i=1..2*K+2*m)]:
mu:=[seq(EvalQPS(H,n,i),i=1..2*K+2*m)]:
if gu<>mu then
print(H, `didn't work out, too bad`):
RETURN(FAIL):
fi:
H:
end:




#GFtoQPS1(Gu,q,n,K): An explicit expression, in terms of n,  as a sum of quasi-polynomials  in n for the coefficient of q^n in the rational function Gu of q that is
#whose denominator has roots that are roots of unity.
#For example, Try:
#GFtoQPS1(1/(1-t)/(1-t^2),t,n,40);
GFtoQPS1:=proc(Gu,q,n,K) local gu,i,lu1,d,H,i1,gu1,T,makh,mu:
option remember:
gu:=convert(Gu,parfrac):

H:={}:

for i from 1 to nops(gu) do
 gu1:=op(i,gu):
 d:=FindExp(denom(gu1))-1:
 lu1:=taylor(gu1,q=0,K+1):
 lu1:=[seq(coeff(lu1,q,i1),i1=1..K)]:
 lu1:= GuessQP2(lu1,d,n):
   if lu1=FAIL then
     RETURN(FAIL):
   else
    H:=H union {lu1}:
   fi:
od:

makh:={seq(nops(H[i1]),i1=1..nops(H))}:
makh:=sort(convert(makh,list)):

for i1 from 1 to nops(makh) do
 T[makh[i1]]:= {}:
od:

for i1 from 1 to nops(H) do
T[nops(H[i1])]:=T[nops(H[i1])] union {H[i1]}:
od:

H:=[seq(T[makh[i1]],i1=1..nops(makh))]:

H:=[seq(KhaberLists(H[i1]),i1=1..nops(H))]:

gu:=taylor(Gu,q=0,6*K+1):
gu:=[seq(coeff(gu,q,i),i=1..6*K)]:
mu:=[seq(EvalQPS(H,n,i),i=1..6*K)]:
if gu<>mu then
print(H, `didn't work out, too bad`):
RETURN(FAIL):
fi:
H:
end:



#GFtoQPS(Gu,q,n,MaxK): 
#An explicit expression, in terms of n,  as a sum of quasi-polynomials  in n for the coefficient of q^n in the rational function Gu of q that is
#whose denominator has roots that are roots of unity.
#For example, Try:
#GFtoQPS(1/(1-t)/(1-t^2),t,n,1000);
GFtoQPS:=proc(Gu,q,n,MaxK) local gu,K:

for  K from 40 by 20 to MaxK do
 gu:=GFtoQPS1(Gu,q,n,K) :
 if gu<>FAIL then
   RETURN(gu):
 fi:
od:
FAIL:
end:


###end from PARTITIONS######################################################################


#ChangeInfo(A,t,n,N0): inputs a set of positive integes A, a continuous variable t, a discrete variable n and a large positive integer N0 outputs
#a list of
#(i) the rational function in t whose coefficient of t^n is the number of ways of having n coins in your pockets, where the denominations are
#from the set A and the total monetary amounts in the two pockets is the same
#(ii) The sum of quasi-polynomial in n of that quantity
#The actual value for n=N0. Try:
#ChangeInfo({1,5,10,25},t,n,10^100);
ChangeInfo:=proc(A,t,n,N0) local lu,z,mu:
lu:=GfPAB(1,z,t,A,A);
mu:=GFtoQPS(lu,t,n,1000):
if mu<>FAIL then
 RETURN([lu,mu,EvalQPS(mu,n,N0)]):
else
 RETURN([lu]):
fi:

end:




#ChangeInfoVV(A,t,n,N0): A Very verbose form of ChangeInfo(A,t,n,N0) (q.v.)
ChangeInfoVV:=proc(A,t,n,N0) local gu,i,j,P,d:
gu:=ChangeInfo(A,t,n,N0):

print(``):

print(` The  Number of Ways of Having n Coins in Your Two pockets, with denominations in the set`, A, ` in such a way`):
print(`That both pockets carry the same amount `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Theorem: Let a(n) be the n the  number of ways of having n Coins in your two pockets, with denominations in the set`, A, ` in such a way`):
print(`that both pockets carry the same amount , then`):
print(``):
print(Sum(a(n)*t^n, n=0..infinity) = gu[1]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu[1]):
print(``):
print(``):
if nops(gu)>1 then
 print(`Furthermore, a(n)  is a quasi-polynomial given as sum of`, nops(gu[2]), ` quasi-polynomials `):
 print(``):
 print(a(n)=Sum(P[i](n),i=1..nops(gu[2]) ))  :
 print(``):
 print(` where `, seq(P[i](n),i=1..nops(gu[2])), ` are defined as followed `):
print(``):

 print(P[1](n), `is the polynomial `);
 print(``):
print( gu[2][1][1]):
 print(``):
print(`and in Maple notation`):
 print(``):
lprint( gu[2][1][1]):
 print(``):
print(`This is the leading term`):
 print(``):
 d:=degree(gu[2][1][1],n):
 print(``):
  print(`in particular, a(n) , is asymptotic to`):
 print(``):

 print(coeff(gu[2][1][1],n,d)*n^d):
 print(``):
  for i from 2 to nops(gu[2]) do
   print(P[i](n), `is defined as follows`):
   print(``):
   for j from 1 to nops(gu[2][i]) do
    print(`if n is `, j, ` modulo`, nops(gu[2][i]) , `then  it equals `):
    print(``):
    print(gu[2][i][j]):
     print(``):
  print(`and in Maple notation`):
    print(``):
    lprint(gu[2][i][j]):
      print(``):
  od:
 print(``):
 od:

print(``):
print(`That makes it very easy to compute a(n) for large n. In particular, the number of ways of having`, N0, `coinss in  yourtwo pocekts`):
print(`so that both pockets have the same amount is`):
print(``):
lprint(gu[3]):
fi:

print(``):
print(`------------------------------------`):
print(``):
end:



#ChangeInfoV(A,t,n,N0): A Very verbose form of ChangeInfo(A,t,n,N0) (q.v.)
ChangeInfoV:=proc(A,t,n,N0) local gu,i,P,d,t0:
t0:=time():
gu:=ChangeInfo(A,t,n,N0):

print(``):

print(` The  Number of Ways of Having n Coins in Your Two pockets, with denominations in the set`, A, ` in such a way`):
print(`That both pockets carry the same amount `):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Theorem: Let a(n) be the  number of ways of having n Coins in your two pockets, with denominations in the set`, A, ` in such a way`):
print(`that both pockets carry the same amount , then`):
print(``):
print(Sum(a(n)*t^n, n=0..infinity) = gu[1]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu[1]):
print(``):
print(`For the sake of Sloane, here are the first 31 term, starting at n=0`):
print(``):
print(seq(coeff(taylor(gu[1],t=0,31),t,i),i=0..30)):
print(``):
if nops(gu)>1 then
 print(`Furthermore, a(n)  is a quasi-polynomial given as sum of`, nops(gu[2]), ` quasi-polynomials `):
 print(``):
 print(a(n)=Sum(P[i](n),i=1..nops(gu[2]) ))  :
 print(``):
 print(` where `, seq(P[i](n),i=1..nops(gu[2])), ` are defined as followed `):
print(``):

 print(P[1](n), `is the polynomial `);
 print(``):
print( gu[2][1][1]):
 print(``):
print(`and in Maple notation`):
 print(``):
lprint( gu[2][1][1]):
 print(``):
print(`This is the leading term`):
 print(``):
 d:=degree(gu[2][1][1],n):
 print(``):
  print(`in particular, a(n) , is asymptotic to`):
 print(``):

 print(coeff(gu[2][1][1],n,d)*n^d):
 print(``):
  for i from 2 to nops(gu[2]) do
   print(P[i](n), `is defined by the following list whose i-th entry is the expression if n is congruent to i mod`, nops(gu[2][i])):
      print(``):
  print(gu[2][i]):
   print(``):
   print(`and in Maple format`):
  lprint(gu[2][i]):
 print(``):
 od:

print(``):
print(`That makes it very easy to compute a(n) for large n. In particular, the number of ways of having`, N0, `coinss in  yourtwo pocekts`):
print(`so that both pockets have the same amount is`):
print(``):
lprint(gu[3]):
fi:

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





#GfPABinfo(P,z,A,B,t,n,N0): inputs (i) a Laurent polynomial in P, in the variable z
#(ii) a variable z
#(iii) sets of positive integers A and B
#(iv) a (continuous) variable t
#(v) a (discrete) variable n
#(vi) a large positive integer N0
#Outputs
#(i) the rational function in t that is the constant term, in z, of P/Product(1-z^a*t, a in A)/Product(1-z^(-b)*t, b in B)
#(ii) the first 31 terms in the Maclaurin expansion in t starting at n=0
#(iii) The sum of quasi-polynomial in n of that quantity
#(iv) The actual value for n=N0. Try:
#GfPABinfo((1+z)^2/2/z,z,t,{1,3},{1,3},n,10^100);

GfPABinfo:=proc(P,z,t,A,B,n,N0) local lu,mu,ku,i:
lu:=GfPAB(P,z,t,A,B);
ku:=[seq(coeff(taylor(lu,t=0,31),t,i),i=0..30)]:
mu:=GFtoQPS(lu,t,n,1000):
if mu<>FAIL then
 RETURN([lu,ku,mu,EvalQPS(mu,n,N0)]):
else
 RETURN([lu,ku]):
fi:

end:


#GfPABinfoV(P,z,t,A,B,n,N0): verbose form of GfPABinfo(P,z,t,A,B,n,N0) (q.v.)
GfPABinfoV:=proc(P,z,t,A,B,n,N0) local gu,a,b,d,i,Q:
gu:=GfPABinfo(P,z,t,A,B,n,N0):

print(``):

print(` The Constant Term in`, z, `of a certain Rational function in `, z, ` and `, t):
print(``):
print(`By Shalosh B. Ekhad `):
print(``):
print(`Theorem: Let f(t) be the constant term of the rational function in z and t`):
print(``):
print(P/mul(1-z^a*t, a in A)/mul(1-t/z^b, b in B)):
print(``):
print(` We have `):
print(``):
print(f(t) = gu[1]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu[1]):
print(``):
print(`For the sake of Sloane, here are the first 31 term, starting at n=0`):
print(``):
print(gu[2]):
print(``):
if nops(gu)>2 then
 print(`Furthermore, a(n)  is a quasi-polynomial given as sum of`, nops(gu[3]), ` quasi-polynomials `):
 print(``):
 print(a(n)=Sum(Q[i](n),i=1..nops(gu[3]) ))  :
 print(``):
 print(` where `, seq(Q[i](n),i=1..nops(gu[3])), ` are defined as followed `):
print(``):

 print(Q[1](n), `is the polynomial `);
 print(``):
print( gu[3][1][1]):
 print(``):
print(`and in Maple notation`):
 print(``):
lprint( gu[3][1][1]):
 print(``):
print(`This is the leading term`):
 print(``):
 d:=degree(gu[3][1][1],n):
 print(``):
  print(`in particular, a(n) , is asymptotic to`):
 print(``):

 print(coeff(gu[3][1][1],n,d)*n^d):
 print(``):
  for i from 2 to nops(gu[3]) do
   print(Q[i](n), `is defined by the following list whose i-th entry is the expression if n is congruent to i mod`, nops(gu[3][i])):
      print(``):
  print(gu[3][i]):
   print(``):
   print(`and in Maple format`):
  lprint(gu[3][i]):
 print(``):
 od:

print(``):
print(`That makes it very easy to compute a(n) for large n. In particular`, a(N0), `equals `):
print(``):
lprint(gu[4]):
fi:

print(``):
print(`------------------------------------`):
print(``):
end:




#GfPABinfoVch(P,z,t,A,B,n,N0,ch): Like GfPABinfoVch(P,z,t,A,B,n,N0) but with chapter number
GfPABinfoVch:=proc(P,z,t,A,B,n,N0,ch) local gu,a,b,d,i,Q:
gu:=GfPABinfo(P,z,t,A,B,n,N0):

print(``):

print(`Theorem Number`, ch, `: Let f(t) be the constant term of the rational function in z and t`):
print(``):
print(P/mul(1-z^a*t, a in A)/mul(1-t/z^b, b in B)):
print(``):
print(` We have `):
print(``):
print(f(t) = gu[1]):
print(``):
print(`and in Maple notation`):
print(``):
lprint(gu[1]):
print(``):
print(`For the sake of Sloane, here are the first 31 term, starting at n=0`):
print(``):
print(gu[2]):
print(``):
if nops(gu)>2 then
 print(`Furthermore, a(n)  is a quasi-polynomial given as sum of`, nops(gu[3]), ` quasi-polynomials `):
 print(``):
 print(a(n)=Sum(Q[i](n),i=1..nops(gu[3]) ))  :
 print(``):
 print(` where `, seq(Q[i](n),i=1..nops(gu[3])), ` are defined as followed `):
print(``):

 print(Q[1](n), `is the polynomial `);
 print(``):
print( gu[3][1][1]):
 print(``):
print(`and in Maple notation`):
 print(``):
lprint( gu[3][1][1]):
 print(``):
print(`This is the leading term`):
 print(``):
 d:=degree(gu[3][1][1],n):
 print(``):
  print(`in particular, a(n) , is asymptotic to`):
 print(``):

 print(coeff(gu[3][1][1],n,d)*n^d):
 print(``):
  for i from 2 to nops(gu[3]) do
   print(Q[i](n), `is defined by the following list whose i-th entry is the expression if n is congruent to i mod`, nops(gu[3][i])):
      print(``):
  print(gu[3][i]):
   print(``):
   print(`and in Maple format`):
  lprint(gu[3][i]):
 print(``):
 od:

print(``):
print(`That makes it very easy to compute a(n) for large n. In particular`, a(N0), `equals `):
print(``):
lprint(gu[4]):
fi:

print(``):
print(`------------------------------------`):
print(``):
end:



#GertPaper(N,N0):  Inputs a positive integer N and a very large positive integer N0 and outputs a book of interest to
#Gert Almkvist (1934-2018). Try:
#GertPaper(6,10^100);
GertPaper:=proc(N,N0) local psi,ch,n1,n,t,z,t0,i1:
t0:=time():
print(`The first `, N, `rational functions `, psi[n],  `discussed in Gert Almkvist's Beautiful article "Invariants, Mostly Old ones" `):
print(``):
print(`Pacific Journal of Mathematics v. 86 (1980), No. 1, 1-13. `):
print(``):
print(`By Shalosh B. Ekhad`):
print(``):
print(`In a beautiful article of Gert Almkvist (1934-2018) , "Invariants, Mostly Old ones" , Pacific Journal of Mathematics v. 86 (1980), No. 1, 1-13, `):
print(``):
print(`he  investigated certaion rational functions, that he called`, psi[n](t), ` that turned out to have been studied `):
print(`by Faa Di Bruno, Cayley, Sylvester, and Fabian Franklin`):
print(``):
print(`Here we present fully rigorous, computer-generated explicit expressions for this quantities, as well as explicit expressions`):
print(``):
print(`as sum of quasi-polynomials for the coefficients.`):
print(``):
print(`These generating functions turned out to enumerate covariantes of binary quadratic forms of degree n.`):
print(``):
print(`and have been computed by James Joseph Sylvester and his student, Fabian Franklin (of Franklin's involution fame) up to n=10 and n=12`):
print(``):
print(`Here we will do it up to n=`, N):

ch:=0:
for n1 from 2 to N do
 print(``):
 print(` For `, psi[n1](t), `we have the following theorem. `):
 print(``):
 ch:=ch+1:

 if n1 mod 2=0 then
  GfPABinfoVch((1+z)^2/2/z,z,t,{seq(2*i1,i1=0..n1/2)} , {seq(2*i1,i1=1..n1/2)} ,n,N0,ch):

else
GfPABinfoVch((1+z)^2/2/z,z,t,{seq(2*i1+1,i1=0..(n1-1)/2)} , {seq(2*i1+1,i1=0..(n1-1)/2)} ,n,N0,ch):
fi:
print(``):
print(`---------------------------------------------------------`):
od:

print(``):

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

end: