#solve1:=readlib(`solve/linear`): solve1:=solve: print(`This is AppsWZ, one of the two Maple packages`): print(`accompanying the article `): print(`Five Enumerative and Probabilstic Applications `): print(` of Wilf-Zeilberger Theory `): print( ` by Moa Apagodu and Doron Zeilberger.`): print(`First posted: Oct. 18, 2006 `): lprint(``): print(`The most current version is available on WWW at:`): print(` http://www.math.rutgers.edu/~zeilberg .`): print(`Please report all bugs to: zeilberg at math dot rutgers dot edu .`): lprint(``): print(`For general help, and a list of the available functions,`): print(` type "ezra();". For specific help type "ezra(procedure_name)" `): lprint(``): #########Begin EKHAD############# ezraEkhad:=proc() if args=NULL then print(`EKHAD`): print(`A Maple package for proving Hypergemetric (Binomial Coeff.)`): print(` and other kinds of identities`): print(`This version (Feb, 25, 1999) is much faster than the previous`): print(`version, thanks to a SLIGHT (yet POWERFUL) modification suggested by`): print(` FREDERIC CHYZAK `): lprint(``): 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 polct:ynomials.`): 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(`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.`): 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: 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: 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,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:=100: 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:=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':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: 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: #########End EKHAD############# ###Begin New Stuff################### ezra:=proc() if args=NULL then print(`AppsWZ`): lprint(``): print(`For help with a specific procedure, type "ezra(procedure_name);"`): print(`Contains procedures: Asy `): print(` AZdL , ExpProbFirstDie `): print(`FairDieStoryE, LoadedDiceEqs, RecProbVisit, RecProbVisitE `): print(` SeqFromRec, TeamEffortMoneyChanging `): elif nops([args])=1 and op(1,[args])=Asy then print(`Asy(ope,n,N,K): the asymptotic expansion of solutions `): print(`to ope(n,N)f(n)=0,where ope(n,N) is a recurrence operator`): print(`up to the K's term`): print(`For example, try: Asy((n+1)*N-(n+2),n,N,3);`): elif nops([args])=1 and op(1,[args])=AsyC then print(`AsyC(ope,n,N,K,Ini,L): the asymptotic expansion `): print(`complete with the estimated constant, of solutions `): print(`to ope(n,N)f(n)=0, with initial conditions`): print(`Ini, where ope(n,N) is a recurrence operator`): print(`up to the K's term. The constant is estimated using up to L terms`): print(`For example, try`): print(`AsyC(N^2-N-1,n,N,1,[1,1],100);`): elif nops([args])=1 and op(1,[args])=AZdL then print(` AZdL(A,y,n,N): an inhomogeneous recurrence in n `): print(`for int(A(y,n),y=a..b), using N as a shift operator`): print(`returns the operator, the right hand side `): print(`For example, try: `): print(`AZdL(x^n*(1-x)^m,x,n,N,0,1);`): elif nops([args])=1 and op(1,[args])=ExpProbFirstDie then print(`ExpProbFirstDie(Lp1,Lp2,a0): in a Hidden-Markov Model`): print(`with two dice (with the same set of outcomes)`): print(`with prob. dist Lp1 and Lp2, if the outcome`): print(`was given by the list a0 what is the "most reasonable"`): print(`(in the Bayesian-Laplace sense, NOT the ML sense) for the`): print(`prob. of the first die being chosen. For example, try:`): print(`ExpProbFirstDie([1/2,1/2],[1/3,2/3],[10,10]):`): elif nops([args])=1 and op(1,[args])=FairDieStoryE then print(`FairDieStoryE(n,a, Gvul):Gives the story of fair dice`): print(`with an even number of faces, up to 2*Gvul`): print(`2n is the number of steps and a is a symbol for a(n)`): print(`Try:`): print(`FairDieStoryE(n,a,2); `): elif nops([args])=1 and op(1,[args])=IntDist then print(`IntDist(Lp1,Lp2,a0): the Integral of the distribution w.r.t. x`): print(`induced from two loaded dice whose faces are 1, 2, .., nops(Lp1)`): print(`with probability distibution Lp1 and Lp2`): print(`For example, try: IntDist([1/2,1/2],[1/3,2/3],[2,2]);`): elif nops([args])=1 and op(1,[args])=IntDist1 then print(`IntDist1(Lp1,Lp2,a0,k): the Integral of the distribution w.r.t. x`): print(`induced from two loaded dice whose faces are 1, 2, .., nops(Lp1)`): print(`with probability distibution Lp1 and Lp2`): print(`times x^k`): print(`For example, try: IntDist1([1/2,1/2],[1/3,2/3],[2,2],1);`): elif nops([args])=1 and op(1,[args])=LoadedDiceEqs then print(`LoadedDiceEqs(Lp1,Lp2,a,A):The non-homogeneous recurrence`): print(`equations for each of the variables i (i=1.. nops(Lp1))`): print(`for the integral of the prob.`): print(`dist induced by two loaded dice with prob. dist.`): print(`in the form [ope,rhs]. For example, try:`): print(`LoadedDiceEqs([1/2,1/2],[1/3,2/3],a,A);`): elif nops([args])=1 and op(1,[args])=LoadedDiceEqs1 then print(`LoadedDiceEqs1(Lp1,Lp2,a,A):The non-homogeneous recurrence`): print(`equations for each of the variables i (i=1.. nops(Lp1))`): print(`for the integral of the prob. of x time `): print(`dist induced by two loaded dice with prob. dist.`): print(`in the form [ope,rhs]. For example, try:`): print(`LoadedDiceEqs1([1/2,1/2],[1/3,2/3],a,A);`): elif nops([args])=1 and op(1,[args])=LoadedDiceP then print(`LoadedDiceP(Lp1,Lp2,x,a): the probability that you get a[i]'s outcome`): print(`of i if you have two loaded dice whose faces are 1, 2, .., nops(Lp1)`): print(`with probability distibution Lp1 and Lp2`): print(`with the probability of choosing the first die being x`): print(`For example, try: LoadedDiceP([1/2,1/2],[1/3,2/3],x,a);`): elif nops([args])=1 and op(1,[args])=pgf then print(`pgf(L,x): Given a list of pairs [i,p] where the prob. of`): print(`of going from x to x+i is p, outputs the prob. generating`): print(`function w.r.t. x for one trial`): print(`For example, try: pgf( [[-1,1/2],[1,1/2]],x);`): elif nops([args])=1 and op(1,[args])=RecProbVisit then print(`RecProbVisit(L,n,N,k): if you roll a die with faces and prob. dist.`): print(`given by a list L of pairs [i,p], n times, let a_k(n) be the`): print(`probability that you would have exactly k dollars after n rolls`): print(`outputs a linear recurrence equation with polynomial coefficients`): print(`for a_k(n) followed by initial conditions`): print(`For example, try: RecProbVisit([[1,1/2],[-1,1/2]],n,N,0):`): elif nops([args])=1 and op(1,[args])=RecProbVisitE then print(`RecProbVisitE(L,n,N,k):`): print(` if you roll a die with faces and prob. dist.`): print(`given by a list L of pairs [i,p], n times, let a_k(n) be the`): print(`probability that you would have exactly k dollars after 2n rolls`): print(`outputs a linear recurrence equation with polynomial coefficients`): print(`for a_k(n) followed by initial conditions`): print(`For example, try: RecProbVisitE([[1,1/2],[-1,1/2]],n,N,0):`): elif nops([args])=1 and op(1,[args])=SeqFromRec then print(`SeqFromRec(ope,n,N,Ini,K): Given the first L-1`): print(`terms of the sequence Ini=[f(1), ..., f(L-1)]`): print(`satisfied by the recurrence ope(n,N)f(n)=0`): print(`extends it to the first K values`): elif nops([args])=1 and op(1,[args])=TeamEffortMoneyChanging then print(`TeamEffortMoneyChanging(S,n,N,r): Given a set of denominations`): print(`S and symbols n and r, outputs the linear`): print(`recurrence equation with polynomial coefficients(in n) of`): print(`A_r(n): the number of ways of coming up with n cents`): print(`using the denominations in S that can come from`): print(`r different people. For example try:`): print(`TeamEffortMoneyChanging({1,2},n,N,r);`): else print(`There is no such thing as`, args): fi: end: ezra1:=proc() if nargs=0 then print(`The supporting procedures are: `): print(` IntDist, IntDist1, LoadedDiceP, pgf `): else ezra(args): fi: end: #AZdL(A,y,n,N): an inhomogeneous recurrence in n #for int(A(y,n),y=a..b), using N as a shift operator #returns the operator, the right hand side #For example, try: #AZdL(x^n*(1-x)^m,x,n,N,0,1); AZdL:= proc(A,y,n,N,a,b) local ORDER,gu,KAMA: option remember: KAMA:=100: for ORDER from 0 to KAMA do gu:=duisdL(A,ORDER,y,n,N,a,b): if gu<>0 then RETURN(gu): fi: od: print(`No recurrence of order <=`, KAMA): FAIL: end: #duisdL(A,ORDER,y,n,N,a,b): an inhomogenous recurrence #for int(A(n), y=a..b): duisdL:= proc(A,ORDER,y,n,N,a,b) local gu,ope,rhs,cert: gu:=duisd(A,ORDER,y,n,N): if gu=0 then RETURN(0): fi: ope:=gu[1]: cert:=gu[2]: rhs:=normal(subs(y=b,cert*A)-subs(y=a,cert*A)): [ope,rhs]: end: #########Begin Bayesian Dice################### #LoadedDiceP(Lp1,Lp2,x,a): the probability that you get a[i]'s outcome #of i if you have two loaded dice whose faces are 1, 2, .., nops(Lp1) #with probability distibution Lp1 and Lp2 #with the probability of choosing the first die being x #For example, try: LoadedDiceP([1/2,1/2],[1/3,2/3],x,a); LoadedDiceP:=proc(Lp1,Lp2,x,a) local i: option remember: if nops(Lp1)<>nops(Lp2) or (type(convert(Lp1,`+`),numeric) and convert(Lp1,`+`)<>1 ) or (type(convert(Lp2,`+`),numeric) and convert(Lp2,`+`)<>1 ) then ERROR(`Bad input`): fi: mul((Lp1[i]*x+Lp2[i]*(1-x))^a[i],i=1..nops(Lp1)): end: #LoadedDiceEqs(Lp1,Lp2,a,A):The non-homogeneous recurrence #equations for each of the variables i (i=1.. nops(Lp1)) #in the form [ope,rhs] for the integral of the prob. #dist induced by two loaded dice with prob. dist. #Lp1,Lp2. For example, try: #LoadedDiceEqs([1/2,1/2],[1/3,2/3],a,A); LoadedDiceEqs:=proc(Lp1,Lp2,a,A) local i,gu,x: option remember: if nops(Lp1)<>nops(Lp2) or (type(convert(Lp1,`+`),numeric) and convert(Lp1,`+`)<>1 ) or (type(convert(Lp2,`+`),numeric) and convert(Lp2,`+`)<>1 ) then ERROR(`Bad input`): fi: gu:=LoadedDiceP(Lp1,Lp2,x,a): gu:=[seq(AZdL(gu,x,a[i],A[i],0,1),i=1..nops(Lp1))]: if member(FAIL,{op(gu)}) then RETURN(FAIL): fi: gu: end: #LoadedDiceEqs1(Lp1,Lp2,a,A,k):The non-homogeneous recurrence #equations for each of the variables i (i=1.. nops(Lp1)) #in the form [ope,rhs] for the integral of x times the prob. #dist induced by two loaded dice with prob. dist. #Lp1,Lp2 times x^k. For example, try: #LoadedDiceEqs1([1/2,1/2],[1/3,2/3],a,A,1); LoadedDiceEqs1:=proc(Lp1,Lp2,a,A,k) local i,gu,x: option remember: if nops(Lp1)<>nops(Lp2) or (type(convert(Lp1,`+`),numeric) and convert(Lp1,`+`)<>1 ) or (type(convert(Lp2,`+`),numeric) and convert(Lp2,`+`)<>1 ) then ERROR(`Bad input`): fi: gu:=x^k*LoadedDiceP(Lp1,Lp2,x,a): [seq(AZdL(gu,x,a[i],A[i],0,1),i=1..nops(Lp1))]: end: #IntDist(Lp1,Lp2,a0): the Integral of the distribution w.r.t. x #induced from two loaded dice whose faces are 1, 2, .., nops(Lp1) #with probability distibution Lp1 and Lp2 #For example, try: IntDist([1/2,1/2],[1/3,2/3],[2,2]); IntDist:=proc(Lp1,Lp2,a0) local mu,gu,x,a,A,a0S,r,lu,a0Sr,i,j,ope,rhs: option remember: if nops(Lp1)<>nops(Lp2) or (type(convert(Lp1,`+`),numeric) and convert(Lp1,`+`)<>1 ) or (type(convert(Lp2,`+`),numeric) and convert(Lp2,`+`)<>1 ) then ERROR(`Bad input`): fi: if nops(Lp1)<>nops(a0) then ERROR(`Bad input`): fi: mu:=LoadedDiceP(Lp1,Lp2,x,a): gu:=LoadedDiceEqs(Lp1,Lp2,a,A): if gu=FAIL then RETURN(FAIL): fi: for i from 1 to nops(a0) do if a0[i]>=degree(gu[i][1],A[i]) then a0S:=[op(1..i-1,a0),a0[i]-degree(gu[i][1],A[i]),op(i+1..nops(a0),a0)]: ope:=gu[i][1]: rhs:=gu[i][2]: lu:=subs({seq(a[j]=a0S[j],j=1..nops(a0S))},rhs): for r from 0 to degree(ope,A[i])-1 do a0Sr:=[op(1..i-1,a0),a0[i]-degree(gu[i][1],A[i])+r,op(i+1..nops(a0),a0)]: lu:= lu-subs({seq(a[j]=a0S[j],j=1..nops(a0S))},coeff(ope,A[i],r))* IntDist(Lp1,Lp2,a0Sr): od: lu:=lu/ subs({seq(a[j]=a0S[j],j=1..nops(a0S))},coeff(ope,A[i],degree(ope,A[i]))) : RETURN(lu): fi: od: mu:=subs({seq(a[j]=a0[j],j=1..nops(a0))},mu): int(mu,x=0..1): end: #IntDist1(Lp1,Lp2,a0,k): the Integral of the distribution w.r.t. x #induced from two loaded dice whose faces are 1, 2, .., nops(Lp1) #with probability distibution Lp1 and Lp2 (times x^k) #For example, try: IntDist1([1/2,1/2],[1/3,2/3],[2,2],1); IntDist1:=proc(Lp1,Lp2,a0,k) local mu,gu,x,a,A,a0S,r,lu,a0Sr,i,j,ope,rhs: option remember: if member(0,{op(a0)}) then mu:=x^k*LoadedDiceP(Lp1,Lp2,x,a): mu:=subs({seq(a[j]=a0[j],j=1..nops(a0))},mu): RETURN(int(mu,x=0..1)): fi: if nops(Lp1)<>nops(Lp2) or (type(convert(Lp1,`+`),numeric) and convert(Lp1,`+`)<>1 ) or (type(convert(Lp2,`+`),numeric) and convert(Lp2,`+`)<>1 ) then ERROR(`Bad input`): fi: if nops(Lp1)<>nops(a0) then ERROR(`Bad input`): fi: mu:=x^k*LoadedDiceP(Lp1,Lp2,x,a): gu:=LoadedDiceEqs1(Lp1,Lp2,a,A,k): if gu=FAIL then RETURN(FAIL): fi: for i from 1 to nops(a0) do if a0[i]>=degree(gu[i][1],A[i]) then a0S:=[op(1..i-1,a0),a0[i]-degree(gu[i][1],A[i]),op(i+1..nops(a0),a0)]: ope:=gu[i][1]: rhs:=gu[i][2]: lu:=subs({seq(a[j]=a0S[j],j=1..nops(a0S))},rhs): for r from 0 to degree(ope,A[i])-1 do a0Sr:=[op(1..i-1,a0),a0[i]-degree(gu[i][1],A[i])+r,op(i+1..nops(a0),a0)]: lu:= lu-subs({seq(a[j]=a0S[j],j=1..nops(a0S))},coeff(ope,A[i],r))* IntDist1(Lp1,Lp2,a0Sr,k): od: lu:=lu/ subs({seq(a[j]=a0S[j],j=1..nops(a0S))},coeff(ope,A[i],degree(ope,A[i]))) : RETURN(lu): fi: od: mu:=subs({seq(a[j]=a0[j],j=1..nops(a0))},mu): int(mu,x=0..1): end: #########End Bayesian Dice###################### ##Begin Prob. of random walk #pgf(L,x): Given a list of pairs [i,p] where the prob. of #of going from x to x+i is p, outputs the prob. generating #function w.r.t. x for one trial #For example, try: pgf( [[-1,1/2],[1,1/2]],x); pgf:=proc(L,x) local i: add(L[i][2]*x^L[i][1],i=1..nops(L)): end: #RecProbVisit(L,n,N,k): if you roll a die with faces and prob. dist. #given by a list L of pairs [i,p], n times, let a_k(n) be the #probability that you would have exactly k dollars after n rolls #outputs a linear recurrence equation with polynomial coefficients #for a_k(n) followed by initial conditions #For example, try: RecProbVisit([[1,1/2],[-1,1/2]],n,N,0): RecProbVisit:=proc(L,n,N,k) local ope,gu,x,lu,n1: lu:=pgf(L,x): gu:=lu^n/x^(k+1): ope:=AZd(gu,x,n,N)[1]: lu:=[seq(coeff(expand(lu^n1),x,k),n1=1..degree(ope,N))]: ope,lu: end: #RecProbVisitE(L,n,N,k): if you roll a die with faces and prob. dist. #given by a list L of pairs [i,p], 2n times, let a_k(n) be the #probability that you would have exactly k dollars after 2n rolls #outputs a linear recurrence equation with polynomial coefficients #for a_k(n) followed by initial conditions #For example, try: RecProbVisitE([[1,1/2],[-1,1/2]],n,N,0): RecProbVisitE:=proc(L,n,N,k) local ope,gu,x,lu,n1: lu:=pgf(L,x): gu:=lu^(2*n)/x^(k+1): ope:=AZd(gu,x,n,N)[1]: lu:=[seq(coeff(expand(lu^(2*n1)),x,k),n1=1..degree(ope,N))]: ope,lu: end: ####Begin Asymptotics with(SolveTools): Yafe:=proc(ope,N) local i,ope1,coe1,L: if ope=0 then RETURN(1,0): fi: ope1:=expand(ope): L:=degree(ope1,N): coe1:=coeff(ope1,N,L): ope1:=normal(ope1/coe1): ope1:=normal(ope1): ope1:= convert( [seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`): factor(coe1),ope1: end: #Aope1(ope,n,N): the asymptotics of the difference operator with poly. coeffs. ope(n,N) Aope1:=proc(ope,n,N) local gu: gu:=expand(numer(normal(ope))): gu:=coeff(gu,n,degree(gu,n)): gu:=expand(gu/coeff(gu,N,degree(gu,N))): factor(gu): end: C3:=proc(a,b,c) option remember: if [a,b,c]=[0,0,0] then RETURN(1): fi: if a<0 or b<0 or c<0 then RETURN(0): fi: C3(a-1,b,c)+C3(a,b-1,c)+C3(a,b,c-1)-4*C3(a-1,b-1,c-1): end: #Aluf1(pit): is it positive dominant? Aluf1:=proc(pit) local aluf,shi,i,makom: shi:=evalf(abs(pit[1])): aluf:={evalf(evalc(pit[1]))}: makom:={1}: for i from 2 to nops(pit) do if evalf(abs(pit[i]))>evalf(shi) then shi:=abs(evalf(evalc(pit[i]))): aluf:={pit[i]}: makom:={i}: elif evalf(abs(pit[i]))=shi then aluf:=aluf union {evalf(evalc(pit[i]))}: makom:=makom union {i}: fi: od: aluf,shi,makom: end: #OneStepAS1(ope1,n,N,alpha,f,S1): Given the partial asymptotic #expansion of solutions of ope1(n,N) with exponent alpha #extends it by one term OneStepAS1:=proc(ope1,n,N,alpha,f,S1) local x,f1,L,F,A,mu,i,A1: f1:=subs(n=1/x,f): L:=degree(f1,x): F:=f1+A*x^(L+S1): mu:=add( subs(n=1/x,coeff(ope1,N,i))*(1+i*x)^alpha*subs(x=x/(1+i*x),F) ,i=ldegree(ope1,N)..degree(ope1,N)): mu:=normal(mu): mu:=taylor(mu,x=0,L+11): mu:=simplify(mu): for i from L+1 to L+9 while coeff(mu,x,i)=0 do od: if i=L+10 then RETURN(FAIL): fi: mu:=coeff(mu,x,i): A1:=subs(Linear({mu},{A}),A): if A1=0 then RETURN(FAIL): fi: subs({x=1/n,A=A1},F): end: #OneStepA(ope1,n,N,alpha,f): Given the partial asymptotic #expansion of solutions of ope1(n,N) with exponent alpha #extends it by one term OneStepA:=proc(ope1,n,N,alpha,f) local S1: for S1 from 1 while OneStepAS1(ope1,n,N,alpha,f,S1)=FAIL do od: OneStepAS1(ope1,n,N,alpha,f,S1): end: #Asy(ope,n,N,K): the asymptotic expansion of solutions #to ope(n,N)f(n)=0,where ope(n,N) is a recurrence operator #up to the K's term Asy:=proc(ope,n,N,K) local gu,pit,lu,makom,alpha,mu,ope1,ku,i,f,x,ka: gu:=Aope1(ope,n,N): if degree(gu,N)=0 then print(`Not of exponential growth, case not implemented`): RETURN(FAIL): fi: if not type(expand(normal(ope/gu)),`+`) then pit:=[solve(gu,N)]: makom:=Aluf1(pit)[3]: RETURN(pit[makom[1]]^n*expand((nops(makom)-1)!*binomial(nops(makom)-1+n,nops(makom)-1))): fi: pit:=[solve(gu,N)]: makom:=Aluf1(pit)[3]: if nops(makom)<>1 then lu:=evalf(evalc(pit[makom[1]])): if coeff(lu,I,1)>(1/10)^(Digits-2) then print(`Not only that but Dominant roots are complex`): RETURN(FAIL): fi: fi: lu:=evalf(evalc(pit[makom[1]])): if coeff(lu,I,1)>(1/10)^(Digits-2) then print(`Dominant root is complex`): RETURN(FAIL): fi: lu:=pit[makom[1]]: ope1:=Yafe(subs(N=lu*N,ope),N)[2]: mu:=add( subs(n=1/x,coeff(ope1,N,i))*(1+i*x)^alpha,i=ldegree(ope1,N)..degree(ope1,N)): mu:=normal(mu): ku:=factor(coeff(taylor(mu,x=0,nops(makom)+1),x,nops(makom))): ka:=simplify([solve(ku,alpha)]): alpha:=max(op(ka)): if normal(subs(N=1,ope1))=0 then RETURN(lu^n*n^alpha): fi: f:=1: for i from 1 to K do f:=OneStepA(ope1,n,N,alpha,f): od: lu^n*n^alpha*f: end: ###EndAsymptotics #SeqFromRec(ope,n,N,Ini,K): Given the first L-1 #terms of the sequence Ini=[f(1), ..., f(L-1)] #satisfied by the recurrence ope(n,N)f(n)=0 #extends it to the first K values SeqFromRec:=proc(ope,n,N,Ini,K) local ope1,gu,L,n1,j1: ope1:=Yafe(ope,N)[2]: L:=degree(ope1,N): if nops(Ini)<>L then ERROR(`Ini should be of length`, L): fi: ope1:=expand(subs(n=n-L,ope1)/N^L): gu:=Ini: for n1 from nops(Ini)+1 to K do gu:=[op(gu), -add(gu[nops(gu)+1-j1]*subs(n=n1,coeff(ope1,N,-j1)), j1=1..L)]: od: gu: end: #AsyC(ope,n,N,K,Ini,L): the asymptotic expansion #complete with the estimated constant, of solutions #to ope(n,N)f(n)=0, with initial conditions #Ini, where ope(n,N) is a recurrence operator #up to the K's term. The constant is estimated using up to L terms #For example, try #AsyC(N^2-N-1,n,N,1,[1,1],100); AsyC:=proc(ope,n,N,K,Ini,L) local lu,gu,C,dig,i: Digits:=100: if nops(Ini)