# file zeilb for summation of terminating hypergeometric series # by Zeilberger's algorithm. # # Copyright 1992, 1996, 1998, 2000 by Tom H. Koornwinder (thk@science.uva.nl) # # The 1992 version was a completely rewritten version of a Maple procedure # ident_prover.maple by D. Zeilberger. # # The 1996 version was adapted to Maple V, Release 3 and 4. # # The 1998 version was adapted to Maple V, Release 4 and 5. # # The present version, 9 December 2000, works in Maple 6, as well as in # Maple V, Release 4 and 5. fac:=proc(a,n) local k: option `Copyright 1992, 1996, 1998 by Tom H. Koornwinder`: if n=0 then 1 elif type(n,posint) then product('a+k','k'=0..n-1) elif type(n,negint) then product('(a-k)^(-1)','k'=1..-n) else RETURN('procname(args)') fi end: comppol:=proc(f1,g1,x) local f,g,n,a,b,c,d,j: option `Copyright 1992, 1996, 1998 by Tom H. Koornwinder`: # Tests for polynomials f1 and g1 in x if f1(x)=g1(x+j) for some nonnegative # integer j, while f1 and g1 are not constant in x f:=collect(f1,x): g:=collect(g1,x): n:=degree(f,x): if n=0 or n<>degree(g,x) then RETURN(-1): fi: a:=coeff(f,x,n): b:=coeff(f,x,n-1): c:=coeff(g,x,n): d:=coeff(g,x,n-1): j:=normal((b*c-a*d)/n/a/c): if not type(j,nonnegint) then RETURN(-1): fi: if collect(c*f-a*subs(x=x+j,g),x)=0 then RETURN(j): else RETURN(-1): fi: end: maxshift:=proc(r1,r2,x) local fr1,fr2,i,j,m,ma,n: option `Copyright 1992, 1996, 1998 by Tom H. Koornwinder`: ma:=-1: fr1:=factors(r1): fr2:=factors(r2): m:=nops(op(2,fr1)): n:=nops(op(2,fr2)): for i from 1 to m do for j from 1 to n do ma:=max(ma,comppol(op(1,op(i,op(2,fr1))),op(1,op(j,op(2,fr2))),x)): od: od: RETURN(ma): end: testnonnegintroots:=proc(p,k) local f,fp,i,k0: option `Copyright 1992, 1996, 1998 by Tom H. Koornwinder`: k0:=-1: fp:=factors(p): for i from 1 to nops(op(2,fp)) do f:=collect(op(1,op(i,op(2,fp))),k): if degree(f,k)=1 and coeff(f,k)=1 and type(-subs(k=0,f),nonnegint) then k0:=max(k0,-subs(k=0,f)) fi: od: RETURN(k0): end: refactors:=proc(f,n) local A,i,c,d,g,pol,deg: option `Copyright 1992, 1996, 1998 by Tom H. Koornwinder`: A:=op(1,f): g:=1: for i from 1 to nops(op(2,f)) do pol:=op(1,op(i,op(2,f))): deg:=op(2,op(i,op(2,f))): if degree(pol,n)<2 then d:=coeff(pol,n,1): c:=coeff(pol,n,0)+d: if degree(pol,n)=1 then A:=A*d^deg: c:=c/d: g:=g*fac(c,n)^deg: else A:=A*c^deg: fi: else RETURN(FAIL): fi: od: A^n*g: end: zeilb:=proc() local A, B, BO, BOdeg, C, D, EE, F, F1, INH, ORDER, P, P0, P1, R1, R2, R3, R4, S, TO, TOdeg, a, b, co, co1, co2, deg, deg1, deg2, eq, g, i, inh, j, k, m, ma, n, n0, nonunique, p, r, ra1, ra2, s, sol, t, talklevel, time0, va, vb, vbn, var, z; options `Copyright 1992, 1996, 1998, 2000 by Tom H. Koornwinder`; time0:=time(): if nargs<5 then ERROR(`too few parameters in function call`): fi: if not type(args[4],function) or nops(args[4])<>1 then ERROR(`fourth argument is not a function of one variable`) fi: n:=op(1,args[4]): if not type(n,name) then ERROR(`fourth argument is not of form f(n) with n being a name`): fi: ORDER:=args[5]: if not type(ORDER,nonnegint) then ERROR(`fifth argument is not a nonnegative integer`): fi: n0:=ORDER-1: if nargs>5 then talklevel:=args[6]: else talklevel:=max(1,printlevel): fi: if not type(talklevel,posint) then ERROR(`sixth argument is not a positive integer`): fi: TO:=args[1]: if not type(TO,list(ratpoly)) then ERROR(`first argument is not a list of rational functions`): fi: BO:=args[2]: if not type(BO,list(ratpoly)) then ERROR(`second argument is not a list of rational functions`): fi: z:=args[3]: if not type(z,ratpoly) then ERROR(`third argument is not a rational function`): fi: TO:=map(normal,TO,expanded): BO:=map(normal,BO,expanded): z:=normal(z): if z=0 then RETURN(1): fi: if has(z,n) then ERROR(cat(`third argument depends on `,n)): fi: for i while i <= nops(TO) do if member(TO[i], BO, 'j') then TO := subsop(i = NULL, TO): BO := subsop(j = NULL, BO): i := i-1: fi: od: r:=nops(TO): s:=nops(BO): if has(map(type,TO,integer),true) then ERROR(`some upper index is integer`): fi: if has(BO,0) or has(map(type,BO,negint),true) then ERROR(`some lower index is nonpositive integer`): fi: if ORDER=0 then if has(TO,n) or has(BO,n) then ERROR(cat(`order=0 while some index depends on `,n)): fi: fi: if ORDER>0 then if not member(-n,TO) then ERROR(cat(`order>0 and no upper index is equal to `,`-n`)): fi: for i from 1 to r do if has(TO[i],n) and not type(TO[i],linear(n)) then ERROR(cat(`some upper index is not polynomial of degree <=1 in `,n)): fi: if not type(coeff(TO[i],n),integer) then ERROR(cat(`some upper index has non-integer coefficient of `,n)): fi: od: for i from 1 to s do if has(BO[i],n) and not type(BO[i],linear(n)) then ERROR(cat(`some lower index is not polynomial of degree <=1 in `,n)): fi: if not type(coeff(BO[i],n),integer) then ERROR(cat(`some lower index has non-integer coefficient of `,n)): fi: if type(coeff(BO[i],n),negint) and type(subs(n=0,BO[i]),integer) then ERROR(cat(`some lower index equals integer minus negative integer \ times `,n)): fi: od: fi: if talklevel>2 then print(`zeilb: computing P, R1 and R2`): fi: if talklevel>3 then print(cat(`time=`,convert(time()-time0,string))): fi: TOdeg:=map(coeff,TO,n): BOdeg:=map(coeff,BO,n): B:=product('fac(TO[t]-TOdeg[t]+k,TOdeg[t])/ fac(TO[t]-TOdeg[t],TOdeg[t])','t'=1..r)* product('fac(BO[t]-BOdeg[t],BOdeg[t])/ fac(BO[t]-BOdeg[t]+k,BOdeg[t])','t'=1..s): B:=normal(B): C:=denom(B): B:=numer(B): D:=product('TO[t]+k-1','t'=1..r)/ product('BO[t]+k-1','t'=1..s)/k*z: D:=normal(D): EE:=denom(D): D:=numer(D): P0:=product('subs(n=n-i,B)','i'=0..ORDER-1): for j from 1 to ORDER do P0:=P0+b[j]*product('subs(n=n-i,C)','i'=0..j-1)* product('subs(n=n-i,B)','i'=j..ORDER-1): od: R1:=D*product('subs(n=n-i,k=k-1,B)','i'=0..ORDER-1)/ EE/product('subs(n=n-i,B)','i'=0..ORDER-1): R1:=normal(R1): R2:=denom(R1): R1:=numer(R1): ma:=maxshift(R1,R2,k): j:=0: P1:=1: while j <= ma do g:=gcd(R1,subs(k=k+j,R2)): if g <> 1 then R1:=normal(R1/g): R2:=normal(R2/subs(k=k-j,g)): P1:=P1*product('subs(k=k-p,g)' , 'p'=0..j-1): fi: j:=j+1: od: P:=expand(P0*P1,k): if talklevel>2 then print(`zeilb: computing maximal degree of P`): fi: if talklevel>3 then print(cat(`time=`,convert(time()-time0,string))): fi: R3:=expand(subs(k=k+1,R1)+R2,k): deg1:=degree(R3,k): co1:=coeff(R3,k,deg1): R4:=expand(subs(k=k+1,R1)-R2,k): if R4=0 then deg2:=-1: else deg2:=degree(R4,k): fi: co2:=coeff(R4,k,deg1-1): if deg1 <= deg2 then deg:=degree(P,k)-deg2 else co:= -2*co2/co1: co:=simplify(co): if type(co,nonnegint) then deg:=max(co,degree(P,k)-deg1+1) else deg:=degree(P,k)-deg1+1: fi: fi: if deg < 0 then print(`Algorithm fails`): if talklevel>1 then print(`degree(F)<0`): fi: RETURN(): fi: if talklevel>2 then print(`zeilb: solving equations`): fi: if talklevel>3 then print(cat(`time=`,convert(time()-time0,string))): fi: F:=sum('a[i]*k^i','i'=0..deg): F1:=expand(subs(k=k+1,R1)*F-R2*subs(k=k-1,F)-P): deg1:=degree(F1,k): eq:=seq(coeff(F1,k,i),i=0..deg1): va := seq(a[i],i=0..deg): vb := seq(b[i],i=1..ORDER): var:=vb,va: sol:=solve({eq},{var}): # eq := Groebner['gbasis']([eq], 'plex'(var) ); # eq := convert(eq, 'set'): # # select equations involving b[i]'s # eqb := select(has,eq,{vb}): # eqa := eq minus eqb: # equation with a[i]'s only # sola := solve(eqa,{va}): # solb := solve(eqb,{vb}): # sol:=sola union solb: if sol=NULL then print(`Algorithm fails`): if talklevel>1 then print(`System of equations has no solutions`): fi: RETURN(): fi: if ORDER=0 then F:=subs(sol,F): if testnonnegintroots(numer(normal(F)),k)>-1 then print(`Algorithm fails`): if talklevel>1 then print(`F(k)=0 at some nonnegative integer`): fi: RETURN(): fi: A:=product('fac(TO[j],k)','j'=1..r)/ product('fac(BO[j],k)','j'=1..s)/fac(1,k)*z^k: RETURN(factor(subs(k=n+1,R1)*subs(k=n,A)*subs(k=n,F)/subs(k=n,P))): fi: nonunique:=has({seq(op(2,sol[i]),i=1..deg+ORDER+1)},{var}); va:=subs(sol,[va]): vb:=subs(sol,[vb]): if normal(1+sum('vb[i]','i'=1..ORDER))=0 then print(`Algorithm fails`): if talklevel>1 then print(`a(0) is equal to 0`): fi: RETURN(): fi: P:=subs(sol,P): if testnonnegintroots(numer(normal(P)),k)>-1 then print(`Algorithm fails`): if talklevel>1 then print(`P(k)=0 at some nonnegative integer`): fi: RETURN(): fi: F:=subs(sol,F): if testnonnegintroots(numer(normal(F)),k)>-1 then print(`Algorithm fails`): if talklevel>1 then print(`F(k)=0 at some nonnegative integer`): fi: RETURN(): fi: for j from 1 to nops(va) do n0:=max(testnonnegintroots(denom(normal(va[j])),n),n0): od: for j from 1 to nops(vb) do n0:=max(testnonnegintroots(denom(normal(vb[j])),n),n0): od: n0:=max(testnonnegintroots(denom(normal(subs(k=n+1,R2))),n),n0): n0:=max(testnonnegintroots(numer(normal(subs(k=n+1,P))),n),n0): if ORDER =1 and not nonunique then vbn:=-normal(vb[1]): ra1:=factors(numer(vbn)): ra2:=factors(denom(vbn)): ra1:=refactors(ra1,n): ra2:=refactors(ra2,n) fi: INH:=op(0,args[4]): inh:='INH': if talklevel>2 then S:=Sum(inh(n,k)+sum('vb[j]*inh(n-j,k)','j'=1..ORDER),'k'=0..m): print(`Write the input series as`): print(INH(n)=Sum(inh(n,k),'k'=0..n)): print(`Then `): print(S=subs(k=m+1,R2)*subs(k=m,F)/subs(k=m+1,P1)/ product('subs(n=n-i,k=m+1,B)','i'=0..ORDER-1)*inh(n,m+1),n>n0): if talklevel>3 then print(cat(`time =`,convert(time()-time0,string))): fi: print(cat(`Hence `,INH,`(`,n,`) equals`)): fi: if talklevel=2 then print(cat(`The input series `,INH,`(`,n,`) equals`)): fi: if ORDER>1 or nonunique or n0>=ORDER or ra1=FAIL or ra2=FAIL then RETURN(-sum('factor(vb[i])*INH(n-i)','i'=1..ORDER),n>n0): fi: if talklevel>1 then print(-sum('factor(vb[i])*INH(n-i)','i'=1..ORDER),n>n0): print(cat(`Hence `,INH,`(`,n,`) equals`)): fi: RETURN(simplify(ra1/ra2)): end: