#######################################################################
## EKHAD: Save this file as EKHAD. To use it, stay in the             #
## same directory, get into Maple (by typing: maple <Enter> )         #
## and then type:  read EKHAD : <Enter>                               #
## Then follow the instructions given there                           #
##                                                                    #
## Written by Doron Zeilberger, Rutgers University ,                  #
##  zeilberg@math.rutgers.edu.                                        # 
#######################################################################

with(SolveTools):
solve1:=Linear:

print(` Last update: April 26, 2018, thanks to Daniel G. DuParc`):
print(` Previous updates:   Dec. 11, 2015 (adding procedures AZdI, AZcI ); May 30, 2015 (adding procedure TerryTao )`):
print(` May 29, 2014 (adding Ekhad )`):
print(`Version of July 2003: adapted to Maple 8 and 9 `):
print(`Many thanks to Drew Sills`):
print(`In the penultimate Version of Feb 25, 1999 a suggestion`):
print(` of Frederic Chyzak was used, with considerable `):
print(`speed-up. We thank him SO MUCH!`):

print(``):
print(`The penpenultimate version, Feb. 1997,`):
print(` corrected a subtle bug discovered by Helmut Prodinger`):
print(``):
 
print(`Previous versions benefited from comments by Paula Cohen, `):
print(`Lyle Ramshaw, and Bob Sulanke.`):
print(``):
print(`This is EKHAD, One of the Maple packages`):
print(`accompanying the book `):
print(` "A=B" `):
print(` (published by A.K. Peters, Wellesley, 1996) `):
print( ` by Marko Petkovsek, Herb Wilf, and Doron Zeilberger.`):
print(``):
print(`The most current version is available on WWW at:`):
print(`http://www.math.rutgers.edu/~zeilberg .`):
print(`Information about the book, and how to order it, can be found in`):
print(`http://www.central.cis.upenn.edu/~wilf/AeqB.html .`):
print(`Please report all bugs to: zeilberg@math.rutgers.edu .`):
print(`All bugs or other comments used will be acknowledged in future`):
print(`versions.`):
print(``):
print(`For general help, and a list of the available functions,`):
print(` type "ezra();". For specific help type "ezra(procedure_name)" `):
print(``):
 
 
 
ezra:=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+2)*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 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:
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:=6:
 
  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
 gu1:=cttry(SUMMAND,ORDER,k,n,resh,N):
 if gu1=1 then
  gu:=ct(SUMMAND,ORDER,k,n,N):
 
   if gu<>0 then
     RETURN(gu):
   fi:
 fi:
od:
 
ERROR(`No recurrence of order<=`,MAXORDER,`was found`):
 
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
 gu1:=cttry(SUMMAND,ORDER,k,n,resh,N):
 if gu1=1 then
  gu:=ct(SUMMAND,ORDER,k,n,N):
 
   if gu<>0 then
     RETURN(gu):
   fi:
 fi:
od:
 
ERROR(`No recurrence of order<=`,MAXORDER,`was found`):
 
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
 gu1:=cttry(SUMMAND,ORDER,k,n,resh,N):
 if gu1=1 then
  gu:=ct(SUMMAND,ORDER,k,n,N):
 
   if gu<>0 then
     RETURN(gu):
   fi:
 fi:
od:
 
ERROR(`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: