print(`This is a copy of the package EKHAD accompanying the book A=B`): print(`by M. Petkovsek, H. Wilf an D. Zeilberger`): print(`Because it was used in Zeilberger's paper published in v.2`): print(`of the New York J. of Math, it is reproduced here for the`): print(`convenience of future readers, and to guarantee its permanence.`): print(`Zeilberger would like to thank Mark Steinberger for his suggestion`): print(`to install a permanent copy of the package in the NYJM's archives`): lprint(``): lprint(``): print(` Version of May 31, 1996`): lprint(``): lprint(``): print(`This is EKHAD, One of the Maple packages`): print(`accompanying the book `): print(` "A=B" `): print(` (published by A.K. Peters, Welesley, 1996) `): print( ` by Marko Petkovsek, Herb Wilf, and Doron Zeilberger.`): lprint(``): print(`The most current version is available on WWW at:`): print(` http://www.math.temple.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.temple.edu .`): print(`All bugs or other comments used will be acknowledged in future`): print(`versions.`): lprint(``): print(`Version of May 31, 1996`): lprint(``): print(`This version corrects a subtle bug discovered by Helmut Prodinger`): lprint(``): print(`Previous versions benefited from comments by Paula Cohen, `): print(`Lyle Ramshaw, and Bob Sulanke.`): lprint(``): print(`For general help, and a list of the available functions,`): print(` type "ezra();". For specific help type "ezra(procedure_name)" `): lprint(``): ezra:=proc() if args=NULL then print(`EKHAD`): print(`A Maple package for proving Hypergemetric (Binomial Coeff.)`): print(` and other kinds of identities`): print(`For help with a specific procedure, type "ezra(procedure_name);"`): print(`Contains procedures: `): print(`findrec,ct,zeil,zeilpap,zeillim,AZd,AZc,AZpapd,AZpapc,celine`): 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.`): lprint(``): 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 AZpapc(1/sqrt(1-2*x*t+t^2)/t^(n+1),t,x) gives the`): print(`the differential equation satisfied by the Legendre polynomials`): 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 AZpapd(1/sqrt(1-2*x*t+t^2)/t^(n+1),t,n) gives the`): print(`the recurrence, and proof, satisfied by the Legendre polynomials`): 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 integrale 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 AZd(1/sqrt(1-2*x*t+t^2)/t^(n+1),t,n,N) gives`): print(`the diff.eq., and proof certificate, for the Legendre polynomials.`): 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 <=8,`): 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 AZc(1/sqrt(1-2*x*t+t^2)/t^(n+1),t,x,D) gives the`): print(`the diff.eq., and proof certificate for the Legendre polynomials.`): 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])=`zeil` then print(`Like ct,`): 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(`But it is not necessary to guess the ORDER`): 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.`): lprint(``): 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: 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: 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:=resultant(Q,subs(k=k+j,R),k): kv:=RootOf1(res1,j): 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: 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:=solve(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: 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:=resultant(Q,subs(k=k+j,R),k): kv:=RootOf1(res1,j): 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: 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:=solve(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: lprint(``): lprint(`A PROOF OF THE`,NAME,SHEM): lprint(``): lprint(`By Shalosh B. Ekhad, Temple University, ekhad@math.temple.edu`): lprint(``): lprint(`I will give a short proof of the following result(`,REF,`).`): lprint(``): if degree(ope1,N)=1 then lprint(`(Note that since the recurrence below is first order, this`): lprint(`means that the sum`, SUM(n), `has closed form,and it is`): lprint(`easily seen to be equivalent.)`): lprint(``): fi: lprint(`Theorem:Let`, F(n,k), `be given by`): lprint(``): print(SUMMAND): lprint(``): lprint(`and let`, SUM(n),`be the sum of`, F(n,k),` with`): lprint(`respect to`, k,`.`): lprint(``): lprint(SUM(n),` satisfies the following linear recurrence equation`): lprint(``): print(yafe(ope1,N,n,SUM(n))): print(`=0.`): lprint(``): lprint(`PROOF: We cleverly construct`, G(n,k),`:=`): lprint(``): print(SDZ1*SUMMAND): lprint(`with the motive that`): lprint(``): print(yafe(ope1,N,n,F(n,k))): lprint(`=`,G(n,k+1)-G(n,k), `(check!)`): lprint(``): lprint(`and the theorem follows upon summing with respect to`, k,`.QED.`): end: paper3:=proc(SUMMAND,k,n,N,ope1,SDZ1): lprint(``): lprint(`A PROOF OF A RECURRENCE`): lprint(``): lprint(`By Shalosh B. Ekhad, Temple University, ekhad@math.temple.edu`): lprint(``): lprint(`Theorem:Let`, F(n,k), `be given by`): lprint(``): print(SUMMAND): lprint(``): lprint(`and let`, SUM(n),`be the sum of`, F(n,k),` with`): lprint(`respect to`, k,`.`): lprint(``): lprint(SUM(n),` satisfies the following linear recurrence equation`): lprint(``): print(yafe(ope1,N,n,SUM(n))): print(`=0.`): lprint(``): lprint(`PROOF: We cleverly construct`, G(n,k),`:=`): lprint(``): print(SDZ1*SUMMAND): lprint(`with the motive that`): lprint(``): print(yafe(ope1,N,n,F(n,k))): lprint(`=`,G(n,k+1)-G(n,k), `(check!)`): lprint(``): lprint(`and the theorem follows upon summing with respect to`, k,`.QED.`): end: paperc:=proc(INTEGRAND,y,x,D,ope1,SDZ1): lprint(``): lprint(``): lprint(`A PROOF OF A DIFFERENTIAL EQUATION SATISFIED BY AN INTEGRAL`): lprint(``): lprint(`By Shalosh B. Ekhad, Temple University, ekhad@math.temple.edu`): lprint(``): lprint(`I will give a short proof of the following result.`): lprint(``): lprint(`Theorem:Let`, F(x,y), `be given by`): lprint(``): print(INTEGRAND): lprint(``): lprint(`and let`, INTEGRAL(x),`be the integral of`, F(x,y),` with`): lprint(`respect to`, y,`.`): lprint(``): lprint(INTEGRAL(x),` satisfies the following linear differential equation`): lprint(``): print(yafec(ope1,D,x,INTEGRAL(x))): print(`=0.`): lprint(``): lprint(`PROOF: We cleverly construct`, G(x,y),`:=`): lprint(``): print(SDZ1*INTEGRAND): lprint(`with the motive that`): lprint(``): print(yafec(ope1,D,x,F(x,y))): lprint(`=`,diff(G(x,y),y), `(check!)`): lprint(``): lprint(`and the theorem follows upon integrating with respect to`, y,`.`): end: paperd:=proc(INTEGRAND,x,n,N,ope1,SDZ1): lprint(``): lprint(`A PROOF OF A LINEAR RECURRENCE SATISFIED BY AN INTEGRAL`): lprint(``): lprint(`By Shalosh B. Ekhad, Temple University, ekhad@math.temple.edu`): lprint(``): lprint(`I will give a short proof of the following result.`): lprint(``): lprint(`Theorem:Let`, F(n,x), `be given by`): lprint(``): print(INTEGRAND): lprint(``): lprint(`and let`, INTEGRAL(n),`be the integral of`, F(n,x),` with`): lprint(`respect to`, x,`.`): lprint(``): lprint(INTEGRAL(n),` satisfies the following linear recurrence equation`): lprint(``): print(yafe(ope1,N,n,INTEGRAL(n))): print(`=0.`): lprint(``): lprint(`PROOF: We cleverly construct`, G(n,x),`:=`): lprint(``): print(SDZ1*INTEGRAND): lprint(`with the motive that`): lprint(``): print(yafe(ope1,N,n,F(n,x))): lprint(`=`,diff(G(n,x),x)): lprint(``): lprint(`and the theorem follows upon integrating with respect to`, x,`.QED.`): end: paperc:=proc(INTEGRAND,y,x,D,ope1,SDZ1): lprint(`A PROOF OF A DIFFERENTIAL EQUATION SATISFIED BY AN INTEGRAL`): lprint(``): lprint(`By Shalosh B. Ekhad, Temple University, ekhad@math.temple.edu`): lprint(``): lprint(`I will give a short proof of the following result.`): lprint(``): lprint(`Theorem:Let`, F(x,y), `be given by`): lprint(``): print(INTEGRAND): lprint(``): lprint(`and let`, INTEGRAL(x),`be the integral of`, F(x,y),` with`): lprint(`respect to`, y,`.`): lprint(``): lprint(INTEGRAL(x),` satisfies the following linear differential equation`): lprint(``): print(yafec(ope1,D,x,INTEGRAL(x))): print(`=0.`): lprint(``): lprint(`PROOF: We cleverly construct`, G(x,y),`:=`): lprint(``): print(SDZ1*INTEGRAND): lprint(`with the motive that`): lprint(``): print(yafec(ope1,D,x,F(x,y))): lprint(`=`,diff(G(x,y),y), `(check!)`): lprint(``): lprint(`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: 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:=8: for ORDER from 0 to KAMA do gu:=duisc(A,ORDER,y,x,D): if gu<>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:=solve(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:=solve(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:=solve(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:=solve({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':lprint(` ` ):lprint(` `):lprint(`The full recurrence is`): lprint(` `): 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: lprint(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))):lprint(` ` ):lprint(` `): lprint(`The telescoped form is`);lprint(` `); zy:=factor(subs(yy,sum(simplify(sum(b[l*(jj+1)+m],l=0..ii))*F(n-m,k),m=0..jj))): lprint(zy,`==G(n,k)-G(n,k-1)`); lprint(` `); lprint(`where G(n,k)=R(n,k)*F(n,k) and the rational function R(n,k) is `): lprint(` `); rr; end: