#####################################################################
##FOOL.txt: Save this file as  FOOL.txt     .                        #
## To use it, stay in the                                            #
##same directory, get into Maple (by typing: maple <Enter> )         #
##and then type:  read  FOOL.txt<Enter>      ,                       #
##then follow the instructions given there                           #
##                                                                   #
##Written by Doron Zeilberger, Rutgers University ,                  #
#DoronZeil at gmail dot com                                          #
######################################################################
 
#Created: June 20, 2018
 
print(`Created:  June 20, 2018`):
print(` This is  FOOL.txt `):
print(`It accompanies Doron Zeilberger's plenary talk at the conference "Combinatory Analysis 2018" in honor of George Andrews' forthcoming`):
print(`80th birthday` ):
print(``):
print(`It the packages that created the book `):
print(` A Fool Can Ask More Questions than a Wise Person Can Answer`):
print(`by Shalosh B. Ekhad `):
print(`dedicated to George Andrews (b. Dec. 4, 1938) `):

print(` It is available from Zeilberger's website`):
print(``):
 print(`http://sites.math.rutgers.edu/~zeilberg/FoolBook.html  .`):
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: ChIpre, ChIIpre, ChIIIpre, ChIVpre, RandBCS,  RandPol, RandPrime, RandqBCS  `):

 print(``):

else
ezra(args):
fi:

end:

ezra:=proc()

if args=NULL then
 print(`The main procedures are:ChI, ChII, ChIII, ChIV  `):
 print(` `):



elif nops([args])=1 and op(1,[args])=ChI then
print(`ChI(N,k): a Problem Chapter with k challenging integer factorization problems, to pactorize 2*N-digits positive integers `):
print(` For example Try: `):
print(` ChI(50,10); `):

elif nops([args])=1 and op(1,[args])=ChII then
print(`ChII(x,K1,D1,N1,R1,k): Outputs Chapter II k challenging problems about proving that a given polynomial is non-negative, i.e. that its global`):
print(`minimum is larger or equal to 0. `):
print(`where the answer is by displaying the polynomial as a sum of squares and there are at least 3 and at most R1 terms. Try:`):
print(`ChII([x,y,z,w],10,10,10000,10,10);`):


elif nops([args])=1 and op(1,[args])=ChIII then
print(`ChIII(R1,L1,KAMA):A chapter with KAMA challenging problems about proving that a given `):
print(`sequence defined by a linear recurrence equation with polynomial coefficients with certain initial conditions all produce integers`):
print(`ChIII(5,10,10);`):

elif nops([args])=1 and op(1,[args])=ChIV then
print(`ChIV(R1,R2,R3,MaxOrder,KAMA): KAMA challenging problems for George Andrews  `):
print(`inspired by his beautiful treatment of the L-M-W conjecture. Try: `):
print(`ChIV(n,j,N,2,3,5,3,10);`):

elif nops([args])=1 and op(1,[args])=ChIpre then
print(`ChIpre(N,k): a list of k random integer factorization challented of the form [n,[p,q]] where n=p*q and p and q are primes with N digits`):
print(` Try: `):
print(` ChIpre(200,10); `):

elif nops([args])=1 and op(1,[args])=ChIIpre then
print(`ChIIpre(x,K1,D1,N1,R1,k): the data for Chapter II. k challenging problems about proving that a given polynomial is non-negative`):
print(`where the answer is by displaying the polynomial as a sum of squares and there are at least 3 and at most R1 terms. Try:`):
print(`ChIIpre([x,y,z,w],10,10,10000,10,10);`):


elif nops([args])=1 and op(1,[args])=ChIIIpre then
print(`ChIIIpre(n,k,N,R1,L1,KAMA): the data for Chapter III. KAMA challenging problems about proving that a given `):
print(`sequence defined by a linear recurrence equation with polynomial coefficients with initial conditions all produce integers`):
print(`ChIIIpre(n,k,N,5,10,10);`):

elif nops([args])=1 and op(1,[args])=ChIVpre then
print(`ChIVpre(n,j,N,R1,R2,R3,MaxOrder,KAMA): the data for Chapter IV. KAMA challenging problems about proving that a given `):
print(`sequence defined by a q-linear recurrence equation with polynomial coefficients with initial conditions`):
print(`has a certain limit as n goes to infinity`):
print(`ChIVpre(n,j,N,2,3,5,3,10);`):

elif nops([args])=1 and op(1,[args])=RandBCS then
print(`RandBCS(n,k,N,R1,L1): a random binomial coefficients sum where the summand is F(n,k) . `):
print(`F(n,k) is binomial(n,k) times something of the form binomial(a1*n+b1*k+c1,k)*g^k`):
print(`where a1,b1,c1 are between 1 and R1 and g is between 1 to L1. The output is`):
print(` [F,[ope,INI]], where ope is the linear recurrence operator with polynomial coefficients (N being the shift operator)`):
print(`annihilating a(n):=Sum(F(n,k),k=0..n) and INI are the initial conditionals [a(1), ..., a(ORDER)]. Try:`):
print(`RandBCS(n,k,N,3,10);`):

elif nops([args])=1 and op(1,[args])=RandqBCS then
print(`RandqBCS(n,j,N,R1,R2,R3,MaxOrder): a random q-binomial coefficients sum inspired by George Andrews' seminal`):
print(`work on the Lusztig-Macdonald-wall conjecture. Tweaking the summand`):
print(`qbinomial(2*n,n-j)*q^(j^2) to be `):
print(`qbinomial(2*n,n-a*j)*q^(b*j*(j-1)/2+c*j) to be `):
print(`where a between 1 and R1, b is between 1 and R2 and c is between -R3 and R3 and the recurrence has at most order MaxOrder.`):
print(`It returns FAIL if the ordrer is larger.`):
print(`The output is `):
print(`[Summand,RecurrenceOperator,InitialConditions,LimitAs n goes to infinity]`):
print(`Try:`):
print(`RandqBCS(n,j,N,2,5,5,3);`):


elif nops([args])=1 and op(1,[args])=RandPrime then
print(`RandPrime(N): a random prime with N (decimal) digits. Try:`):
print(`RandPrime(200);`):

elif nops([args])=1 and op(1,[args])=RandPol then
print(`RandPol(x,K1,D1,N1): a random polynomial in the list of variables x with K1 monomials of degree <=D1, with coefficients between -N1 to N1.`):
print(`Try:`):
print(`RandPol([x,y,z],10,20,10000);`):



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


#RandPrime(N): a random prime with N (decimal) digits. Try:
#RandPrime(200);
RandPrime:=proc(N) local a,i:
a:=rand(10^N..9*10^N)():
for i from a while not isprime(i) do od:
i:
end:


#ChIpre(N,k): a list of k random integer factorization challented of the form [n,[p,q]] where n=p*q and p and q are primes with N digits
#Try:
#ChIpre(200,10);
ChIpre:=proc(N,k) local i, gu,p,q,n:
gu:=[]:

for i from 1 to k do
p:=RandPrime(N):
q:=RandPrime(N):
n:=p*q:
gu:=[ op(gu),[n,[min(p,q),max(p,q)]] ]:
od:

gu:

end:

#ChI(k,N): k challenging factorization problems
ChI:=proc(k,N) local gu,lu,i,p,q:
gu:=ChIpre(k,N):
print(``):

gu:=ChIpre(k,N):
gu:=convert(gu,list):
print(`Chapter I: `, nops(gu) , `Challenging Factorization Problems`):
print(``):



print(`Each the following positive integers is a product of two primes. Find them.`):

for i from 1 to nops(gu) do

print(` Problem `, i, `: Factorize the following integer`):
print(``):
print(gu[i][1]):
print(``):
od:
print(`--------------------------------------------`):
print(``):

print(`Answers : `):

for i from 1 to nops(gu) do
print(``):
print(`Ans. to Problem `, i, `: `):
print(``):
print(`The factorization of`, gu[i][1], ` is `):
lu:=gu[i][2]:
print(``):
print(p=lu[1]):
print(``):
print(q=lu[1]):
print(``):
od:


end:

#RandPol(x,K1,D1,N1): a random polynomial in the list of variables x with K1 monomials of degree <=D1, with coefficients between -N1 to N1.
#Try:
#RandPol[x,y,z],10,20,10000);
RandPol:=proc(x,K1,D1,N1) local ra1,ra2,gu,lu,i1,i,mu:
ra1:=rand(0..D1):
ra2:=rand(-N1..N1):

gu:=0:

for i from 1 to K1 do
 lu:=[seq(ra1(),i1=1..nops(x))]:
 mu:=ra2():
  gu:=gu+mu*mul(x[i1]^lu[i1],i1=1..nops(x)):
od:

gu:

end:


#ChIIpre(x,K1,D1,N1,R1,k): the data for Chapter II. k challenging problems about proving that a given polynomial is non-negative
#where the answer is by displaying the polynomial as a sum of squares and there are at least 3 and at most R1 terms. Try:
#ChIIpre(x,10,10,10000,10,10);
ChIIpre:=proc(x,K1,D1,N1,R1,k)
local gu,lu,P,r,i,ra,i1:
ra:=rand(3..R1):

gu:=[]:

for i from 1 to k do
r:=ra():
lu:=[seq(RandPol(x,K1,D1,N1),i1=1..r)]:
P:=expand(add(lu[i1]^2,i1=1..nops(lu))):
gu:=[op(gu),[P,lu]]:
od:

gu:

end:







#ChII(x,K1,D1,N1,R1,k): Inputs a list of variables x, positive integers K1,D1,N1,R1,k, and outputs
#a list of k challenging  problems about polynomials being non-negative. Try:
#where the answer is by displaying the polynomial as a sum of squares and there are at least 3 and at most R1 terms. Try:
#ChII([x,y,z],10,10,10000,10,10);
ChII:=proc(x,K1,D1,N1,R1,k) local gu,i,P,lu,i1:
gu:=ChIIpre(x,K1,D1,N1,R1,k):
gu:=convert(gu,list):
print(``):
print(k, `Challenging Positivity Problems for Multivariable polynomials`):
print(``):
for i from 1 to nops(gu) do
 P:=gu[i][1]:
 print(`Problem Number`, i, `Prove that the following polynomial is non-negative for all real values of the variables`, op(x)):
print(``):
print(P):
print(``):
od:
print(``):
print(`Solutions to the Problems`):
print(``):
for i from 1 to nops(gu) do
 print(``):
 print(`Proof of Problem`, i):
 print(``):
 print(`The polynomial`):
 print(``):
 print(gu[i][1]):
 print(``):
 print(`is non-negative for all its argument, since it can be expressed as a sum of squares as follows`):
 print(``):
 lu:=gu[i][2]:
 print(add(lu[i1]^2,i1=1..nops(lu))):
print(``):
 print(`Check! `):
print(``):
od:
print(``):
print(`----------------------------------`):

end:


###Start EKHAD
with(SolveTools):
solve1:=Linear:

 
 
ezraE:=proc()
if args=NULL then
print(`EKHAD`):
print(`A Maple package for proving Hypergeometric (Binomial Coeff.)`):
print(` and other kinds of identities`):
print(`This version is for  Maple 8 (thanks to Drew Sills for noticing).`):
print(`The version (Feb, 25, 1999) is much faster than the previous`):
print(`version, thanks to a SLIGHT (yet POWERFUL) modification suggested by`):
print(` FREDERIC CHYZAK `):
print(``):
print(`For help with a specific procedure, type "ezra(procedure_name);"`):
print(`Contains procedures: `):
print(`findrec,ct,zeil,zeilpap,zeillim,AZd,AZdI, AZc,AZcI, AZpapd,AZpapc,celine, TerryTao `):
fi:
 
 
if nops([args])=1 and op(1,[args])=`celine` then
print(`celine(FUNCTION,ORDER_n,ORDER_k) applies Sister Celine's technique`):
print(`It inputs a function (n,k)->FUNCTION, and the guessed orders in`):
print(`n and k, ORDER_n,ORDER_k respectively, and outputs a partial`):
print(`k-free recurrence`):
print(`it  also finds the telescoped form of the recurrence.`):   
print(` In input, do not raise products of factorials to powers;  `):  
print(`instead raise each factorial to the power. `):
print(`For example:celine((n,k)->binomial(n,k),1,1);`):
print(`celine ((n,k)->n!*(2*k)!*(-2)^(n-k)/(k!^3*(n-k)!),2,2);`):
 
fi:
 
 
if nops([args])=1 and op(1,[args])=`findrec` then
print(`findrec(f,DEG,ORDER,n,N) finds empirically an ordi. linear recurrence`):
print(` with polynomial coeffs. The input is a sequence f given as a list`):
print(`STARTING at f[1],i.e. f[0] is not considered`):
print(` where DEG:=the maximal degree of the coefficients`):
print(`and ORDER:=the order of the recurrence.  The output is the operator`):
print(` in n and N, where N is the forward unit shift: Nf(n):=f(n+1).`):
print(`For example findrec([1,1,2,3,5,8,13,21,34],0,2,n,N) should yield`):
print(`N^2-N-1 , and findrec([1,2,5,14,42,132,429],1,1,m,M) should yield`):
print(`(m+1)*M-(4*m-2). If there is not enough data, you will get an`):
print(`an error message. If there is no operator, you would get 0`):
fi:
 
if nops([args])=1 and op(1,[args])=`AZpapc` then
print(`AZpapc(INTEGRAND,y,x) inputs a hypergeometric integrand`):
print(`in the continuous variables y and x (i.e. the logarithmic derivatives`):
print(`diff(INTEGRAND,x)/INTEGRAND and diff(INTEGRAND,y)/INTEGRAND are`):
print(`rational functions in x and y)`):
print(`and outputs a paper with a proof of the linear differential equation`):
print(`that the definite integral w.r.t. to y (which is a function of x)`):
print(`satisfies. It uses the method of`):
print(`Gert Almkvist and Doron Zeilberger, "The method of differentiating`):
print(`under the integral sign", J. Symbolic Computation 10(1990), 571-591.`):
print(``):
print(`It could be used to establish the diff. eq. of the `):
print(`classical orthogonal polynomials, when they are defined in terms`):
print(`or their generating funtion.`):
print(`For example for `):
print(`the differential equation satisfied by the Legendre polynomials. Try:`):
print(`For example AZpapc(1/sqrt(1-2*x*t+t^2)/t^(n+1),t,x); `):
 
fi:
 
 
 
if nops([args])=1 and op(1,[args])=`AZpapd` then
print(`AZpapd(INTEGRAND,x,n) inputs a hypergeometric integrand`):
print(`in the continuous variable x and the discrete variable n`):
print(`i.e. i.e. A(x,n+1)/A(x,n) and A'(x,n)/A(x,n) are rational functions`):
print(`of (x,n)),`):
print(`and outputs a paper with a proof of the linear recurrence equation`):
print(`that the definite integral w.r.t. to y (which is a function of n)`):
print(`assuming that the integrand (or rather it times the certificate`):
print(`vanishes at the endpoints, or it is a  contour integrals`):
print(`satisfies. It assumes the following: A(x,n) is not a product of`):
print(`of a function of n and a function of x`):
print(`It uses the method of`):
print(`Gert Almkvist and Doron Zeilberger, "The method of differentiating`):
print(`under the integral sign", J. Symbolic Computation 10(1990), 571-591`):
print(`It could be used to establish the recurrences of the `):
print(`classical orthogonal polynomials, when they are defined in terms`):
print(`or their generating funtion`):
print(`For example  for `):
print(`the recurrence, and proof, satisfied by the Legendre polynomials, try: `):
print(`AZpapd(1/sqrt(1-2*x*t+t^2)/t^(n+1),t,n) ; `):
 
fi:
 
 
if nops([args])=1 and op(1,[args])=`AZd` then
print(`AZd(A,x,n,N) finds a recurrence of order ORDER<=8`):
print(`phrased in terms of n, and the shift in n, N`):
print(`for the integral of A(x,n) with respect to x, whenever A(x,n) is`):
print(`hypergeometric in (x,n),i.e. A(x,n+1)/A(x,n) and A'(x,n)/A(x,n) are`):
print(`rational funtions of x and n. It follows the algorithn of`):
print(`Gert Almkvist and Doron Zeilberger, "The method of differentiating`):
print(`under the integral sign", J. Symbolic Computation 10(1990), 571-591`):
print(`A should not be a product of a function of x and a function of n.`):
print(``):
print(`AZd(A,x,n,N) returns the expression seq. ope(n,N),cert(x,n)`):
print(`satisfying ope(n,N)A(x,n)=diff(cert(x,n)*A(x,n),x).`):
print(`If no recurrence is found, it returns 0.`):
print(``):
print(`A verbose version is AZpapd(A,x,n), type ezra(AZpapd) for details.`):
print(``):
print(`For example for `):
print(`the recurrence equation and proof certificate, for the Legendre polct:ynomials. Try: `):
print(`AZd(1/sqrt(1-2*x*t+t^2)/t^(n+1),t,n,N) ; `):
fi:

if nops([args])=1 and op(1,[args])=`AZdI` then
print(`AZdI(A,x,n,N,a,b) finds an INHOMOGENEOUS recurrence of order ORDER<=8`):
print(`phrased in terms of n, and the shift in n, N`):
print(`for the integral of A(x,n) with respect to x, from x=a to x=b, whenever A(x,n) is`):
print(`hypergeometric in (x,n),i.e. A(x,n+1)/A(x,n) and A'(x,n)/A(x,n) are`):
print(`rational funtions of x and n. It follows the algorithn of`):
print(`Gert Almkvist and Doron Zeilberger, "The method of differentiating`):
print(`under the integral sign", J. Symbolic Computation 10(1990), 571-591`):
print(`A should not be a product of a function of x and a function of n.`):
print(``):
print(`AZd(A,x,n,N) returns the pair [ope(n,N),rhs(n)],cert(x,n)]`):
print(`satisfying ope(n,N)A(x,n)=diff(cert(x,n)*A(x,n),x).`):
print(`and rhs(n) is the right hand side of the inhomogeneous recurrence`):
print(`If no recurrence is found, it returns 0.`):
print(``):
print(``):
print(`For example try: `):
print(`AZdI(x^n/(1+x^2),x,n,N,-1,1) ; `):
fi:
 
 
if nops([args])=1 and op(1,[args])=`AZc` then
print(`AZc(A,y,x,D) tries to finds a linear diff.eq. of order <=10,`):
print(` phrased in terms of x, and diff.w.r.t x, D`):
print(`for the integrale of A(x,y) with respect to x, whenever A(x,y) is`):
print(`hypergeometric in (x,y),i.e.  A_x(x,y)/A(x,y) and A_y(x,y)/A(x,y) are`):
print(`rational funtions of x and y. It follows the algorithn of`):
print(`Gert Almkvist and Doron Zeilberger, "The method of differentiating`):
print(`under the integral sign", J. Symbolic Computation 10(1990), 571-591`):
print(`AZc(A,y,x,D) returns the expression seq. ope(x,D),cert(x,n)`):
print(`satisfying ope(x,D)A(x,y)=diff(cert(x,y)*A(x,y),y).`):
print(`If no linear diff. eq. of order<=8 is found, it returns 0`):
print(`see AZpapc for a verbose vsersion`):
print(`For example, for `): 
print(`the diff.eq., and proof certificate for the Legendre polynomials.`):
print(` AZc(1/sqrt(1-2*x*t+t^2)/t^(n+1),t,x,D); `):
fi:

if nops([args])=1 and op(1,[args])=`AZcI` then
print(`AZcI(A,y,x,D,a,b) tries to finds an inhomogeneous linear diff.eq. of order <=10,`):
print(` phrased in terms of x, and diff.w.r.t x, D`):
print(`for the integrale of A(x,y) with respect to x from x=a to x=b, whenever A(x,y) is`):
print(`hypergeometric in (x,y),i.e.  A_x(x,y)/A(x,y) and A_y(x,y)/A(x,y) are`):
print(`rational funtions of x and y. It follows the algorithn of`):
print(`Gert Almkvist and Doron Zeilberger, "The method of differentiating`):
print(`under the integral sign", J. Symbolic Computation 10(1990), 571-591`):
print(`AZc(A,y,x,D) returns the expression seq. ope(x,D),cert(x,n)`):
print(`satisfying ope(x,D)A(x,y)=diff(cert(x,y)*A(x,y),y).`):
print(`If no linear diff. eq. of order<=10 is found, it returns 0`):
print(`For example Try:`):
print(` AZcI(1/(1-x*t-(x+1)*t^2),x,t,D,0,1); `):
print(` AZcI(x*(2-x)/(1-x*t-(x+1)*t^2),x,t,D,0,2); `):
fi:
 
 
 if nops([args])=1 and op(1,[args])=`ct` then
 
 print(` ct(SUMMAND,ORDER,k,n,N)`): 
print(`This is a Maple inplementation of the algorithm described in Ch. 6`):
print(`of the book A=B, first proposed in : D. Zeilberger, "The method of`):
print(`of creative telescoping", J.Symbolic Computation 11(1991) 195-204`):
print(``):
print(`finds a recurrence for SUMMAND in the parameters k and n,   `):
print(` of order ORDER, with N is the chosen symbol for the shift in n.`):
 print(``):
 print(`SUMMAND should be a product of factorials and/or binomial coeffs`):
 print(`and/or rising factorials, where (a)_k is denoted by rf(a,k)`):
 print(`and/or powers of k and n, and, optionally, a polynomial factor`):
print(`The output consists of an operator ope(N,n) and a certificate R(n,k)`):
print(`with the properties that if we define G(n,k):=R(n,k)*SUMMAND then`):
print(`ope(N,n)SUMMAND(n,k)=G(n,k+1)-G(n,k)`):
print(`which is a routinely verifiable identity.`):
 print(`For example "ct(binomial(n,k),1,k,n,N);" would yield the output`):
print(`  N-2, k/(k-n-1) `):
 fi:

if nops([args])=1 and op(1,[args])=TerryTao then
print(`TerryTao(P,Q,x,z,k,L,K): inputs hyper-exponential functions P and Q of the variable x, `):
print(`a continuous variable z (for the certificate),`):
print(`a discrete variable`):
print(`k, and an affine-linear expression L in k of the form a0+b0*k for some non-negative integers a0 and b0`):
print(`and a symbol,K, for the shift operator in k.`):
print(`Outputs a recurrence, and a certificate (using the Almkvist-Zeilberger algorithm) for the`):
print(`sequence , in k, of the functions (of x)`):
print(`D_x^L(Q(x)*P(x)^k). For example, for the problem in Terry Tao's blog, May 30, 2015`):
print(` A differentiation identity, type `):
print(` TerryTao(1+x^2,(1+x^2)^(-1/2),x,z,k,2*k,K); `):
fi:
 
if nops([args])=1 and op(1,[args])=`zeil` then
print(`The syntax is:`):
 print(` zeil(SUMMAND,k,n,N) or zeil(SUMMAND,k,n,N,MAXORDER) or  `):
 print(` zeil(SUMMAND,k,n,N,MAXORDER,parameter_list) `):
 print(` finds a linear recurrence equation for SUMMAND, with`):
 print(` polynomial coefficients`):
 print(`of ORDER<=MAXORDER, where the default of MAXORDER is 6`):
 print(`in the parameter n, the shift operator in n being N`):
print(`of the form ope(N,n)SUMMAND=G(n,k+1)-G(n,k)`):
print(`where G(n,k):=R(n,k)*SUMMAND, and R(n,k) is the 2nd item of output.`):
print(`The output is ope(N,n),R(n,k) .`):
 print(`For example zeil(binomial(n,k),k,n,N) would yield`):
print(`  N-2, k/(k-n-1)  `);
print(` in which N-2 is the "ope" operator, and k/(k-n-1) is R(n,k) `);
 print(`SUMMAND should be a product of factorials and/or binomial coeffs`):
 print(`and/or rising factorials, where (a)_k is denoted by rf(a,k)`):
 print(`and/or powers in k and n, and, optionally, a polynomial factor.`):
 print(``):
 print(`The last optional parameter, is the list of other parameters,`):
 print(`if present. Giving them causes considerable speedup. For example`):
 print(` zeil(binomial(n,k)*binomial(a,k)*binomial(b,k),k,n,N,6,[a,b])`):
 
fi:
 
if nops([args])=1 and op(1,[args])=`zeilpap` then
print(` zeilpap(SUMMAND,k,n) or zeilpap(SUMMAND,k,n,NAME,REF)`):
print(`Just like zeil but writes a paper with the proof`):
print(`NAME and REF are optional name and reference`):
print(`Warning: It assumes that the definite summation w.r.t. k`):
print(`extends over all k where it is non-zero, and that it is zero`):
print(`for other k`):
print(`For non-natural summation limits, use zeillim`):
fi:
 
if nops([args])=1 and op(1,[args])=`zeillim` then
 print(` zeillim(SUMMAND,k,n,N,alpha,beta) `):
print(`Similar to zeil(SUMMAND,k,n,N) but outputs a recurrence for `):
print(` the sum of SUMMAND from k=alpha to k=n-beta .`):
print(`Outputs the recurrence operator, certificate and right hand side.`):
print(`Note carefully: The upper limit of the sum will be n-beta, not beta. `):
print(`For example, "zeillim(binomial(n,k),k,n,N,0,1);" gives output of `):
print(` N-2, k/(k-n-1),1 `):
print(` which means that SUM(n):=2^n-1 satisfies the recurrence `):
print(` (N-2)SUM(n)=1, as certified by R(n,k):=k/(k-n-1) `):
fi:
 
 
end:
 
#yafe just translates from operator notation to everyday notation
 
yafe:=proc(ope,N,n,SUM)
local gu,i:
gu:=0:
 
for i from 0 to degree(ope,N) do
 gu:=gu+coeff(ope,N,i)*subs(n=n+i,SUM):
od:
 
gu:
end:
 
 
 
yafec:=proc(ope,D,x,INTEGRAND)
local gu,i:
 
 gu:=coeff(ope,D,0)*INTEGRAND:
 
 
for i from 1 to degree(ope,D) do
 gu:=gu+coeff(ope,D,i)*diff(INTEGRAND,x$i):
od:
 
gu:
end:
 
  
simplify1:=proc(bitu,k,a)
local gu,gu1,i,khez,sp:
sp:=1:
gu:=bitu:
 
if type(gu,`*`) then
 
 for i from 1 to nops(gu) do
  gu1:=op(i,gu):
 
 if type(gu1,`^`) and type(op(2,gu1), integer) then
 
     khez:=op(2,gu1):
     gu1:=op(1,gu1):
     sp:=sp*simplify(subs(k=k+a,gu1)/gu1)^khez:
  
   else 
     sp:=sp*simplify(subs(k=k+a,gu1)/gu1):
 
   fi:
 
od:
 
elif type(gu,`^`) and type(op(2,gu), integer) then
   khez:=op(2,gu):
   gu1:=op(1,gu):
   sp:=sp*simplify(subs(k=k+a,gu1)/gu1)^khez:  
 
 
else
 
 sp:=simplify(subs(k=k+a,gu)/gu):
 
fi:
 
sp:
 
end:
 
 
rf:=proc(a,b):
 (a+b-1)!/(a-1)!:
end:
 
 
 
RootOf1:=proc(f,x)
local kv,kvi,i:
kv:=[solve(f=0,x)]:
kvi:={}:
for i from 1 to nops(kv) do
 
  if type(kv[i],integer) and kv[i]>0 then
   kvi:=kvi union {kv[i]}
  fi:
 
od:
 
RETURN(kvi):
 
end:
 
pashet:=proc(p,N)
local i,gu1,gu,p1:
p1:=normal(p):
gu1:=denom(p1):
p1:=numer(p1):
p1:=expand(p1):
gu:=0:
for i from 0 to degree(p,N) do
gu:=gu+factor(coeff(p1,N,i))*N^i:
od:
RETURN(gu/gu1):
end:
 
 
 
hafel:=proc(POL,SUMMAND,ope,n,N)
local i,FAC,degg,bi,rat,POL1,SUMMAND1,zub:
degg:=degree(ope,N):
FAC:=0:
 
for i from 0 to degg do
 bi:=coeff(ope,N,i):
 rat:=simplify1(SUMMAND,n,i):
 FAC:=FAC+bi*subs(n=n+i,POL)*rat:
 FAC:=normal(FAC):
od:
 
POL1:=numer(FAC):
zub:=denom(FAC):
SUMMAND1:=SUMMAND/zub:
RETURN(POL1,SUMMAND1,zub):
end:
 
ctold:=proc(SUMMAND,ORDER,k,n,N)
local gu,i,ope,POL1,SUMMAND1,yakhas,P,Q,R,j,res1,kv,hakhi,g,
A1, B1,  P1, SDZ, SDZ1,
degg, eq, fu, gugu, i1, ia1, k1, kvutsa, l1, l2, meka1, mekb1, mumu, 
va1, va2,ope1,denFAC,ope1a:
 
if nargs<>5 then
ERROR(`Syntax: ct(SUMMAND,ORDER,summation_variable,auxiliary_var, Shift_n)`):
fi:
 
 
ope:=0:
 
for i from 0 to ORDER do
 ope:=ope+bpc[i]*N^i:
od: 
gu:=hafel(1,convert(SUMMAND,factorial),ope,n,N):
POL1:=gu[1]:
SUMMAND1:=gu[2]:
denFAC:=gu[3]:
yakhas:=simplify1(SUMMAND1,k,-1):
yakhas:=normal(1/yakhas):
 
print(`yakhas is`,yakhas):
P1:=1:
Q:=numer(yakhas):
R:=denom(yakhas):
res1:=resultant(Q,subs(k=k+j,R),k):
print(`res1 is`,res1):
kv:=RootOf1(res1,j):
 
 
 print(`kv is`,kv):
while nops(kv) > 0 do
 hakhi:=max(op(kv)):
 g:=gcd(Q,subs(k=k+hakhi,R)):
 Q:=normal(Q/g):
 R:=normal(R/subs(k=k-hakhi,g)):
 P1:=P1*product(subs(k=k-i1,g),i1=0..hakhi-1):
 res1:=resultant(Q,subs(k=k+j,R),k):
 kv:=RootOf1(res1,j):
od:
 
P:=POL1*P1:
 
A1:=subs(k=k+1,Q)+R:
A1:=expand(A1):
B1:=subs(k=k+1,Q)-R:
B1:=expand(B1):
l1:=degree(A1,k):
if B1=0 then
l2:=-100:
else
l2:=degree(B1,k):
fi:
meka1:=coeff(A1,k,l1):
mekb1:=coeff(B1,k,l2):
if l1<=l2 
then
k1:=degree(P,k)-l2:
elif l2=l1-1 
then
 mumu:= (-2)*mekb1/meka1:
  if type(mumu,integer) 
  then
   k1:=max(mumu, degree(P,k)-l1+1):
  else
   k1:=degree(P,k)-l1+1:
  fi:
elif l2 < l1-1
then
 k1:= max( 0, degree(P,k)-l1+1 ):
fi:
fu:=0:
 
 
print(`k1 is`,k1):
if k1 < 0 then
RETURN(0):
fi:
if k1 >= 0 then
for ia1 from 0 to k1 do
fu:=fu+apc[ia1]*k^ia1:
od:
gugu:=subs(k=k+1,Q)*fu-R*subs(k=k-1,fu)-P:
gugu:=expand(gugu):
degg:=degree(gugu,k):
 
for ia1 from 0 to degg do
eq[ia1+1]:=coeff(gugu,k,ia1)=0:
od:
for ia1 from 0 to k1 do
va1[ia1+1]:=apc[ia1]:
od:
for ia1 from 0 to ORDER do
va2[ia1+1]:=bpc[ia1]:
od:
eq:=convert(eq,set):
va1:=convert(va1,set):
va2:=convert(va2,set):
va1:=va1 union va2:
va1:=solve1(eq,va1):
kvutsa:={va1}:
fu:=subs(va1,fu):
ope:=subs(va1,ope):
fi:
 
 
 
if ope=0 or kvutsa={} or fu=0 then
RETURN(0):
fi:
gu:={}:
for i1 from 0 to k1 do
gu:=gu union {apc[i1]=1}:
od:
for i1 from 0 to ORDER do
gu:=gu union {bpc[i1]=1}:
od:
fu:=subs(gu,fu):
ope:=subs(gu,ope):
 
ope:=pashet(ope,N):
ope1:=ope*denom(ope):
 
SDZ:=denom(ope)*subs(k=k+1,Q)*fu/P1/denFAC :
SDZ1:=subs(k=k-1,SDZ)*simplify1(convert(SUMMAND,factorial),k,-1):
SDZ1:=simplify(SDZ1):
 
ope1a:=0:
for i from 0 to degree(ope1,N) do
ope1a:=ope1a+sort(coeff(ope1,N,i)*N^i):
od:
 
RETURN(ope1a,SDZ1):
end:
 
 
 
 
ct:=proc(SUMMAND,ORDER,k,n,N)
local gu,i,ope,POL1,SUMMAND1,yakhas,P,Q,R,j,res1,kv,hakhi,g,
A1, B1,  P1, SDZ, SDZ1,
degg, eq, fu, gugu, i1, ia1, k1, kvutsa, l1, l2, meka1, mekb1, mumu, 
va1, va2,ope1,denFAC,ope1a:
 
if nargs<>5 then
ERROR(`Syntax: ct(SUMMAND,ORDER,summation_variable,auxiliary_var, Shift_n)`):
fi:
 
if tovu(convert(SUMMAND,factorial),k,n)=0 then
ERROR(`The summand can be separated into f(`,k,`)g(`,n,`)`): 
fi:
 
ope:=0:
 
for i from 0 to ORDER do
 ope:=ope+bpc[i]*N^i:
od:
 
gu:=hafel(1,convert(SUMMAND,factorial),ope,n,N):
POL1:=gu[1]:
SUMMAND1:=gu[2]:
denFAC:=gu[3]:
yakhas:=simplify1(SUMMAND1,k,-1):
yakhas:=normal(1/yakhas):
 
P1:=1:
Q:=numer(yakhas):
R:=denom(yakhas):
 
 
res1:=findgQR(Q,R,k,100):
 
while res1[2]<>0 do
 g:=res1[1]:
 hakhi:=res1[2]:
 Q:=normal(Q/g):
 R:=normal(R/subs(k=k-hakhi,g)):
 P1:=P1*product(subs(k=k-i1,g),i1=0..hakhi-1):
 res1:=findgQR(Q,R,k,100):
od:
 
P:=POL1*P1:
 
A1:=subs(k=k+1,Q)+R:
A1:=expand(A1):
B1:=subs(k=k+1,Q)-R:
B1:=expand(B1):
l1:=degree(A1,k):
if B1=0 then
l2:=-100:
else
l2:=degree(B1,k):
fi:
meka1:=coeff(A1,k,l1):
mekb1:=coeff(B1,k,l2):
if l1<=l2 
then
k1:=degree(P,k)-l2:
elif l2=l1-1 
then
 mumu:= (-2)*mekb1/meka1:
  if type(mumu,integer) 
  then
   k1:=max(mumu, degree(P,k)-l1+1):
  else
   k1:=degree(P,k)-l1+1:
  fi:
elif l2 < l1-1
then
 k1:= max( 0, degree(P,k)-l1+1 ):
fi:
fu:=0:
 
 
if k1 < 0 or k1=FAIL then
RETURN(0):
fi:
if k1 >= 0 then
for ia1 from 0 to k1 do
fu:=fu+apc[ia1]*k^ia1:
od:
gugu:=subs(k=k+1,Q)*fu-R*subs(k=k-1,fu)-P:
gugu:=expand(gugu):
degg:=degree(gugu,k):
 
for ia1 from 0 to degg do
eq[ia1+1]:=coeff(gugu,k,ia1)=0:
od:
for ia1 from 0 to k1 do
va1[ia1+1]:=apc[ia1]:
od:
for ia1 from 0 to ORDER do
va2[ia1+1]:=bpc[ia1]:
od:
eq:=convert(eq,set):
va1:=convert(va1,set):
va2:=convert(va2,set):
va1:=va1 union va2:
va1:=solve1(eq,va1):
kvutsa:={va1}:
fu:=subs(va1,fu):
ope:=subs(va1,ope):
fi:
 
 
 
if ope=0 or kvutsa={} or fu=0 then
RETURN(0):
fi:
gu:={}:
for i1 from 0 to k1 do
gu:=gu union {apc[i1]=1}:
od:
for i1 from 0 to ORDER do
gu:=gu union {bpc[i1]=1}:
od:
fu:=subs(gu,fu):
ope:=subs(gu,ope):
 
ope:=pashet(ope,N):
ope1:=ope*denom(ope):
 
SDZ:=denom(ope)*subs(k=k+1,Q)*fu/P1/denFAC :
SDZ1:=subs(k=k-1,SDZ)*simplify1(convert(SUMMAND,factorial),k,-1):
SDZ1:=simplify(SDZ1):
 
ope1a:=0:
for i from 0 to degree(ope1,N) do
ope1a:=ope1a+sort(coeff(ope1,N,i)*N^i):
od:
 
RETURN(ope1a,SDZ1):
end:
 
 
 
cttry:=proc(SUMMAND,ORDER,k,n,resh,N)
local gu,i,ope,POL1,SUMMAND1,yakhas,P,Q,R,j,res1,kv,hakhi,g,
A1, B1,  P1, 
degg, eq, fu, gugu, i1, ia1, k1, kvutsa, l1, l2, meka1, mekb1, mumu, 
va1, va2,eqg:
 
if nargs<>6 then
ERROR(`The syntax is cttry(term,ORDER,sum_variable,aux_var,para_list,Shift)`):
fi:
 
if tovu(convert(SUMMAND,factorial),k,n)=0 then
ERROR(`The summand can be separated into f(`,k,`)g(`,n,`)`): 
fi:
 
ope:=0:
 
for i from 0 to ORDER do
 ope:=ope+bpc[i]*N^i:
od:
 
gu:=hafel(1,convert(SUMMAND,factorial),ope,n,N):
POL1:=gu[1]:
SUMMAND1:=gu[2]:
yakhas:=simplify1(SUMMAND1,k,-1):
yakhas:=normal(1/yakhas):
 
P1:=1:
Q:=numer(yakhas):
R:=denom(yakhas):
 
res1:=findgQR(Q,R,k,100):
 
while res1[2]<>0 do
 g:=res1[1]:
 hakhi:=res1[2]:
 Q:=normal(Q/g):
 R:=normal(R/subs(k=k-hakhi,g)):
 P1:=P1*product(subs(k=k-i1,g),i1=0..hakhi-1):
 res1:=findgQR(Q,R,k,100):
od:
 
P:=POL1*P1:
 
A1:=subs(k=k+1,Q)+R:
A1:=expand(A1):
B1:=subs(k=k+1,Q)-R:
B1:=expand(B1):
l1:=degree(A1,k):
if B1=0 then
l2:=-100:
else
l2:=degree(B1,k):
fi:
meka1:=coeff(A1,k,l1):
mekb1:=coeff(B1,k,l2):
if l1<=l2 
then
k1:=degree(P,k)-l2:
if k1=FAIL then
 RETURN(FAIL):
fi:
elif l2=l1-1 
then
 mumu:= (-2)*mekb1/meka1:
  if type(mumu,integer) 
  then
   k1:=max(mumu, degree(P,k)-l1+1):
  else
   k1:=degree(P,k)-l1+1:
  fi:
elif l2 < l1-1
then
 k1:= max( 0, degree(P,k)-l1+1 ):
fi:
fu:=0:
 
 
if k1 < 0 then
RETURN(0):
fi:
if k1 >= 0 then
for ia1 from 0 to k1 do
fu:=fu+apc[ia1]*k^ia1:
od:
gugu:=subs(k=k+1,Q)*fu-R*subs(k=k-1,fu)-P:
gugu:=expand(gugu):
degg:=degree(gugu,k):
for ia1 from 0 to degg do
eq[ia1+1]:=coeff(gugu,k,ia1)=0:
od:
for ia1 from 0 to k1 do
va1[ia1+1]:=apc[ia1]:
od:
for ia1 from 0 to ORDER do
va2[ia1+1]:=bpc[ia1]:
od:
eq:=convert(eq,set):
va1:=convert(va1,set):
va2:=convert(va2,set):
va1:=va1 union va2:
eqg:=subs(n=1/2,eq):
 
for i from 1 to nops(resh) do
 eqg:=subs(op(i,resh)=i^2+3,eqg):
od:
 
va1:=solve1(eq,va1):
kvutsa:={va1}:
fu:=subs(va1,fu):
ope:=subs(va1,ope):
fi:
 
 
 
if ope=0 or kvutsa={} or fu=0 then
RETURN(0):
 else
  RETURN(1):
fi:
 
end:
 
 
paper:=proc(SUMMAND,k,n,N,ope1,SDZ1,NAME,REF)
local SHEM,IDENTITY,RECURRENCE:
 
if degree(ope1,N)=1 then
SHEM:=IDENTITY:
else
SHEM:=RECURRENCE:
fi:
print(``):
print(`A PROOF OF THE`,NAME,SHEM):
print(``):
print(`By Shalosh B. Ekhad,Rutgers University, c/o zeilberg@math.rutgers.edu`):
print(``):
print(`I will give a short proof of the following result(`,REF,`).`):
print(``):
if degree(ope1,N)=1 then
print(`(Note that since the recurrence below is first order, this`):
print(`means that the sum`, SUM(n), `has closed form,and it is`):
print(`easily seen to be equivalent.)`):
print(``):
fi:
print(`Theorem:Let`, F(n,k), `be given by`):
print(``):
print(SUMMAND):
print(``):
print(`and let`, SUM(n),`be the sum of`, F(n,k),` with`):
print(`respect to`, k):
print(``):
print(SUM(n),` satisfies the following linear recurrence equation`):
print(``):
print(yafe(ope1,N,n,SUM(n))):
print(`=0.`):
print(``):
print(`PROOF: We cleverly construct`, G(n,k),`:=`):
print(``):
print(SDZ1*SUMMAND):
print(`with the motive that`):
print(``):
print(yafe(ope1,N,n,F(n,k)) = G(n,k+1)-G(n,k), `(Check!)`):
print(``):
print(`and the theorem follows upon summing with respect to`, k,`. QED.`):
end:
 
 
 
paper3:=proc(SUMMAND,k,n,N,ope1,SDZ1):
 
print(``):
print(`A PROOF OF A RECURRENCE`):
print(``):
print(`By Shalosh B. Ekhad, Rutgers University `):
print(``):
 
print(`Theorem: Let`, F(n,k), `be given by`):
print(``):
print(SUMMAND):
print(``):
print(`and let`, SUM(n),`be the sum of`, F(n,k),` with`):
print(`respect to`, k):
print(``):
print(SUM(n),` satisfies the following linear recurrence equation`):
print(``):
print(yafe(ope1,N,n,SUM(n))):
print(`=0.`):
print(``):
print(`PROOF: We cleverly construct`, G(n,k),`:=`):
print(``):
print(SDZ1*SUMMAND):
print(`with the motive that`):
print(``):
print(yafe(ope1,N,n,F(n,k))=G(n,k+1) - G(n,k), `(Check!)`):
print(``):
print(`and the theorem follows upon summing with respect to`, k,`. QED.`):
end:
 
 
paperc:=proc(INTEGRAND,y,x,D,ope1,SDZ1):
print(``):
print(``):
print(`A PROOF OF A DIFFERENTIAL EQUATION SATISFIED BY AN INTEGRAL`):
print(``):
print(`By Shalosh B. Ekhad, Rutgers University, c/o zeilberg@math.rutgers.edu`):
print(``):
print(`I will give a short proof of the following result.`):
print(``):
 
print(`Theorem:Let`, F(x,y), `be given by`):
print(``):
print(INTEGRAND):
print(``):
print(`and let`, INTEGRAL(x),`be the integral of`, F(x,y),` with`):
print(`respect to`, y):
print(``):
print(INTEGRAL(x),` satisfies the following linear differential equation`):
print(``):
print(yafec(ope1,D,x,INTEGRAL(x))):
print(`=0.`):
print(``):
print(`PROOF: We cleverly construct`, G(x,y),`:=`):
print(``):
print(SDZ1*INTEGRAND):
print(`with the motive that`):
print(``):
print(yafec(ope1,D,x,F(x,y)) = diff(G(x,y),y),  `(Check!)`):
print(``):
print(`and the theorem follows upon integrating with respect to`, y):
end:
 
paperd:=proc(INTEGRAND,x,n,N,ope1,SDZ1):
print(``):
print(`A PROOF OF A LINEAR RECURRENCE SATISFIED BY AN INTEGRAL`):
print(``):
print(`By Shalosh B. Ekhad, Rutgers University, c/o zeilberg@math.rutgers.edu`):
print(``):
print(`I will give a short proof of the following result.`):
print(``):
 
print(`Theorem:Let`, F(n,x), `be given by`):
print(``):
print(INTEGRAND):
print(``):
print(`and let`, INTEGRAL(n),`be the integral of`, F(n,x),` with`):
print(`respect to`, x):
print(``):
print(INTEGRAL(n),` satisfies the following linear recurrence equation`):
print(``):
print(yafe(ope1,N,n,INTEGRAL(n))):
print(`=0.`):
print(``):
print(`PROOF: We cleverly construct`, G(n,x),`:=`):
print(``):
print(SDZ1*INTEGRAND):
print(`with the motive that`):
print(``):
print(yafe(ope1,N,n,F(n,x)) = diff(G(n,x),x)):
print(``):
print(`and the theorem follows upon integrating with respect to`, x,`. QED.`):
end:
 
 
paperc:=proc(INTEGRAND,y,x,D,ope1,SDZ1):
print(`A PROOF OF A DIFFERENTIAL EQUATION SATISFIED BY AN INTEGRAL`):
print(``):
print(`By Shalosh B. Ekhad, Rutgers University`):
print(``):
print(`I will give a short proof of the following result.`):
print(``):
 
print(`Theorem:Let`, F(x,y), `be given by`):
print(``):
print(INTEGRAND):
print(``):
print(`and let`, INTEGRAL(x),`be the integral of`, F(x,y),` with`):
print(`respect to`, y):
print(``):
print(INTEGRAL(x),` satisfies the following linear differential equation`):
print(``):
print(yafec(ope1,D,x,INTEGRAL(x))):
print(`=0.`):
print(``):
print(`PROOF: We cleverly construct`, G(x,y),`:=`):
print(``):
print(SDZ1*INTEGRAND):
print(`with the motive that`):
print(``):
print(yafec(ope1,D,x,F(x,y)) = diff(G(x,y),y), `(Check!)`):
print(``):
print(`and the theorem follows upon integrating with respect to`, y,`. QED.`):
end:
 
 
 
zeil4:=proc(SUMMAND,k,n,N)
local ORDER,MAXORDER,gu,gu1,resh:
 
 
 MAXORDER:=5:
 
  resh:=[]:
 
 
if simplify1(SUMMAND,n,1)=1 then
ERROR(`Summand does not depend on`,n,` hence the sum is identically constant`):
fi:
for ORDER from 0 to MAXORDER do
  gu:=ct(SUMMAND,ORDER,k,n,N):
 
   if gu<>0 and gu<>FAIL then
     RETURN(gu):
  fi:

od:
 

FAIL:
 
end:
 
 
 
 
zeil5:=proc(SUMMAND,k,n,N,MAXORDER)
local ORDER,gu,gu1,resh:
 
resh:=[]:
 
if simplify1(SUMMAND,n,1)=1 then
ERROR(`Summand does not depend on`,n,` hence the sum is identically constant`):
fi:
for ORDER from 0 to MAXORDER do
  gu:=ct(SUMMAND,ORDER,k,n,N):
 
   if gu<>0 and gu<>FAIL then
     RETURN(gu):
   fi:

od:
 
print(`No recurrence of order<=`,MAXORDER,`was found`):
RETURN(FAIL):
 
end:
 
 
 
 
zeil6:=proc(SUMMAND,k,n,N,MAXORDER,resh)
local ORDER,gu,gu1:
 
 
if simplify1(SUMMAND,n,1)=1 then
ERROR(`Summand does not depend on`,n,` hence the sum is identically constant`):
fi:
for ORDER from 0 to MAXORDER do
  gu:=ct(SUMMAND,ORDER,k,n,N):
 
   if gu<>0 then
     RETURN(gu):
   fi:

od:
 
print(`No recurrence of order<=`,MAXORDER,`was found`):
 
end:
 
 
zeil:=proc():
 
if nargs=4 then
 zeil4(args):
elif nargs=5 then
 zeil5(args):
elif nargs=6 then
 zeil6(args):
else
  print(`zeil(SUMMAND,k,n,N) or zeil(SUMMAND,k,n,N,MAXORDER) or`):
  ERROR(`zeil(SUMMAND,k,n,N,MAXORDER,param_list)`):
fi:
 
 
 
end:
 
 
 
 
 
zeillim:=proc(SUMMAND,k,n,N,alpha,beta)
local ope,CERT,lu,k1,i,gu,lu1,lu2:
 
gu:=Zeillim(SUMMAND,k,n,N,alpha+1,beta+1):
ope:=gu[1]:
CERT:=gu[2]:
lu:=gu[3]: 
lu1:=subs(k=alpha,SUMMAND)+subs(k=n-beta,SUMMAND):
lu1:=simplify(lu1):
lu1:=normal(lu1):
lu2:=0:
for i from 0 to degree(ope,N) do
   lu2:=lu2+coeff(ope,N,i)*simplify(subs(n=n+i,lu1)):
 od:
 lu2:=normal(lu2):
ope,CERT,normal(expand(normal(simplify(expand(normal(lu+lu2)))))):
 
end:
 
Zeillim:=proc(SUMMAND,k,n,N,alpha,beta)
local ope,CERT,lu,k1,i,gu:
 
gu:=zeil(SUMMAND,k,n,N):
ope:=gu[1]:
CERT:=gu[2]:
 
lu:=simplify(subs(k=n-beta+1,CERT)*subs(k=n-beta+1,SUMMAND))
-simplify(subs(k=alpha,CERT)*subs(k=alpha,SUMMAND)):
 
for i from 0 to degree(ope,N) do
 for k1 from 1 to i do
   lu:=lu+coeff(ope,N,i)*simplify(subs({n=n+i,k=n-beta+k1},SUMMAND)):
 od:
od:
 
gu,normal(expand(lu)):
 
end:
 
 
 
zeilpap3:=proc(SUMMAND,k,n)
local  SDZ1, gu, ope1,N:
 
gu:=zeil4(SUMMAND,k,n,N):
ope1:=gu[1]:
SDZ1:=gu[2]:
paper3(SUMMAND,k,n,N,ope1,SDZ1):
end:
 
 
 
 
zeilpap5:=proc(SUMMAND,k,n,NAME,REF)
local SDZ1, gu, ope1,N:
 
gu:=zeil4(SUMMAND,k,n,N):
ope1:=gu[1]:
SDZ1:=gu[2]:
paper(SUMMAND,k,n,N,ope1,SDZ1,NAME,REF):
end:
 
 zeilpap:=proc()
 
if nargs=5 then
zeilpap5(args):
elif nargs=3 then
zeilpap3(args):
else
 ERROR(`zeilpap(SUMMAND,k,n,NAME,REF) or  zeilpap(SUMMAND,k,n)`):
fi:
end:
 
 
 
 
AZpapc:=proc(INTEGRAND,y,x)
local D,SDZ1,gu,ope1:
gu:=AZc(INTEGRAND,y,x,D):
ope1:=gu[1]:
SDZ1:=gu[2]:
paperc(INTEGRAND,y,x,D,ope1,SDZ1):
end:
 
 
 
 
AZpapd:=proc(INTEGRAND,x,n)
local D,SDZ1, gu, ope1:
gu:=AZd(INTEGRAND,x,n,D):
ope1:=gu[1]:
SDZ1:=gu[2]:
paperd(INTEGRAND,x,n,D,ope1,SDZ1):
end:
 
 
 
 
 
goremd:=proc(N,ORDER)
local i, gu:
gu:=0:
for i from 0 to ORDER do
gu:=gu+(bpc[i])*N^i:
od:
RETURN(gu):
end:
 
goremc:=proc(D,ORDER)
local i,gu:
gu:=0:
for i from 0 to ORDER do
gu:=gu+(bpc[i])*D^i:
od:
RETURN(gu):
end:
 
gzor:=proc(f,x,n)
local i,gu:
gu:=f:
for i from 1 to n do
gu:=diff(gu,x):
od:
RETURN(gu):
end:
 
gzor1:=proc(a,D,x)
local i,gu:
gu:=0:
for i from 0 to degree(a,D) do
gu:=gu+diff(coeff(a,D,i),x)*D^i+coeff(a,D,i)*D^(i+1):
od:
end:
 
pashetd:=proc(p,N)
local gu,p1,mekh,i:
p1:=normal(p):
mekh:=denom(p1):
p1:=numer(p1):
p1:=expand(p1):
gu:=0:
for i from 0 to degree(p1,N) do
gu:=gu+factor(coeff(p1,N,i))*N^i:
od:
RETURN(gu,mekh):
end:
 
pashetc:=proc(p,D)
local gu,p1,i,mekh:
p1:=normal(p):
mekh:=denom(p1):
p1:=numer(p1):
p1:=expand(p1):
gu:=0:
for i from 0 to degree(p1,D) do
gu:=gu+factor(coeff(p1,D,i))*D^i:
od:
RETURN(gu,mekh):
end:
 
 
AZd:= proc(A,y,n,N)
local ORDER,gu,KAMA:
KAMA:=8:
 
for ORDER from 0 to KAMA do
gu:=duisd(A,ORDER,y,n,N):
if gu<>0 then
 RETURN(gu):
fi:
od:
0:
end:
 
 
AZc:=proc(A,y,x,D)
local ORDER,gu,KAMA:
KAMA:=10:
 
for ORDER from 1 to KAMA do
gu:=duisc(A,ORDER,y,x,D):
if gu<>0 and gu[1]<>0 then
 RETURN(gu):
fi:
od:
0:
end:
 
 
duisd:= proc(A,ORDER,y,n,N)
local gorem, conj, yakhas,lu1a,LU1,P,Q,R,S,j1,g,Q1,Q2,l1,l2,mumu,
mekb1,fu,meka1,k1,gugu,eq,ia1,va1,va2,degg,shad,KAMA,i1:
KAMA:=10:
gorem:=goremd(N,ORDER):
 
 
conj:=gorem:
yakhas:=0:
 
for i1 from 0 to degree(conj,N) do
 
yakhas:=yakhas+coeff(conj,N,i1)*simplify(subs(n=n+i1,A)/A):
od:
 
lu1a:=yakhas:
LU1:=numer(yakhas):
yakhas:=1/denom(yakhas):
yakhas:=normal(diff(yakhas,y)/yakhas+diff(A,y)/A):
P:=LU1:
Q:=numer(yakhas):
R:=denom(yakhas):
j1:=0:
while j1 <= KAMA do
g:=gcd(R,Q-j1*diff(R,y)):
if g <> 1 then
Q2:=(Q-j1*diff(R,y))/g:
R:=normal(R/g):
Q:=normal(Q2+j1*diff(R,y)):
P:=P*g^j1:
j1:=-1:
fi:
j1:=j1+1:
od:
P:=expand(P):
R:=expand(R):
Q:=expand(Q):
Q1:=Q+diff(R,y):
Q1:=expand(Q1):
l1:=degree(R,y):
l2:=degree(Q1,y):
meka1:=coeff(R,y,l1):
mekb1:=coeff(Q1,y,l2):
if l1-1 <>l2 
then
k1:=degree(P,y)-max(l1-1,l2):
else
 mumu:= -mekb1/meka1:
  if type(mumu,integer) and mumu > 0 
  then
   k1:=max(mumu, degree(P,y)-l2):
  else
   k1:=degree(P,y)-l2:
  fi:
fi:
fu:=0:
if k1 < 0 then
RETURN(0):
fi:
for ia1 from 0 to k1 do
fu:=fu+apc[ia1]*y^ia1:
od:
gugu:=Q1*fu+R*diff(fu,y)-P:
gugu:=expand(gugu):
degg:=degree(gugu,y):
for ia1 from 0 to degg do
eq[ia1+1]:=coeff(gugu,y,ia1)=0:
od:
for ia1 from 0 to k1 do
va1[ia1+1]:=apc[ia1]:
od:
for ia1 from 0 to ORDER do
va2[ia1+1]:=bpc[ia1]:
od:
eq:=convert(eq,set):
va1:=convert(va1,set):
va2:=convert(va2,set):
va1:=va1 union va2 :
va1:=solve1(eq,va1):
fu:=subs(va1,fu):
gorem:=subs(va1,gorem):
if fu=0 and gorem=0 then
 RETURN(0):
fi:
for ia1 from 0 to k1 do
gorem:=subs(apc[ia1]=1,gorem):
od:
for ia1 from 0 to ORDER do
gorem:=subs(bpc[ia1]=1,gorem):
od:
fu:=normal(fu):
 
 
shad:=pashetd(gorem,N):
S:=lu1a*R*fu*shad[2]/P:
S:=subs(va1,S):
S:=normal(S):
S:=factor(S):
for ia1 from 0 to k1 do
S:=subs(apc[ia1]=1,S):
od:
for ia1 from 0 to ORDER do
S:=subs(bpc[ia1]=1,S):
od:
 
RETURN(shad[1],S):
end:
 
 
duisc:= proc(A,ORDER,y,x,D)
local KAMA,gorem,conj, yakhas,lu1a,LU1,P,Q,R,S,j1,g,Q1,Q2,l1,l2,mumu,
mekb1,fu,meka1,k1,gugu,eq,ia1,va1,va2,degg,i,shad:
 
KAMA:=10:
gorem:=goremc(D,ORDER):
 
conj:=gorem:
yakhas:=0:
 
for i from 0 to degree(conj,D) do
yakhas:=yakhas+normal(coeff(conj,D,i)*gzor(A,x,i)/A):
yakhas:=normal(yakhas):
od:
 
lu1a:=yakhas:
LU1:=numer(yakhas):
yakhas:=1/denom(yakhas):
yakhas:=normal(diff(yakhas,y)/yakhas+diff(A,y)/A):
P:=LU1:
Q:=numer(yakhas):
R:=denom(yakhas):
j1:=0:
while j1 <= KAMA do
g:=gcd(R,Q-j1*diff(R,y)):
if g <> 1 then
Q2:=(Q-j1*diff(R,y))/g:
R:=normal(R/g):
Q:=normal(Q2+j1*diff(R,y)):
P:=P*g^j1:
j1:=-1:
fi:
j1:=j1+1:
od:
P:=expand(P):
R:=expand(R):
Q:=expand(Q):
Q1:=Q+diff(R,y):
Q1:=expand(Q1):
l1:=degree(R,y):
l2:=degree(Q1,y):
meka1:=coeff(R,y,l1):
mekb1:=coeff(Q1,y,l2):
if l1-1 <>l2 
then
k1:=degree(P,y)-max(l1-1,l2):
else
 mumu:= -mekb1/meka1:
  if type(mumu,integer) and mumu > 0 
  then
   k1:=max(mumu, degree(P,y)-l2):
  else
   k1:=degree(P,y)-l2:
  fi:
fi:
fu:=0:
if k1 < 0 then
RETURN(0):
fi:
for ia1 from 0 to k1 do
fu:=fu+apc[ia1]*y^ia1:
od:
gugu:=Q1*fu+R*diff(fu,y)-P:
gugu:=expand(gugu):
degg:=degree(gugu,y):
for ia1 from 0 to degg do
eq[ia1+1]:=coeff(gugu,y,ia1)=0:
od:
for ia1 from 0 to k1 do
va1[ia1+1]:=apc[ia1]:
od:
for ia1 from 0 to ORDER do
va2[ia1+1]:=bpc[ia1]:
od:
eq:=convert(eq,set):
va1:=convert(va1,set):
va2:=convert(va2,set):
va1:=va1 union va2 :
va1:=solve1(eq,va1):
fu:=subs(va1,fu):
gorem:=subs(va1,gorem):
if fu=0 and gorem=0 then
 RETURN(0):
fi:
for ia1 from 0 to k1 do
gorem:=subs(apc[ia1]=1,gorem):
od:
for ia1 from 0 to ORDER do
gorem:=subs(bpc[ia1]=1,gorem):
od:
fu:=normal(fu):
shad:=pashetc(gorem,D):
S:=lu1a*R*fu*shad[2]/P:
S:=subs(va1,S):
S:=normal(S):
S:=factor(S):
for ia1 from 0 to k1 do
S:=subs(apc[ia1]=1,S):
od:
for ia1 from 0 to ORDER do
S:=subs(bpc[ia1]=1,S):
od:
shad[1],S:
end:
 
 
 
bdokcertc:=proc(A,y,x,D,ope,S)
local gu,i:
gu:=0:
 
for i from 0 to degree(ope,D) do
gu:=gu+coeff(ope,D,i)*simplify(gzor(A,x,i)/A):
gu:=normal(gu):
od:
 
gu:=gu/simplify(diff(S*A,y)/A):
normal(gu);
end:
 
 
bdokcertd:=proc(A,y,n,N,ope,S)
local gu,i:
gu:=0:
 
for i from 0 to degree(ope,N) do
gu:=gu+coeff(ope,N,i)*simplify( subs(n=n+i,A)/A):
gu:=normal(gu):
od:
 
gu:=gu/simplify(diff(S*A,y)/A):
normal(gu);
end:
 
 
 
bdokcto:=proc(SUMMAND1,ORDER,k,n,N)
local mu,gu,i,G1,ope,lu,SUMMAND:
 
SUMMAND:=convert(SUMMAND1,factorial):
mu:=ct(SUMMAND,ORDER,k,n,N):
 
if mu=0 then
 RETURN(0):
fi:
 
ope:=mu[1]:
G1:=mu[2]*SUMMAND:
gu:=0:
 
for i from 0 to degree(ope,N) do
gu:=gu+coeff(ope,N,i)*simplify( subs(n=n+i,SUMMAND)/SUMMAND):
gu:=normal(gu):
od:
 
lu:=simplify(subs(k=k+1,G1)/SUMMAND)-mu[2]:
lu:=normal(lu):
normal(gu/lu);
end:
 
 
 
 
bdokct:=proc(SUMMAND1,ORDER,k,n,N)
local mu,gu,i,G1,ope,lu,SUMMAND:
 
SUMMAND:=convert(SUMMAND1,factorial):
mu:=ct(SUMMAND,ORDER,k,n,N):
 
if mu=0 then
 RETURN(0):
fi:
 
ope:=mu[1]:
G1:=mu[2]*SUMMAND:
gu:=0:
 
for i from 0 to degree(ope,N) do
gu:=gu+coeff(ope,N,i)*simplify1(SUMMAND,n,i):
gu:=normal(gu):
od:
 
lu:=subs(k=k+1,mu[2])*simplify1(SUMMAND,k,1)-mu[2]:
lu:=normal(lu):
normal(gu/lu);
end:
 
 
findrec:=proc(f,DEGREE,ORDER,n,N)
local ope,var,eq,a,i,j,n0,kv,var1,eq1,mu:
if (1+DEGREE)*(1+ORDER)+2+ORDER>nops(f) then
ERROR(`Insufficient data for a recurrence of order`,ORDER, `degree`,DEGREE):
fi:
ope:=0:
var:={}:
 
for i from 0 to ORDER do
 for j from 0 to DEGREE do
  ope:=ope+a[i,j]*n^j*N^i:
  var:=var union {a[i,j]}:
 od:
od:
    
eq:={}:
 
for n0 from 1 to (1+DEGREE)*(1+ORDER)+2 do
  eq1:=0:
 
  for i from 0 to ORDER do
     eq1:=eq1+subs(n=n0,coeff(ope,N,i))*op(n0+i,f):
  od:
 
   eq:= eq union {eq1}:
od:
 
var1:=solve1(eq,var):
 
kv:={}:
 
for i from 1 to nops(var1) do
 mu:=op(i,var1):
 
if op(1,mu)=op(2,mu) then
   kv:= kv union {op(1,mu)}:
 fi:
 
od:
ope:=subs(var1,ope):
 
 
if nops(kv)>1 then
  print(` either DEGREE or ORDER are too high`):
   print(`The output is not the minimal possible operator`):
fi:
 
for i from 1 to nops(kv) do
  ope:=subs(op(i,kv)=1,ope):
od:
 
ope:
 
end:
        
celine:=proc(f,ii,jj)
local m,j,l,i,tt,yy,zz,rr,xx,ss,zy,
 b, iii, k, n, r, zzz:
 
m:='m':j:='j':l:='l':i:='i':
iii:=ii*(jj+1)+jj:
xx:=seq(b[i],i=0..iii):i:='i':
ss:=numer(simplify(expand(sum(sum(b[i*(jj+1)+j]*f(n-j,k-i)/f(n,k), 
            i=0..ii),j=0..jj)))):
tt:=coeffs(collect(ss,k),k):
yy:=solve1({tt},{xx}):
zz:=simplify(sum(sum(b[l*(jj+1)+m]*F(n-m,k-l),l=0..ii),m=0..jj)): 
i:='i':j:='j':print(` ` ):print(` `):print(`The full  recurrence is`):
print(` `):
zz:=numer(simplify(subs(yy,zz))):
for i from 0 to ii do for j from 0 to jj do                
zz:=collect(zz,F(n-j,k-i)) od od:
print(zz,`==0`);zzz:=zz:
l:='l':j:='j':m:='m':
r:=subs(yy,simplify(expand(sum(sum(sum(b[l*(jj+1)+m],l=j+1..ii)*simplify(expand(f(n-m,k-j)/f(n,k))),j=0..ii),m=0..jj)))): 
 
rr:=simplify(factor(simplify(r))):print(` ` ):print(` `):
print(`The telescoped form is`);print(` `);
zy:=factor(subs(yy,sum(simplify(sum(b[l*(jj+1)+m],l=0..ii))*F(n-m,k),m=0..jj))):
print(zy,`==G(n,k)-G(n,k-1)`); 
 
print(` `);
print(`where G(n,k)=R(n,k)*F(n,k) and the rational function  R(n,k) is `):
print(` `);
rr;
end:
 
findgQR:=proc(Q,R,k,L)
local j,g:
 
for j from 1 to L do
 g:=gcd(Q,subs(k=k+j,R)):
 if degree(g,k)>0 then
  RETURN(g,j):
 fi:
od:
1,0:
end:
 
  
 
 
tovu:=proc(SU,k,n)
local gu:
 
gu:=simplify1(simplify1(SU,n,1),k,1):
 
if gu=1 then
RETURN(0):
else
RETURN(1):
fi:
end:
 
#TerryTao(P,Q,x,z,k,L,K): inputs hyper-exponential functions P and Q of the variable x, 
#a continuous variable z (for the certificate),
#a discrete variable
#k, and an affine-linear expression L in k of the form a0+b0*k for some non-negative integers a0 and b0
#and a symbol,K, for the shift operator in k.
#Outputs a recurrence, and a certificate (using the Almkvist-Zeilberger algorithm) for the
#sequence , in k, of the functions (of x)
#D_x^L(Q(x)*P(x)^k). For example, for the problem in Terry Tao's blog, May 30, 2015
#A differentiation identity, type
TerryTao(1+x^2,(1+x^2)^(-1/2),x,z,k,2*k,K);

TerryTao:=proc(P,Q,x,z,k,L,K) :

AZd(L!*subs(x=z,P)^k*subs(x=z,Q)/(z-x)^(L+1),z,k,K):


end:




####start Dec. 11, 2015 paper
AZdI:= proc(A,y,n,N,a,b)
local ORDER,gu,KAMA,ope,cert,yemin1:
KAMA:=12:
 
for ORDER from 0 to KAMA do
gu:=duisd(A,ORDER,y,n,N):
if gu<>0 then

ope:=gu[1]:
cert:=gu[2]:
yemin1:=normal(subs(y=b, normal(cert*A))-subs(y=a, normal(cert*A))):
RETURN([[ope,yemin1],cert]):

fi:
od:
0:
end:



AZcI:=proc(A,y,x,D,a,b) local ORDER,gu,KAMA,ope,cert,yemin1:
KAMA:=8:
 
for ORDER from 1 to KAMA do
gu:=duisc(A,ORDER,y,x,D):
if gu<>0 and gu[1]<>0 then
ope:=gu[1]:
cert:=gu[2]:
yemin1:=normal(subs(y=b, normal(cert*A))-subs(y=a, normal(cert*A))):
RETURN([[ope,yemin1],cert]):

fi:
od:
0:
end:
###End EKHAD


#RandBCS(n,k,N,R1,L1): a random binomial coefficients sum where the summand is F(n,k) . 
#F(n,k) is binomial(n,k) times something of the form binomial(a1*n+b1*k+c1,k)*g^k
#where a1,b1,c1 are between 1 and R1 and g is between 1 to L1. The output is
# [F,[ope,INI]], where ope is the linear recurrence operator with polynomial coefficients (N being the shift operator)
#annihilating a(n):=Sum(F(n,k),k=0..n) and INI are the initial conditionals [a(1), ..., a(ORDER)]. 
#RandBCS(n,k,2,3,5);

RandBCS:=proc(n,k,N,R1,L1) local ra1,ra2,F,lu,g,ope,INI,k1,n1:
ra1:=rand(1..R1):
ra2:=rand(1..L1):

F:=binomial(n,k):

lu:=ra1()*n+ra1()*k+ra1():
F:=F*binomial(lu,k):

g:=ra2():
F:=convert(F,factorial)*g^k:


ope:=zeil(F,k,n,N):

if ope=FAIL then
  RETURN(FAIL):
fi:

ope:=ope[1]:

INI:=[seq(add(eval(subs({n=n1,k=k1},F)),k1=0..n1),n1=1..degree(ope,N))]:

[F,ope,INI]:

end:




#ChIIIpre(n,k,N,R1,L1,KAMA): the data for Chapter III. KAMA challenging problems about proving that a given 
#sequence defined by a linear recurrence equation with polynomial coefficients with initial conditions all produce integers
#ChIIIpre(n,k,N,5,10,10);
ChIIIpre:=proc(n,k,N,R1,L1,KAMA)
local lu,gu:

gu:={}:

while nops(gu)<KAMA do
lu:=RandBCS(n,k,N,R1,L1):
 if lu<>FAIL then
  gu:=gu union {lu}:
 fi:
od:
gu:

end:




#ChIII(R1,L1,KAMA): Inputs positive integers R1, L1, and KAMA and outputs a list of KAMA Chellenging problems
#of proving that solutions of certain linear recurrences with polynomial coefficients and initial conditions
#a list of k challenging  problems about polynomials being non-negative. Try: 
#ChIII(4,10,10);
ChIII:=proc(R1,L1,KAMA) local gu,n,N,INI,i,gu1,F,ope:
gu:=ChIIIpre(n,k,N,R1,L1,KAMA):
gu:=convert(gu,list):
print(``):
print(KAMA, `Challenging Problems about Surprising Integrality of Sequences Defined by Linear Recurrences with Polynomial coefficients`):
print(``):
for i from 1 to KAMA do
 ope:=gu[i][2]:
 INI:=gu[i][3]:
 print(`Problem Number`, i,`: Let `, a(n), `be the sequence of rational numbers satisfying the linear recurrence`):
 print(``):
 print(add(coeff(ope,N,i1)*a(n+i1),i1=0..degree(ope,N))=0 ):
print(``):
print(`subject to the initial conditions`):
print(``):
print(seq(a(i1)=INI[i1],i1=1..nops(INI))):
print(``):
print(`Prove that , surprisingly,  that for EVERY n, a(n) an integer!`):
print(``):
print(`------------------------------------------------------------------------`):
print(``):
od:
print(``):
print(`Solutions to the Problems`):
print(``):
for i from 1 to nops(gu) do
 print(``):
 print(`Proof of Problem`, i):
 print(``):
 F:=gu[i][1]:
 print(`It so happened that a(n) equals`):
 print(``):
 print(Sum(F,k=0..n)):
 print(``):
 print(`By the Zeilberger algorithm this sum satisfies the given recurrence and since the initial conditions match (check!) and the sum is always an integer, the result follows by induction.`):
 print(``):
print(``):
print(`------------------------------------------------------------------------`):
print(``):

print(``):
od:
print(``):
print(`----------------------------------`):

end:


####From qEKHAD
ezra1qE:=proc()
if args=NULL then
print(`The supporting procedures are: `):
print(` hafokh `):
else
ezra(args):
fi:

end:
 
ezraqE:=proc()
if args=NULL then
print(`q-EKHAD: A Maple package for proving terminating basic hypergemetric`):
print(` (alias q-binomial coeff.) and other q-identities`):
print(`Handles only single-summation`):
print(`For help with a specific prodecure, type "ezra(procedure_name);"`):
print(`Contains procedures: `):
print(`qzeil,qzeilpap`):
print(`Warning: q is a global variable`):
print(`The summand can be expressed in terms of qfac(k):=(q)_k,`):
print(`or the Gaussian polynomials (alias q-binomial coefficients)`):
print(`qbin(a,b):=qfac(a)/(qfac(b)*qfac(b-a)), `):
print(`or qrf(q^(something),k):=(q^(something))_k, `):
print(`where the "something" must be affine linear in its variables`):
fi:



if nops([args])=1 and op(1,[args])=`hafokh` then
print(`hafokh(EXPRESSION,k,n,q,t,s): converts a q-expression in qfac(a*n+b*k+c) and rf(a,k)`):
print(`to the format required in procedure qctp`):
print(`hafokh  converts a q-expression involving powers of qfac(a*n+b*k+c)`):
print(`and qrf(a,k)'s to the format required in procedure qctp`):
print(`Try:`):
print(`hafokh(q^((5*k^2-k)/2)*(-1)^k*(1+q^k)/qfac(n+k)/qfac(n-k),k,n,q,t,s);`):
fi:
  
if nops([args])=1 and op(1,[args])=`qzeil` then
print(`qzeil(SUMMAND,N,k, n,para,MaxOrder)`):
print(`Inputs SUMMAND in terms of q-factorials of the form`):
print(`qfac(r):=(1-q)(1-q^2)...(1-q^r)`):
print(`qbin(a,b):=qfac(a)/(qfac(b)qfac(a-b),`):
print(` and the familiar q-rising factorial,`):
print(` (something)_k:=(1-something)*(1-q*something)*...*(1-q^(k-1)*something),`):
print(` which should be entered as  qrf(something,k).`):
print(` times q raised to a quadratic in k and n. `):
print(`It also inputs an positive integer, MaxOrder`):
print(`It outputs the recurrence operator ope(N,n), of order <.MxOrder, or FAIL, followed by`):
print(`the certificate R(n,k) such that  G(n,k):=F(n,k)*R(n,k)`):
print(`satisfies ope(N,n)F(n,k)=G(n,k)-G(n,k-1)`):
print(`All parameters are discrete, so if you have, say (a)_k, you must`):
print(`change a to q^(alpha), say. k is the summation variable, n is the`):
print(`the auxilliary variable, with respect to which the recurrence`):
print(`is given, and para is a list of all the other (necessarily) discrete`):
print(`variables present. On output, arbitrary constants are denoted by b1, b2, etc. `):
print(`Examples: The left side of the Finite-Form Rogers-Ramanujan is given`):
print(`by qzeil(q^(k^2)/qfac(k)/qfac(n-k),N,k,n,[],4);,`):
print(`( while the right side (after Peter Paule's symmetrization) is:`):
print(` qzeil(q^((5*k^2-k)/2)*(-1)^k*(1+q^k)/qfac(n+k)/qfac(n-k),N,k,n,[],4);`):
print(`Since the second-order recurrences outputted are identical (try it!)`):
print(`we have a proof of Rogers-Ramanujan.`):

print(`Another example: To prove q-Chu do:`):
print(` qzeil(q^(k^2)*qbin(n,k)*qbin(m,k)/qbin(n+m,m),N,k,n,[m],4);`):

print(`Another example: To prove q-Dixon (after the Paule symmetrization) do:`):
print(` qzeil((-1)^k*(1+q^k)*q^((3*k^2-k)/2)*qbin(2*n,n+k)^3/qfac(3*n)*qfac(n)^3,N,k,n,[],4);`):

print(`and with three parameters`):
print(`qzeil((-1)^k*(1+q^k)*q^((3*k^2-k)/2)*qbin(a+b,a+k)*qbin(a+c,c+k)*qbin(b+c,b+k)/qfac(a+b+c)*qfac(a)*qfac(b)*qfac(c),A,k,a,[b,c]);`):

print(`The syntax, once more, is: qzeil(SUMMAND,N,k, n,para,MaxOrder); `):
fi:
 
if nops([args])=1 and op(1,[args])=`qzeilpap` then
print(` qzeilpap(SUMMAND,k,n,para) or qzeilpap(SUMMAND,k,n,para,NAME)`):
print(`Just like qzeil(SUMMAND,N,k,n,para) but writes a paper with the proof`):
print(`NAME is optional`):
print(`Examples: The left side of the Finite-Form Rogers-Ramanujan is given`):
print(`by qzeilpap(q^(k^2)/qfac(k)/qfac(n-k),k,n,[]);`):
print(`( while the right side (after Peter Paule's symmetrization) is:`):
print(` qzeilpap(q^((5*k^2-k)/2)*(-1)^k*(1+q^k)/qfac(n+k)/qfac(n-k),k,n,[]);`):
print(`Since the second-order recurrences outputted are identical (try it!)`):
print(`we have a proof of the First Rogers-Ramanujan identity.`):
print(`For the two above examples, with the nake indicated, try:`):
print(` qzeilpap(q^(k^2)/qfac(k)/qfac(n-k),k,n,[],"Left Side of Finite Rogers-Ramanujan"); `);
print(``):
print(` qzeilpap(q^((5*k^2-k)/2)*(-1)^k*(1+q^k)/qfac(n+k)/qfac(n-k),k,n,[], "Right Side of Finite Rogers-Ramanujan" );`):


print(`Another example: To prove q-Chu do:`):
print(``):
print(` qzeilpap(q^(k^2)*qbin(n,k)*qbin(m,k)/qbin(n+m,m),k,n,[m]);`):
print(``):
print(`and with the name`):
print(``):
print(` qzeilpap(q^(k^2)*qbin(n,k)*qbin(m,k)/qbin(n+m,m),k,n,[m], "Chu-Vandermonde" );`):
print(``):
print(`Another example: for q-Dixon (after the Paule-symmetrization)`):
print(``):
print(` qzeilpap((-1)^k*(1+q^k)*q^((3*k^2-k)/2)*qbin(2*n,n+k)^3/qfac(3*n)*qfac(n)^3,k,n,[]);`):


print(`Yet another example: for q-Dixon with three parameters(after the Paule-symmetrization)`):
print(``):
print(`qzeilpap((-1)^k*(1+q^k)*q^((3*k^2-k)/2)*qbin(a+b,a+k)*qbin(a+c,c+k)*qbin(b+c,b+k)/qfac(a+b+c)*qfac(a)*qfac(b)*qfac(c),k,a,[b,c]);`):
fi:




 
end:
 
 
 



qpaper4:=proc(SUMMAND,k,n,para) local SHEM,IDENTITY,RECURRENCE, gu, N,ope,SDZ1:


gu:=qzeil(SUMMAND,N,k, n,para,6):

if gu=FAIL then 
  RETURN(FAIL):
fi:

ope:=gu[1]:
SDZ1:=gu[2]:

if degree(ope,N)=1 then
SHEM:=IDENTITY:
else
SHEM:=RECURRENCE:
fi:

print(``):
print(`A PROOF OF a q-Series `, SHEM):
lprint(``):
print(`by Shalosh B. Ekhad, Rutgers University`):
lprint(``):
print(`I will give a short proof of the following result.`):
print(``):
if degree(ope,N)=1 then
print(`(Note that since the recurrence below is first order, this`):
print(`means that the sum`, SUM(n), `has closed form,and it is`):
print(`easily seen to be equivalent.)`):
print(``):
fi:
print(`Theorem: Let `, F(n,k), ` be given by `):
lprint(``):
print(SUMMAND):
print(``):
print(`and let`, SUM(n),` be the sum of `, F(n,k),` with respect to `, k):
print(``):
print(SUM(n),` satisfies the following linear recurrence equation`):
print(``):
print(yafe(ope,N,n,SUM(n))=0):
print(``):
print(`and in Maple format`):
print(``):
lprint(yafe(ope,N,n,SUM(n))=0):
print(``):
print(`PROOF: We cleverly construct`, G(n,k),`:=`):
print(``):
print(SDZ1*SUMMAND):
print(``):
print(`and in Maple format `):
print(``):
lprint(SDZ1*SUMMAND):
print(``):
lprint(`with the motive that`):
print(``):
print(yafe(ope,N,n,F(n,k))=G(n,k)-G(n,k-1) ):
print(`(check!)`):
print(``):
print(`and the theorem follows upon summing with respect to `, k, ` . `):
end:
 


qpaper5:=proc(SUMMAND,k,n,para,NAME) local SHEM,IDENTITY,RECURRENCE, gu, N,ope,SDZ1:


gu:=qzeil(SUMMAND,N,k, n,para,4):

if gu=FAIL then 
  RETURN(FAIL):
fi:

ope:=gu[1]:
SDZ1:=gu[2]:

if degree(ope,N)=1 then
SHEM:=IDENTITY:
else
SHEM:=RECURRENCE:
fi:

print(``):
print(`A PROOF OF The`, NAME , ` q-Series  `, SHEM):
lprint(``):
print(`by Shalosh B. Ekhad, Rutgers University`):
lprint(``):
print(`I will give a short proof of the following result.`):
print(``):
if degree(ope,N)=1 then
print(`(Note that since the recurrence below is first order, this`):
print(`means that the sum`, SUM(n), `has closed form,and it is`):
print(`easily seen to be equivalent.)`):
print(``):
fi:
print(`Theorem: Let `, F(n,k), ` be given by `):
lprint(``):
print(SUMMAND):
print(``):
print(`and let`, SUM(n),` be the sum of `, F(n,k),` with respect to `, k):
print(``):
print(SUM(n),` satisfies the following linear recurrence equation`):
print(``):
print(yafe(ope,N,n,SUM(n))=0):
print(``):
print(`and in Maple format`):
print(``):
lprint(yafe(ope,N,n,SUM(n))=0):
print(``):
print(`PROOF: We cleverly construct`, G(n,k),`:=`):
print(``):
print(SDZ1*SUMMAND):
print(``):
print(`and in Maple format `):
print(``):
lprint(SDZ1*SUMMAND):
print(``):
lprint(`with the motive that`):
print(``):
print(yafe(ope,N,n,F(n,k))=G(n,k)-G(n,k-1) ):
print(`(check!)`):
print(``):
print(`and the theorem follows upon summing with respect to `, k, ` . `):
end:
 


 




 
qzeilpap:=proc(SUMMAND,k,n,para,NAME):
 
if not (nargs=4 or nargs=5) then
 ERROR(` qzeilpap(SUMMAND,k,n,para) or qzeilpap(SUMMAND,k,n,para,NAME)`):
 fi:
 
if nargs=5 then
qpaper5(SUMMAND,k,n,para,NAME):
else
qpaper4(SUMMAND,k,n,para):
fi:
end:
 
 
 
#qrf:converts (q^(linear expressions of variables and parameters))_k
#linear expression in parameters to an expression in qfac
 
qrf:=proc(qal,k)
local qal1:
 
qal1:=simplify(qal):
 
if qal=q then
 RETURN(qfac(k)):
fi:
 
if not type(qal1,`^`) then
 ERROR(`Argument must be a power of q`):
fi:
 
if not op(1,qal1)=q then
  ERROR(`Argument must be a power of q`):
fi:
 
qfac(op(2,qal1)+k-1)/qfac(op(2,qal1)-1):
end:
 
pashet:=proc(p,N)
local i,gu,p1:
p1:=expand(p):
 
gu:=0:
for i from 0 to degree(p1,N) do
gu:=gu+factor(coeff(p1,N,i))*N^i:
od:
RETURN(gu):
end:
 
 
ratn:=proc(alpha,khe,q,s)
local i1,lu,dg:
dg:=degree(alpha,s)*khe:
if dg=0 then
RETURN(1):
fi:
if dg>0 then
lu:=product(1-alpha*q^i1,i1=1..dg):
fi:
if dg <0 then
lu:=1/product(1-alpha/q^i1,i1=0..(-dg)-1):
fi:
lu:=normal(lu):
RETURN(lu):
end:
 
 
ratk:=proc(alpha,q,t)
local lu,dg,i1:
dg:=degree(alpha,t):
if dg=0 then
RETURN(1):
elif
dg > 0 then
lu:=product(1-alpha/q^i1,i1=0..dg-1):
else
lu:=1/product(1-alpha*q^i1,i1=1..(-dg)):
fi:
lu:=expand(lu):
lu:=normal(lu):
RETURN(lu):
end:
 
 
 
qctp:= proc(POL,NUM,DEN,ARG,QPOWER1,ORDER1,N, k, n, q, s, t,para)
local  ARG1,  DEN1,  KAMA, 
      KHE, L, L1, LU, LU1, NUM1, ORDER, P, Q, QPOWER, R, RAT, 
      RATB,  S, degg, eq, fu, g, gorem, gu, 
      gugu, i, ia1, j1, ja1, kak, khe, l1, l2, lu1a, meka1, mekb1, 
      mumu, va1, va2, y1, yakhas,k1,j2,kapi1,kapi2,misht,va11,VAK:
KAMA:=5;
NUM1:=NUM:DEN1:=DEN:ARG1:=ARG:QPOWER:=QPOWER1:
ORDER:=ORDER1:
gorem:=0:
for i from 0 to ORDER do
gorem:=gorem+b[i]*N^i:
od:
 
 
 
LU:=0:
KHE:=degree(gorem,N):
 
for khe from 0 to KHE do
RAT:=product('ratn(NUM1[y1],khe,q,s)','y1'=1..nops(NUM1)):
RATB:=product('ratn(DEN1[y1],khe,q,s)','y1'=1..nops(DEN1)):
RAT:=RAT/RATB:
kak:=subs(n=n+khe,QPOWER)-QPOWER:
kak:=expand(kak):
gu:=t^coeff(kak,k,1)*s^coeff(kak,n,1)*q^(kak-k*coeff(kak,k,1)-
n*coeff(kak,n,1)):
RAT:=RAT*gu*subs(s=s*q^khe,POL):
RAT:=normal(RAT):
LU:=coeff(gorem,N,khe)*RAT+LU:
LU:=normal(LU):
od:
lu1a:=LU:
LU1:=numer(LU):
LU:=1/denom(LU):
LU:=LU/subs(t=t/q,LU):
LU:=normal(LU):
yakhas:=product('ratk(NUM1[y1],q,t)','y1'=1..nops(NUM1))/
product('ratk(DEN1[y1],q,t)','y1'=1..nops(DEN1)):
yakhas:=normal(yakhas):
kak:=QPOWER-subs(k=k-1,QPOWER):
kak:=expand(kak):
kak:=t^coeff(kak,k,1)*s^coeff(kak,n,1)*q^(kak-k*coeff(kak,k,1)-
n*coeff(kak,n,1)):
yakhas:=ARG1*LU*yakhas*kak:
yakhas:=normal(yakhas):
P:=LU1:
Q:=numer(yakhas):
R:=denom(yakhas):
Q:=expand(Q):
R:=expand(R):
j1:=0:
while j1 <= KAMA do
kapi1:=numer(Q):
kapi2:=numer(subs(t=t*q^j1,R)):
kapi1:=expand(kapi1):
kapi2:=expand(kapi2):
 
for j2 from 1 to nops(para) do
 kapi1:=subs(q^op(j2,para)=misht[j2],kapi1):
 kapi2:=subs(q^op(j2,para)=misht[j2],kapi2):
od:
 
g:=gcd(kapi1,kapi2):
if g <> 1 then
 
for j2 from 1 to nops(para) do
 g:=subs(misht[j2]=q^op(j2,para),g):
od:
 
Q:=normal(Q/g):
R:=normal(R/subs(t=t/q^j1,g)):
P:=P*product(subs(t=t/q^ja1,g) , ja1=0..j1-1):
j1:=-1:
fi:
j1:=j1+1:
od:
P:=expand(P):
R:=expand(R):
Q:=expand(Q):
if(coeff(Q,t,0)*coeff(R,t,0)=0) then
L1:=0:
else
L:=ln(normal(coeff(Q,t,0)/coeff(R,t,0)))/ln(q):
if type(L,integer) then
L1:=max(0,L):
else
L1:=0:
fi:
fi:
P:=expand(t^L1*P):
R:=expand(q^L1*R):
l1:=degree(R,t):
l2:=degree(Q,t):
meka1:=coeff(R,t,l1):
mekb1:=coeff(Q,t,l2):
if l1 <>l2 
then
k1:=degree(P,t)-max(l1,l2):
else
 mumu:= meka1/(mekb1*q^l1):
mumu:=ln(mumu)/ln(q):
  if type(mumu,integer) 
  then
   k1:=max(mumu, degree(P,t)-l1):
  else
   k1:=degree(P,t)-l1:
  fi:
fi:
fu:=0:
if k1 < 0 then
RETURN(0):
fi:
for ia1 from 0 to k1 do
fu:=fu+a[ia1]*t^ia1:
od:
gugu:=subs(t=t*q,Q)*fu-R*subs(t=t/q,fu)-P:
gugu:=expand(gugu):
degg:=degree(gugu,t):
for ia1 from 0 to degg do
eq[ia1+1]:=coeff(gugu,t,ia1)=0:
od:
 
for ia1 from 0 to k1 do
va1[ia1+1]:=a[ia1]:
od:
for ia1 from 0 to ORDER do
va2[ia1+1]:=b[ia1]:
od:
eq:=convert(eq,set):
va1:=convert(va1,set):
va2:=convert(va2,set):
va1:=va1 union va2 :
va1:=solve(eq,va1):

VAK:={}:

for va11 in va1 do
if op(1,va11)=op(2,va11) then
 VAK:=VAK union {op(1,va11)}:
fi:
od:

if VAK={} then
  RETURN(FAIL):
fi:

fu:=subs(va1,fu):

fu:=subs({seq(va11=1,va11 in VAK)},fu):

gorem:=subs(va1,gorem):

gorem:=subs({seq(va11=1,va11 in VAK)},gorem):

fu:=normal(fu):
 
S:=subs(t=t*q,Q)*fu*lu1a/P :
S:=normal(S):
S:=factor(S):


gorem,S:
end:
 
qctpnis:= proc(POL,NUM,DEN,ARG,QPOWER1,ORDER1,N, k, n, q, s, t,para)
local  ARG1,  DEN1,  KAMA, 
      KHE, L, L1, LU, LU1, NUM1, ORDER, P, Q, QPOWER, R, RAT, 
      RATB,  S, degg, eq, fu, g, gorem, gu, 
      gugu, i, ia1, j1, ja1, kak, khe, l1, l2, lu1a, meka1, mekb1, 
      mumu, va1, va2, y1, yakhas,k1,j2,kapi1,kapi2,misht,a,b,va1q,goremq,fuq:
KAMA:=5;
NUM1:=NUM:DEN1:=DEN:ARG1:=ARG:QPOWER:=QPOWER1:
ORDER:=ORDER1:
gorem:=0:
for i from 0 to ORDER do
gorem:=gorem+(b[i])*N^i:
od:
 
 
 
LU:=0:
KHE:=degree(gorem,N):
 
for khe from 0 to KHE do
RAT:=product('ratn(NUM1[y1],khe,q,s)','y1'=1..nops(NUM1)):
RATB:=product('ratn(DEN1[y1],khe,q,s)','y1'=1..nops(DEN1)):
RAT:=RAT/RATB:
kak:=subs(n=n+khe,QPOWER)-QPOWER:
kak:=expand(kak):
gu:=t^coeff(kak,k,1)*s^coeff(kak,n,1)*q^(kak-k*coeff(kak,k,1)-
n*coeff(kak,n,1)):
RAT:=RAT*gu*subs(s=s*q^khe,POL):
RAT:=normal(RAT):
LU:=coeff(gorem,N,khe)*RAT+LU:
LU:=normal(LU):
od:
lu1a:=LU:
LU1:=numer(LU):
LU:=1/denom(LU):
LU:=LU/subs(t=t/q,LU):
LU:=normal(LU):
yakhas:=product('ratk(NUM1[y1],q,t)','y1'=1..nops(NUM1))/
product('ratk(DEN1[y1],q,t)','y1'=1..nops(DEN1)):
yakhas:=normal(yakhas):
kak:=QPOWER-subs(k=k-1,QPOWER):
kak:=expand(kak):
kak:=t^coeff(kak,k,1)*s^coeff(kak,n,1)*q^(kak-k*coeff(kak,k,1)-
n*coeff(kak,n,1)):
yakhas:=ARG1*LU*yakhas*kak:
yakhas:=normal(yakhas):
P:=LU1:
Q:=numer(yakhas):
R:=denom(yakhas):
Q:=expand(Q):
R:=expand(R):
j1:=0:
while j1 <= KAMA do

kapi1:=numer(Q):
kapi2:=numer(subs(t=t*q^j1,R)):
kapi1:=expand(kapi1):
kapi2:=expand(kapi2):
 
for j2 from 1 to nops(para) do
 kapi1:=subs(q^op(j2,para)=misht[j2],kapi1):
 kapi2:=subs(q^op(j2,para)=misht[j2],kapi2):
od:
 
g:=gcd(kapi1,kapi2):
if g <> 1 then
 
for j2 from 1 to nops(para) do
 g:=subs(misht[j2]=q^op(j2,para),g):
od:
 
Q:=normal(Q/g):
R:=normal(R/subs(t=t/q^j1,g)):
P:=P*product(subs(t=t/q^ja1,g) , ja1=0..j1-1):
j1:=-1:
fi:
j1:=j1+1:
od:
P:=expand(P):
R:=expand(R):
Q:=expand(Q):
if(coeff(Q,t,0)*coeff(R,t,0)=0) then
L1:=0:
else
L:=ln(normal(coeff(Q,t,0)/coeff(R,t,0)))/ln(q):
if type(L,integer) then
L1:=max(0,L):
else
L1:=0:
fi:
fi:
P:=expand(t^L1*P):
R:=expand(q^L1*R):
l1:=degree(R,t):
l2:=degree(Q,t):
meka1:=coeff(R,t,l1):
mekb1:=coeff(Q,t,l2):
if l1 <>l2 
then
k1:=degree(P,t)-max(l1,l2):
else
 mumu:= meka1/(mekb1*q^l1):
mumu:=ln(mumu)/ln(q):
  if type(mumu,integer) 
  then
   k1:=max(mumu, degree(P,t)-l1):
  else
   k1:=degree(P,t)-l1:
  fi:
fi:
fu:=0:
if k1 < 0 then
RETURN(0):
fi:
for ia1 from 0 to k1 do
fu:=fu+a[ia1]*t^ia1:
od:
gugu:=subs(t=t*q,Q)*fu-R*subs(t=t/q,fu)-P:
gugu:=expand(gugu):
degg:=degree(gugu,t):
for ia1 from 0 to degg do
eq[ia1+1]:=coeff(gugu,t,ia1)=0:
od:
 
for ia1 from 0 to k1 do
va1[ia1+1]:=a[ia1]:
od:
for ia1 from 0 to ORDER do
va2[ia1+1]:=b[ia1]:
od:
eq:=convert(eq,set):
va1:=convert(va1,set):
va2:=convert(va2,set):
va1:=va1 union va2 :
 
 
for j2 from 1 to nops(para) do
eq:=subs(misht[j2]=j2^2+1,eq):
od:
 


va1q:=solve(subs(q=2,eq),va1):

goremq:=subs(va1q,gorem):

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

va1:=solve(eq,va1):
fu:=subs(va1,fu):
gorem:=subs(va1,gorem):
if gorem=0 then
 RETURN(0):
 else
 RETURN(1):
fi:
 
end:
 
#qbin converts the q-binomial coefficient to q-factorials
 
qbin:=proc(a,b):
qfac(a)/qfac(b)/qfac(a-b):
end:
 
qbinomial1:=proc(a,b) local i:
if b=0 then
 RETURN(1):
fi:
if b<0 then
 RETURN(0):
fi:
expand(normal(mul(1-q^(a-i),i=0..b-1)/mul(1-q^i,i=1..b))):
end:

hafokhpol:=proc(POL,n,k,q,t,s):
subs({q^k=t,q^n=s},expand(POL)):
end:
 
 
#hafokh(EXPRESSION,k,n,q,t,s): converts a q-expression in qfac(a*n+b*k+c) and rf(a,k)
#to the format required in procedure qctp
#hafokh  converts a q-expression involving powers of qfac(a*n+b*k+c)
#and qrf(a,k)'s to the format required in procedure qctp
#Try:
#hafokh(q^((5*k^2-k)/2)*(-1)^k*(1+q^k)/qfac(n+k)/qfac(n-k),k,n,q,t,s);
 
hafokh:=proc(bitui1,k,n,q,t,s)
local gu,bitui,POL,NUM,DEN,ARG,QPOWER1,gu1,khezka,gu11,gu11n,gu11k,gu110,
i,r:
 
bitui:=factor(bitui1):
 
NUM:=[]:
DEN:=[]:
ARG:=1:
POL:=1:
QPOWER1:=0:
 
if nops(bitui)>1 then
for i from 1 to nops(bitui) do
 gu:=op(i,bitui):
 
  if type(gu,`+`) then
     POL:=POL*gu:
   POL:=hafokhpol(POL,n,k,q,t,s):
   next:
   fi:
 
 if type(gu,`^`) then 
  gu1:=op(1,gu):
  khezka:=op(2,gu):
 
    if gu1=q then
     if not  degree(khezka,k)<=2  then
       ERROR(q,`can be only raised to a power a quadratic in`,k):
     fi:
     QPOWER1:=QPOWER1+khezka:
     next:
    fi:
 
    if not (type(normal(khezka/k),integer) or type(khezka,integer)) then
     ERROR(`Wrong exponent`):
    fi:
 
    if not type(khezka, integer) then
      if not type(khezka/k ,integer) then
       ERROR(`Exponents must be in`,k):
      else
       ARG:=ARG*gu1^(khezka/k):
      fi:
      next:
    fi:
 
    if type(khezka,integer) then
      if type(gu1,function) then
        if op(0,gu1)<>'qfac' then
         ERROR(`The only functions allowed are qfac and qrf and qbin`):
        fi:
     if op(0,gu1)='qfac' then
      gu11:=op(1,gu1):
 
      gu11n:=coeff(gu11,n,1):
      gu11k:=coeff(gu11,k,1):
      gu110:=expand(gu11-gu11n*n-gu11k*k):
 
      if not type(gu11n,integer) or not type(gu11k,integer)
                    then
       ERROR(`Arguments to qfac must be affine-linear expressions in vars`):
       fi:
 
       gu11:=t^gu11k*s^gu11n*q^gu110:             
      fi:
 
        if khezka>0 then
          for r from 1 to khezka do
           NUM:=[op(NUM),gu11]:
          od:
        else
         for r from 1 to -khezka do
           DEN:=[op(DEN),gu11]:
         od:
        fi: 
     fi:
    next:
  fi:
 
fi:
 
   if type(gu,function) then
      if op(0,gu)<>'qfac'  then
        ERROR(`The only functions allowed are qfac and qrf and qbin`):
      fi:
 
  if  op(0,gu)='qfac' then
      gu11:=op(1,gu):
 
      gu11n:=coeff(gu11,n,1):
      gu11k:=coeff(gu11,k,1):
      gu110:=expand(gu11-gu11n*n-gu11k*k):
 
      if not type(gu11n,integer) or not type(gu11k,integer) then
       ERROR(`Arguments to qfac must be affine-linear expressions in vars`):
       fi:
 
      gu11:=t^gu11k*s^gu11n*q^gu110:             
       NUM:=[op(NUM),gu11]:
    fi:
 
 
   
    fi:
 
od:
 
fi:
 
 gu:=bitui:
 
  if type(gu,`+`) then
     POL:=POL*gu:
   POL:=hafokhpol(POL,n,k,q,t,s):
      fi:
 
 if type(gu,`^`) then 
  gu1:=op(1,gu):
  khezka:=op(2,gu):
 
    if gu1=q then
     if not  degree(khezka,k)<=2  then
       ERROR(q,`can be only raised to a power a quadratic in`,k):
     fi:
     QPOWER1:=QPOWER1+khezka:
     
    fi:
 
    if not (type(normal(khezka/k),integer) or type(khezka,integer)) then
     ERROR(`Wrong exponent`):
    fi:
 
    if not type(khezka, integer) then
      if not type(khezka/k ,integer) then
       ERROR(`Exponents must be in`,k):
      else
       ARG:=ARG*gu1^(khezka/k):
      fi:
          fi:
 
    if type(khezka,integer) then
      if type(gu1,function) then
 
        if op(0,gu1)<>'qfac' then
         ERROR(`The only function allowed are qfac and qrf`):
        fi:
 
      if op(0,gu1)='qfac' then
       gu11:=op(1,gu1):
 
      gu11n:=coeff(gu11,n,1):
      gu11k:=coeff(gu11,k,1):
      gu110:=expand(gu11-gu11n*n-gu11k*k):
 
      if not type(gu11n,integer) or not type(gu11k,integer)  then
       ERROR(`Arguments to qfac must be affine-linear expressions in vars`):
       fi:
 
       gu11:=t^gu11k*s^gu11n*q^gu110:             
     fi:
 
        if khezka>0 then
          for r from 1 to khezka do
           NUM:=[op(NUM),gu11]:
          od:
        else
         for r from 1 to -khezka do
           DEN:=[op(DEN),gu11]:
         od:
        fi: 
     fi:
  fi:
 
fi:
 
   if type(gu,function) then
      if op(0,gu)<>'qfac'  then
        ERROR(`The only functions allowed are qfac and qrf and qbin`):
      fi:
 
      if op(0,gu)='qfac' then
      gu11:=op(1,gu):
 
      gu11n:=coeff(gu11,n,1):
      gu11k:=coeff(gu11,k,1):
      gu110:=expand(gu11-gu11n*n-gu11k*k):
 
      if not type(gu11n,integer) or not type(gu11k,integer) then
       ERROR(`Arguments to qfac must be affine-linear expressions in vars`):
       fi:
 
      gu11:=t^gu11k*s^gu11n*q^gu110:             
     fi:
 
 
       NUM:=[op(NUM),gu11]:
 
   
    fi:
    
POL,NUM,DEN,ARG,QPOWER1:         
end:
 
 
qzeilsp:= proc(SUMMAND,ORDER1,N, k, n,para)
local s,t,kak,POL,NUM,DEN,ARG,QPOWER1,pip,pip1,pip2,i,coe,i1:
 
kak:=hafokh(SUMMAND,k,n,q,t,s):
POL:=kak[1]:NUM:=kak[2]:DEN:=kak[3]:ARG:=kak[4]:QPOWER1:=kak[5]:
pip:=qctpnis(POL,NUM,DEN,ARG,QPOWER1,ORDER1,N, k, n, q, s, t,para):
 
if pip=0 then
 RETURN(0):
fi:
pip:=qctp(POL,NUM,DEN,ARG,QPOWER1,ORDER1,N, k, n, q, s, t,para):
if pip=0 or [pip]=[0,0] then
 RETURN(0):
fi:
pip:=subs({s=q^n,t=q^k},pip[1]),subs({s=q^n,t=q^k},pip[2]):

pip1:=simplify(pip[1]):

pip1:=add(coeff(pip1,N,i)*N^i,i=ldegree(pip1,N)..degree(pip1,N)):

pip2:=expand(normal(simplify(pip[2]))):

coe:=lcoeff(pip1,N):

pip1:=add( factor(expand(coeff(pip1,N,i1)/coe))*N^i1,i1=ldegree(pip1,N)..degree(pip1,N)):

pip2:=factor(expand(pip2/coe)):

pip1,pip2:
end:
 
qzeil:= proc(SUMMAND,N,k, n,para,MaxOrder)
local ORDER1,ope,mu:
 
mu:=subs(n=n+1,SUMMAND)/SUMMAND:
mu:=normal(mu):
mu:=subs(k=k+1,mu)/mu:
mu:=normal(mu):
 
if mu=1 then
ERROR(`Summand is seperable`):
fi:
 
 
ope:=0:
 
for ORDER1 from 1 to MaxOrder while ope=0 do
ope:=qzeilsp(SUMMAND,ORDER1,N, k, n,para):
od:
 
if ope=0 then
  RETURN(FAIL):
else
 RETURN(ope):
fi:
 
end:
 
 
 
qfac1:=proc(r)
option remember:
 
if r=0 then
 1:
elif r<0 then
 0:
else
 expand(qfac1(r-1)*normal((1-q^r)/(1-q))):
fi:
end:
 
 
 
 
 
yafe:=proc(ope,N,n,SUM)
local gu,i:
gu:=0:
 
for i from 0 to degree(ope,N) do
 gu:=gu+coeff(ope,N,i)*subs(n=n+i,SUM):
od:
 
gu:
end:
 



####End from qEKHAD


#RandqBCS(n,k,N,R1,R2,R3,MaxOrder): a random q-binomial coefficients sum inspired by George Andrews' seminal
#work on the Lusztig-Macdonald-wall conjecture. Tweaking the summand
#qbinomial(2*n,n-j)*q^(j^2) to be 
#qbinomial(2*n,n-a*j)*q^(b*j*(j-1)/2+c*j) to be 
#where a between 1 and R1, b is between 1 and R2 and c is between -R3 and R3 and the recurrence has at most order MaxOrder.
#It returns FAIL if the ordrer is larger.
#The output is
#[Summand,RecurrenceOperator,InitialConditions,LimitAs n goes to infinity]
#Try:
#RandqBCS(n,j,N,2,5,5,3);

RandqBCS:=proc(n,j,N,R1,R2,R3,MaxOrder) local F1, ra1,ra2,ra3,po,F,ope,INI,j1,n1,PP,fu,lu:
ra1:=rand(1..R1):
ra2:=rand(1..R2):
ra3:=rand(1..R3):
po:=expand(ra2()*j*(j-1)/2+ra3()*j):
lu:=ra1():

F:=qbinomial(2*n,n-lu*j)*q^po:
F1:=eval(subs(qbinomial=qbin,F)):

PP:=1+q^(-2*coeff(po,j)):

ope:=qzeil(F1*(1+PP),N,j,n,[],MaxOrder):
if ope=FAIL then
  RETURN(FAIL):
fi:


ope:=ope[1]:

if degree(ope,N)=1 then
 RETURN(FAIL):
fi:

INI:=[seq(add(expand(subs({n=n1,j=j1},F)),j1=-trunc(n1/coeff(po,j,2))..trunc(n1/coeff(po,j,2))),n1=1..degree(ope,N))]:
INI:=expand(eval(subs(qbinomial=qbinomial1,INI))):

fu:=q^po:
[F,ope,INI,fu]:

end:


#ChIVpre(n,j,N,R1,R2,R3,MaxOrder,KAMA): the data for Chapter IV. KAMA challenging problems about proving that a given 
#sequence defined by a q-linear recurrence equation with polynomial coefficients with initial conditionn
#has a certain limit as n goes to infinity
#ChIVpre(n,k,N,2,3,5,3,10);
ChIVpre:=proc(n,j,N,R1,R2,R3,MaxOrder,KAMA)
local lu,gu:

gu:={}:

while nops(gu)<KAMA do
lu:=RandqBCS(n,j,N,R1,R2,R3,MaxOrder):
 if lu<>FAIL then
  gu:=gu union {lu}:
 fi:
od:
gu:

end:

#ChIV(R1,R2,R3,MaxOrder,KAMA): creates a challeging book for George Andrews, inspired by his beautiful approach to the L-M-W conjecture.
#Try:
#ChIV(2,3,4,3,4);
ChIV:=proc(n,j,N,R1,R2,R3,MaxOrder,KAMA) local gu,INI,fu,i1,F,i,ope:
gu:=ChIVpre(n,j,N,R1,R2,R3,MaxOrder,KAMA):

print(`gu is`, gu):

gu:=convert(gu,list):
print(``):
print(KAMA, `Challenging Problems about The limit as n goes to infinity to Solutions of q-Recurrence equations with polynomial coefficients in q^n`):
print(``):
for i from 1 to KAMA do
 ope:=gu[i][2]:
 INI:=gu[i][3]:
 fu:=gu[i][4]:
 print(`Problem Number`, i,`: Let `, a(n), `be the sequence of polynomials in q^n defined by the linear recurrence`):
 print(``):
 print(add(coeff(ope,N,i1)*a(n+i1),i1=0..degree(ope,N))=0 ):
print(``):
print(`subject to the initial conditions`):
print(``):
print(seq(a(i1)=INI[i1],i1=1..nops(INI))):
print(``):
print(`Prove that , surprisingly,  the limit of a(n) as n goes to infinity equals`):
print(``):
print(Sum(fu,j=infinity..infinity)/Product(1-q^i1,i1=1..infinity)):
print(`------------------------------------------------------------------------`):
print(``):
od:
print(``):
print(`Solutions to the Problems`):
print(``):
for i from 1 to nops(gu) do
 print(``):
 print(`Proof of Problem`, i):
 print(``):
 F:=gu[i][1]:
 print(`It so happened that a(n) equals`):
 print(``):
 print(Sum(F,j=-n..n)):
 print(``):
 print(`By the q-Zeilberger algorithm (pre-processed with the amazing Peter Paule symmetrization)`):
 print(`this sum satisfies the given recurrence and since the initial conditions match (check!) and the sum is always an integer, the result follows by induction.`):
 print(``):
 print(`Now take the limit as n goes to infinity and you are done.`):
print(``):
print(`------------------------------------------------------------------------`):
print(``):

print(``):
od:
print(``):
print(`----------------------------------`):

end: