####################################################################### ## qZEILBERGER Save this file as qZEILBERGER. To use it, stay in # ## the same directory, get into Maple (by typing: maple ) # ## and then type: read qZEILBERGER : # ## Then follow the instructions given there # ## # ## Written by Mohamud Mohammed and Doron Zeilberger, # ## Rutgers University , # ## mohamudm at math dot rutgers dot edu. # ## zeilberg at math dot rutgers dot edu. # ####################################################################### print(`Written by Mohamud Mohammed and Doron Zeilberger.`): print(`First Written: July 22, 2004; tested for Maple 8 and 9. `): print(`Version of Aug. 16, 2004 (Thanks to Drew Sills!). `): print(): print(`This is qZEILBERGER, A Maple package`): print(`that implements the new-improved (and simplified) q-Zeilberger`): print(`algorithm for single-sum definite hypergeometric summation.`): print(`It does what qzeil does in qEKHAD, and what is done in `): print(`In Axel Riese's Mathematica package qZeil. `): print(`This is one of the two packages`): print(`that accompany the article: Sharp Upper Bounds for the `): print(`Orders of the Recurrences Outputted by the Zeilberger`): print(` and q-Zeilberger algorithms`): print(`by Mohamud Mohammed and Doron Zeilberger`): print(`available from the authors' websites `): 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 .`): print(): print(`For general help, and a list of the available functions,`): print(` type "ezra();". For specific help type "ezra(procedure_name)" . `): print(): print(`Warning: q is a global variable `): ezra:=proc() if args=NULL then print(` qZEILBERGER: `): print(` A Maple package for `): print(` finding recurrences satisfied by definite q-hypergeomeric sums . `): print(` For help with a specific procedure, type ezra(procedure_name); . `): print(`The main procedure is: qDoron . `): print(): fi: if nargs=1 and args[1]=qDoron then print(` qDoron(F,k,n,N,POL,qQuad): Given a pure`): print(`(i.e without a polynomial part), q-hypergeometric term F in the `): print(`variables k and n, and a symbol N for the shift operator in n:`): print(` [i.e. Nf(n):=f(n+1)]`): print(`POL,The polynomial part (a polynomial in q,q^k,q^n) and `): print(`qQuad, a quadratic polynomial in n and k`): print(`Let SUMMAND(n,k):=F*POL*q^qQuad`): print(` returns an operator ope(n,N) and`): print(` the certificate, R(n,k), a rational function in q^k,q^n `): print(` such that G(n,k):=R(n,k)*SUMMAND(n,k) satisfies `): print(` ope(N,n)SUMMMAND(k,n)=G(n,k+1)-G(n,k)`): print(``): print(` A pure q-hypergeometric term is a product of factorials of the form`): print(`[an+bk+c]! where a is a NON-negative integer, b is an integer,`): print(`and c is symbolic, times a polynomial in (n,k) times x^k for some`): print(`numeric or symbolic x`): print(`One can also use binomial(n,k) `): print(``): print(`WARNING: In this package a! means (q)_a:=(1-q)(1-q^2)...(1-q^a)`): print(`and binomial(n,k) means the q-binomial: (q)_n/((q)_k*(q)_{n-k})`): print(``): print(` For example, try:`): print(` qDoron(1/k!^2/(n-k)!,k,n,N,1,k*(k-1)/2);`): print(` qDoron(binomial(n,k)*x^k,k,n,N,1,k*(k-1)/2);`): print(` qDoron(binomial(n,k)^2,k,n,N,1,k^2);`): print(` qDoron((-1)^k*binomial(2*n,n+k)^3,k,n,N,1,k*(3*k-1)/2);`): print(`This is q-Dixon w/o Paule symmetrization`): print(` qDoron((-1)^k*binomial(2*n,n+k)^3,k,n,N,1+q^k,k*(3*k-1)/2);`): print(`This is q-Dixon with Paule symmetrization`): fi: if nargs=1 and args[1]=PqDoron1 then print(` PqDoron1(F,k,n,N,POL,qQuad): Given a pure`): print(`(i.e without a polynomial part), q-hypergeometric term F in the `): print(`variables k and n, and a symbol N for the shift operator in n:`): print(` [i.e. Nf(n):=f(n+1)]`): print(`POL,The polynomial part (a polynomial in q,q^k,q^n) and `): print(`qQuad, a quadratic polynomial in n and k`): print(`Let SUMMAND(n,k):=F*POL*q^qQuad`): print(` returns an operator ope(n,N) and`): print(` a pre-certificate, a polynomial from which one can construct`): print(` G(n,k), such that `): print(` ope(N,n)SUMMMAND(k,n)=G(n,k+1)-G(n,k)`): print(`Note: to get the certificate itself use qDoron)`): print(``): print(` A pure q-hypergeometric term is a product of factorials of the form`): print(`[an+bk+c]! where a is a NON-negative integer, b is an integer,`): print(`and c is symbolic, times a polynomial in (n,k) times x^k for some`): print(`numeric or symbolic x`): print(`One can also use binomial(n,k) `): print(``): print(`WARNING: In this package a! means (1-q)(1-q^2)...(1-q^a)`): print(`and binomial(n,k) means the q-binomial`): print(``): print(` For example, try:`): print(` PqDoron1(binomial(n,k)*x^k,k,n,N,1,k*(k-1)/2);`): print(` PqDoron1(binomial(n,k)^2,k,n,N,1,k^2);`): print(` PqDoron1((-1)^k*binomial(2*n,n+k)^3,k,n,N,1+q^k,k*(3*k-1)/2);`): fi: end: ez:=proc(): print(`c(Ali,Bli,Cli,Dli,POL,e,n,L)`): print(` a(Num,Den,z,k,n,L) , b(Num,Den,k,n,L) :`): print(`Equ1(Num,Den,POL,z,k,n,L,e,X)`): print(`Z1a(Ali,Bli,Cli,Dli,POL,z,k,n,N)`): print(`Z1(F,k,n,N,POL,qQuad),(F=[Num,Den,Pol,z]) `): print(`Z1L(F,k,n,L,N,POL,qQuad),(F=[Num,Den,Pol,z]) `): print(`Z1aL(Ali,Bli,Cli,Dli,POL,z,k,n,L,N)`): print(`hafokh(F,k,n)`): print(`PqDoron1(F,k,n,N,POL,qQuad)`): print(`RecEq(F,k,n,N,low1,high1)`): print(`qkTot(Pol,k,t)`): end: rf:=proc(a,n) local i: if n<0 then 0 elif n=0 then 1 else convert([seq(1-q^(a+i),i=0..n-1)],`*`) fi: end: RatioH:=proc(Num,Den,n,L,i,qQuad) local j,gu: gu:= convert([seq( rf(Num[j]+1,i*coeff(Num[j],n,1)),j=1..nops(Num))],`*`)* convert([seq( rf( subs(n=n+i,Den[j])+1,(L-i)*coeff(Den[j],n,1)),j=1..nops(Den))],`*`): gu:=gu*expand(q^(expand(subs(n=n+i,qQuad)-qQuad))): end: c:=proc(Num,Den,POL,e,n,L,qQuad) local i:convert( [seq(e[i]*subs(n=n+i,POL)*RatioH(Num,Den,n,L,i,qQuad),i=0..L)],`+`): end: a:=proc(Num,Den,z,k,n,L,qQuad) local j,gu: gu:=z: for j from 1 to nops(Num) do if coeff(Num[j],k,1) > 0 then gu:=gu*rf(Num[j]+1,coeff(Num[j],k,1)): fi: od: for j from 1 to nops(Den) do if coeff(Den[j],k,1) < 0 then gu:=gu*rf(subs({n=n+L,k=k+1},Den[j])+1,-coeff(Den[j],k,1)): fi: od: gu:=gu*q^(expand(subs(k=k+1,qQuad)-qQuad)): gu: end: b:=proc(Num,Den,k,n,L) local j,gu: gu:=1: for j from 1 to nops(Num) do if coeff(Num[j],k,1) < 0 then gu:=gu*rf(subs(k=k+1,Num[j])+1,-coeff(Num[j],k,1)): fi: od: for j from 1 to nops(Den) do if coeff(Den[j],k,1) > 0 then gu:=gu*rf(subs({n=n+L},Den[j])+1,coeff(Den[j],k,1)): fi: od: gu: end: Equ1:=proc(Num,Den,POL,z,k,n,L,e,X,qQuad): a(Num,Den,z,k,n,L,qQuad)*X(k+1)- subs(k=k-1,b(Num,Den,k,n,L))*X(k)- c(Num,Den,POL,e,n,L,qQuad): end: yafe:=proc(ope,N) local i,ope1:ope1:=expand(ope): convert( [seq(factor(coeff(ope1,N,i))*N^i,i=ldegree(ope1,N)..degree(ope1,N))],`+`): end: FindDeg:=proc(a1,b1,c1,k) local degX,lam,A,B,d: degX:=degree(c1,k)-max(degree(a1,k),degree(b1,k)): if degree(a1,k)=degree(b1,k) and coeff(expand(a1),k,degree(a1,k))=coeff(expand(b1),k,degree(b1,k)) then lam:=coeff(expand(b1),k,degree(b1,k)): B:=coeff(expand(b1),k,degree(b1,k)-1): A:=coeff(expand(a1),k,degree(a1,k)-1): d:=normal((B-A)/lam): if type(d,integer) then degX:=max(degX+1,d): else degX:=degX+1: fi: fi: degX: end: Z1L:=proc(F,k,n,L,N,POL,qQuad) local e,x,j,degX,X,a1,b1,c1,gu,var1, var,ope,lu,mu,cert,eq,eqTry,varTry,Num,Den,z,t,a1t,b1t,c1t,Y,Y1,Z,Eq1, LowD,HighD,ope1,cert1: Num:=F[1]: Den:=F[2]: z:=F[3]: ope:=convert([seq(e[j]*N^j,j=0..L)],`+`): a1:=a(Num,Den,z,k,n,L,qQuad): b1:=subs(k=k-1,b(Num,Den,k,n,L)): c1:=c(Num,Den,POL,e,n,L,qQuad): a1t:=qkTot(a1,k,t): b1t:=qkTot(b1,k,t): c1t:=qkTot(c1,k,t): Eq1:=numer(normal(a1t*Y1-b1t*Y-c1t*Z)): a1t:=coeff(Eq1,Y1): b1t:=-coeff(Eq1,Y): c1t:=-coeff(Eq1,Z): LowD:= min(ldegree(coeff(Eq1,Z,1),t)-ldegree(coeff(Eq1,Y,1),t), ldegree(coeff(Eq1,Z,1),t)-ldegree(coeff(Eq1,Y1,1),t))-1: HighD:= min(degree(coeff(Eq1,Z,1),t)-degree(coeff(Eq1,Y,1),t), degree(coeff(Eq1,Z,1),t)-degree(coeff(Eq1,Y1,1),t))+1: X:=convert([seq(x[j]*t^j,j=LowD..HighD)],`+`): gu:=expand(a1t*subs(t=t*q,X)-b1t*X-c1t): var1:={seq(x[j],j=LowD..HighD),seq(e[j],j=0..L)}: eq:={seq(coeff(gu,t,j),j=ldegree(gu,t)..degree(gu,t))}: if type(z,symbol) then eqTry:=subs({z=2, n=5/3},eq): else eqTry:=subs({q=2,n=3},eq): fi: varTry:=solve(eqTry,var1): if subs(varTry,ope)=0 then RETURN(FAIL): fi: if type(z,symbol) then eqTry:=subs({z=3,n=7,q=2},eq): else eqTry:=subs({q=2,n=7},eq): fi: if subs(varTry,ope)=0 then RETURN(FAIL): fi: print(`Most likely there is an operator of order`, L, `please be patinet.`): var:=solve(eq,var1): mu:=subs(var,ope),subs(var,X): if mu[1]=0 then RETURN(FAIL): fi: if mu[2]=0 then RETURN(mu): fi: lu:=normal(mu[1]/mu[2]): ope:=yafe(numer(lu),N): cert:=factor(denom(lu)): ope1:=yafe(expand(normal(ope/t^degree(ope,t))),N): cert1:=normal(cert/t^degree(ope,t)): ope1,subs(t=q^k,cert1): end: Z1:=proc(F,k,n,N,POL,qQuad) local L,gu: gu:=Z1L(F,k,n,0,N,POL,qQuad): for L from 1 while gu=FAIL do gu:=Z1L(F,k,n,L,N,POL,qQuad): od: gu: end: ####converting to the input needed by Z1### FindFactorial:=proc(Prod1) local Prod2,i,lu,lu1,lu2: Prod2:=factor(Prod1): if type(Prod2,function) and op(0,Prod2)=factorial then RETURN(Prod1): fi: for i from 1 to nops(Prod2) do lu:=op(i,Prod2): if type(lu,function) and op(0,lu)=factorial then RETURN(lu): fi: if type(lu,`^`) then lu1:=op(1,lu): lu2:=op(2,lu): if type(lu1,function) and op(0,lu1)=factorial then if not type(lu2,integer) and lu2>0 then RETURN(FAIL): else RETURN(lu1): fi: fi: fi: od: FAIL: end: FindPower:=proc(Prod1,k) local Prod2,i,lu: Prod2:=factor(Prod1): if Prod2=1 then RETURN(FAIL): fi: if type(Prod2,`^`) and degree(op(2,Prod1),k)=1 then RETURN(Prod1): fi: for i from 1 to nops(Prod2) do lu:=op(i,Prod2): if type(lu,`^`) and degree(op(2,lu),k)=1 then RETURN(lu): fi: od: FAIL: end: hafokh:=proc(F,k,n) local F1,Num,Den,POL,z,F1t,F1b,lu,lu1,lu2,i: F1:=factor(normal(F)): F1t:=factor(numer(F1)): F1b:=factor(denom(F1)): Num:=[]: lu:=FindFactorial(F1t): while lu<>FAIL do lu1:=op(lu): if degree(lu1,{n,k})>1 or not type(coeff(lu1,k,1),integer) or not type(coeff(lu1,n,1),integer) then RETURN(FAIL): fi: Num:=[op(Num),lu1]: F1t:=normal(F1t/lu): lu:=FindFactorial(F1t): od: z:=1: lu:=FindPower(F1t,k): while lu<>FAIL do lu1:=op(1,lu): lu2:=op(2,lu): if degree(lu2,k)>1 then RETURN(FAIL): else z:=z*lu1^coeff(lu2,k,1): F1t:=normal(simplify(normal(F1t/lu1^(k*coeff(lu2,k,1))))): fi: lu:=FindPower(F1t,k): od: POL:=F1t: Den:=[]: lu:=FindFactorial(F1b): while lu<>FAIL do lu1:=op(lu): if degree(lu1,{n,k})>1 or not type(coeff(lu1,k,1),integer) or not type(coeff(lu1,n,1),integer) then RETURN(FAIL): fi: Den:=[op(Den),lu1]: F1b:=normal(F1b/lu): lu:=FindFactorial(F1b): od: lu:=FindPower(F1b,k): while lu<>FAIL do lu1:=op(1,lu): lu2:=op(2,lu): if degree(lu2,k)>1 then RETURN(FAIL): else z:=z/lu1^coeff(lu2,k,1): F1b:=normal(simplify(normal(F1b/lu1^(k*coeff(lu2,k,1))))): fi: lu:=FindPower(F1b,k): od: POL:=normal(POL/F1b): for i from 1 to nops(Num) do if coeff(Num[i],n,1)<0 then RETURN(FAIL): fi: od: for i from 1 to nops(Den) do if coeff(Den[i],n,1)<0 then RETURN(FAIL): fi: od: if POL<>1 then RETURN(FAIL): fi: [Num,Den,z]: end: ##End conversion part PqDoron1:=proc(F,k,n,N,POL,qQuad) local F1: F1:=normal(convert(F,factorial)): F1:=hafokh(F1,k,n): if F1=FAIL then print(`Bad input`): RETURN(FAIL): fi: Z1(F1,k,n,N,POL,qQuad): end: qDoron:=proc(F,k,n,N,POL,qQuad) local lu, F1,b1,kadima,L: F1:=normal(convert(F,factorial)): F1:=hafokh(F1,k,n): lu:=PqDoron1(F,k,n,N,POL,qQuad): L:=degree(lu[1],N): b1:=b(F1[1],F1[2],k,n,L): kadima:=RatioH(F1[1],F1[2],n,L,0,qQuad): combine( [lu[1],normal(lu[2]/kadima*subs(k=k-1,b1)/F1[3])] ,power,symbolic): end: #RecEq(F,k,n,N,low1,high1): a recurrence equation for #sum(F(n,k),k=low1..high1) where low1 and high1 have the form #a+b*n, with b a non0negative integer RecEq:=proc(F1,k,n,N,low1,high1) local lu,ope,cert,i,i1,gu,F: F:=convert(F1,factorial): lu:=qDoron(F,k,n,N): ope:=lu[1]: cert:=lu[2]: gu:=0: for i from 0 to degree(ope,N) do gu:=coeff(ope,N,i)* (-convert([seq(subs({n=n+i,k=low1+i1},F),i1=0..coeff(low1,n,1)*i-1)],`+`) +convert([seq(subs({n=n+i,k=high1+i1},F),i1=1..coeff(high1,n,1)*i)],`+`)): od: gu:=gu+ subs(k=high1+1,cert)*subs(k=high1+1,F)- subs(k=low1,cert)*subs(k=low1,F): gu:=normal(simplify(normal(gu))): ope,gu: end: #qkTot(Pol,k,t): Given a polynomial Pol in q^k subsitutes q^k by t qkTot:=proc(Pol,k,t) local P1: P1:=expand(Pol): subs(q^k=t,P1): end: