####################################################################### ## qMultiZeilberger: Save this file as qMultiZeilberger. To use # ## it, stay in the # ## same directory, get into Maple (by typing: maple ) # ## and then type: read qMultiZeilberger : # ## Then follow the instructions given there # ## # ## Written by Mohamud Mohammed and Doron Zeilberger, # ## Rutgers University , # ## [mohamudm, zeilberg] at math.rutgers.edu. # ####################################################################### print(`First Written: July 27,2004: tested for Maple 8 `): print(`Version of July 27, 2004: `): print(): print(`This is qMultiZeilberger, A Maple package`): print(`accompanying the article `): print(`Multi-Variable Zeilberger and Almkvist-Zeilberger Algorithms`): print(`and the Sharpening of Wilf-Zeilberger Theory`): print(`By Mohamud Mohammed and Doron Zeilberger`): 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(): ezra:=proc() if args=NULL then print(`qMultiZeilberger`): print(`A Maple package for `): print(`finding recurrences for hypergeometric multi-sums`): print(`Using Zeilberger's Extension to q-Multi-sums of Zeilberger's `): print(`q-Zeilberger algorithm.`): print(``): print(`NOTE: q is a GLOBAL variable, do not use it for anything else!`): print(`qfac(n) is short for (1-q)(1-q^2)...(1-q^n)`): print(`qbin(n,k) is short for qfac(n)/qfac(k)/qfac(n-k) `): print(`qrf(k,n) is short for (1-q^k)(1-q^(k+1))...(1+q^(k+n-1))`): print(``): print(`For help with a specific procedure, type "ezra(procedure_name);"`): print(`Contains procedures: qMulZeil , qMulZeil1, qMulZeilL`): print(): fi: if nargs=1 and args[1]=qMulZeil1 then print(`qMulZeil1(F,POL,qQuad,kList,n,N,para): Given a`): print(`pure q-hypergeometric term F`): print(`in the variable n and the variables of the list kList,`): print(`a symbol N, designating the shift operator in the variable n`): print(`POL, a polynomial in q^n and the q^kList[i]`): print(`and qQuad, a quadratic polynomial in n and the variables of kList`): print(`and a set of auxiliary parameters para, of other variables`): print(`(para may always be set to be the empty set {} )`): print(`finds a recurrence operator ope(n,N) and a list of discrete functions`): print(`that are the pre-certificates (which are polynomials in q^kList[i]'s)`): print(`For example, try:`): print(`Type: F:=qfac(n)/qfac(i)/qfac(j)/qfac(n-i-j)`): print(`qMulZeil1(F,1,0,[i,j],n,N,{}); and `): print(`qMulZeil1(F,q^i+q^j,binomial(i,2)+binomial(j,2),[i,j],n,N,{}); and `): fi: if nargs=1 and args[1]=qMulZeil then print(`qMulZeil(F,POL,qQuad,kList,n,N,para): Given a`): print(`pure q-hypergeometric term F`): print(`in the variable n and the variables of the list kList,`): print(`a symbol N, designating the shift operator in the variable n`): print(`POL, a polynomial in q^n and the q^kList[i]`): print(`and qQuad, a quadratic polynomial in n and the variables of kList`): print(`and a set of auxiliary parameters para, of other variables`): print(`(para may always be set to be the empty set {} )`): print(`finds a recurrence operator ope(n,N) and a list of discrete functions`): print(`that are the certificates (which are polynomials in q^kList[i]'s)`): print(`i.e. a list of length nops(kList) (let's call it r) [R1,.., Rr]`): print(`such that ope(n,N)(POL*q^qQuad*F(n,kList))=`): print(`(K_1-1)(R1*F)+...+(K_r-1)(Rr*F)`): print(`For example, try:`): print(`Type: F:=qfac(n)/qfac(i)/qfac(j)/qfac(n-i-j)`): print(`qMulZeil(F,1,0,[i,j],n,N,{}); and `): print(`qMulZeil(F,q^i+q^j,binomial(i,2)+binomial(j,2),[i,j],n,N,{}); and `): fi: if nargs=1 and args[1]=qMulZeilL then print(`qMulZeilL(F,POL,qQuad,kList,kLimits,n,N,para)`): print(`First does what qMulZeil does, but then tries to find`): print(` a smaller-order recurrence, empirically`): print(`by summing F over the kLimits (i.e. it is tacitly assumed that`): print(`F vanishes outside the summation region`): print(`If there is no lower-order recurrence, the last output is 1`): print(`otherwise, the first output is the empirically-found recurrence`): print(`the last output is the operator that if you left-multiply it`): print(`by the first output you would get the operator given by qMulZeil`): print(`and the second output is like in qMulZeil`): print(`For example, try:`): print(`Type: F:=qfac(n)/qfac(i)/qfac(j)/qfac(n-i-j)`): print(`qMulZeilL(F,0,[i,j],[[0,n],[0,n-i]],n,N,{}); and `): print(`F:=qfac(i+j)*qfac(n-i)*qfac(n-j)/qfac(i)^2/qfac(j)^2/qfac(n-i-j)^2;`): print(`qMulZeilL(F,1,0,[i,j],[[0,n],[0,n-i]],n,N,{}); `): fi: end: qrf1:=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) local j: convert([seq( qrf1(Num[j]+1,i*coeff(Num[j],n,1)),j=1..nops(Num))],`*`)* convert([seq( qrf1( subs(n=n+i,Den[j])+1,(L-i)*coeff(Den[j],n,1)),j=1..nops(Den))],`*`): end: c:=proc(Num,Den,POL,qQuad,e,n,L) local i:convert( [seq(e[i]*q^(expand(subs(n=n+i,qQuad)-qQuad))* subs(n=n+i,POL)*RatioH(Num,Den,n,L,i),i=0..L)],`+`): end: a:=proc(Num,Den,zList,qQuad,kList,n,L) local j,gu,i: for i from 1 to nops(kList) do gu[i]:=zList[i]: od: for j from 1 to nops(Num) do for i from 1 to nops(kList) do if coeff(Num[j],kList[i],1) > 0 then gu[i]:=gu[i]*qrf1(Num[j]+1,coeff(Num[j],kList[i],1)): fi: od: od: for j from 1 to nops(Den) do for i from 1 to nops(kList) do if coeff(Den[j],kList[i],1) < 0 then gu[i]:=gu[i]*qrf1(subs({n=n+L,kList[i]=kList[i]+1}, Den[j])+1,-coeff(Den[j],kList[i],1)): fi: od: od: for i from 1 to nops(kList) do gu[i]:=gu[i]*q^(expand(subs(kList[i]=kList[i]+1,qQuad)-qQuad)): od: [seq(gu[i],i=1..nops(kList))]: end: b:=proc(Num,Den,kList,n,L) local i,j,gu: for i from 1 to nops(kList) do gu[i]:=1: od: for j from 1 to nops(Num) do for i from 1 to nops(kList) do if coeff(Num[j],kList[i],1) < 0 then gu[i]:=gu[i]* qrf1(subs(kList[i]=kList[i]+1,Num[j])+1,-coeff(Num[j],kList[i],1)): fi: od: od: for j from 1 to nops(Den) do for i from 1 to nops(kList) do if coeff(Den[j],kList[i],1) > 0 then gu[i]:=gu[i]*qrf1(subs({n=n+L},Den[j])+1,coeff(Den[j],kList[i],1)): fi: od: od: [seq(gu[i],i=1..nops(kList))]: end: Equ1:=proc(F1,POL,qQuad,kList,n,L,e,X1,X,t) local Num,Den,zList,aList, bList,i: Num:=F1[1]: Den:=F1[2]: zList:=F1[3]: aList:=a(Num,Den,zList,qQuad,kList,n,L): bList:=b(Num,Den,kList,n,L): qkTot(add( aList[i]*X1[i]- subs(kList[i]=kList[i]-1,bList[i])*X[i],i=1..nops(kList))- c(Num,Den,POL,qQuad,e,n,L),kList,t): end: 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: #IV1(d,n): all the vectors of non-negative integers of length d whose #sum is n IV1:=proc(d,n) local gu,i,gu1,i1: if n<0 then gu:=IV1(d,-n): RETURN({seq([seq(-gu[i][i1],i1=1..nops(gu[i]))],i=1..nops(gu))}): fi: if d=0 then if n=0 then RETURN({[]}): else RETURN({}): fi: fi: gu:={}: for i from 0 to n do gu1:=IV1(d-1,n-i): gu:=gu union {seq([op(gu1[i1]),i],i1=1..nops(gu1))}: od: gu: end: #IV(d,n): all the vectors of integers of length d whose #sum is between n[1] and n[2] IV:=proc(d,n) local i: {seq(op(IV1(d,i)),i=n[1]..n[2])}: end: #GenPol(kList,a,deg): The generic polynomial in kList of #ldegree deg[1], and degree deg[2] #using the indexed variable a, followed by the set #of coeffs. GenPol:=proc(kList,a,deg) local gu,i,i1: gu:=IV(nops(kList),deg): convert([seq(a[i]*convert([seq(kList[i1]^gu[i][i1],i1=1..nops(kList))],`*`), i=1..nops(gu))],`+`),{seq(a[i],i=1..nops(gu))}: end: #Extract1(POL1,kList): extracts all the coeffs. of a POL in kList Extract1:=proc(POL1,kList) local k1,kList1,i,POL: POL:=expand(POL1): k1:=kList[nops(kList)]: kList1:=[op(1..nops(kList)-1,kList)]: if nops(kList)=1 then RETURN({seq(coeff(POL,k1,i),i=ldegree(POL,k1)..degree(POL,k1))}): else RETURN({seq( op(Extract1(coeff(POL,k1,i),kList1)), i=ldegree(POL,k1)..degree(POL,k1))}): fi: end: Z1Ld:=proc(F,POL,qQuad,kList,n,L,N,degX,para) local e,x,j,X,gu,var1, var,ope,eq,zList,ku,i,ope1,kap, i1,Sols,j1,Sols1,X1,Sols2,kvar,Y,Y1,eqTry,varTry,t,T,kur,i2: global t0: T:=[seq(t[i1],i1=1..nops(kList))]: ope:=convert([seq(e[j]*N^j,j=0..L)],`+`): zList:=F[3]: var1:=[seq(e[j],j=0..L)]: X:=[]: for i from 1 to nops(kList) do ku:=GenPol(T,x[i],degX[i]): X:=[op(X),ku[1]]: var1:=[op(var1), op(ku[2])]: od: gu:=Equ1(F,POL,qQuad,kList,n,L,e,Y1,Y,t) : for i from 1 to nops(kList) do gu:=subs({Y1[i]=subs(t[i]=q*t[i],X[i]),Y[i]=X[i]},gu): od: eq:=Extract1(gu,T): eqTry:=eq: eqTry:=subs({q=3/2, n=11},eq): for i from 1 to nops(zList) do if type(zList[i],symbol) then eqTry:=subs(zList[i]=5/(i^2+3),eqTry): fi: od: for i from 1 to nops(para) do eqTry:=subs(para[i]=7/(i^3+5),eqTry): od: varTry:=LS(eqTry,var1): if varTry={} or {seq([op(1..L+1,varTry[i])],i=1..nops(varTry))}={[0$(L+1)]} then RETURN(FAIL): fi: print(`The time is:`, time()-t0): print(`There is great hope for a recurrence of order`,L,` please be patient`): var:=LS(eq,var1): Sols:={ seq( [subs({seq(var1[j]=var[i1][j],j=1..nops(var1))},ope), [seq(subs({seq(var1[j]=var[i1][j],j=1..nops(var1))},X[j1]),j1=1..nops(X))]] ,i1=1..nops(var))}: Sols1:={}: for i from 1 to nops(Sols) do ope1:=normal(Sols[i][1]): X1:=normal(Sols[i][2]): if ope1<>0 then kap:=yafe(ope1,N): ope1:=kap[2]: kur:=[seq(normal(X1[i1]/kap[1]),i1=1..nops(X1))]: for i2 from 1 to nops(kList) do kur:=subs(t[i2]=q^kList[i2],kur): od: Sols1:=Sols1 union {[ope1,kur]}: fi: od: if Sols1={} then RETURN(FAIL) fi: if nops(Sols1)>1 then Sols2:={Sols1[1]}: kvar:={Sols1[1][1]}: for i from 2 to nops(Sols1) do if not member(Sols1[i][1],kvar) then kvar:=kvar union {Sols1[i][1]}: Sols2:=Sols2 union {Sols1[i]}: fi: od: if nops(Sols2)>1 then print(`There are more than one operator`): else print(`There is just one operator but several certificates`): fi: print(`The whole thing took:`, time()-t0): RETURN(Sols2): else print(`The whole thing took:`, time()-t0): RETURN(Sols1): fi: end: Z1L:=proc(F,POL,qQuad,kList,n,L,N,para) local degX,gu,i,j,aList,bList,c1,Num,Den,zList,e,t: Num:=F[1]: Den:=F[2]: zList:=F[3]: aList:=a(Num,Den,zList,qQuad,kList,n,L): bList:=b(Num,Den,kList,n,L): c1:=c(Num,Den,POL,qQuad,e,n,L): aList:=[seq(qkTot(aList[i],kList,t),i=1..nops(aList))]: bList:=[seq(qkTot(bList[i],kList,t),i=1..nops(aList))]: c1:=qkTot(c1,kList,t): degX:=FindDeg(aList,bList,c1,t): gu:=Z1Ld(F,POL,qQuad,kList,n,L,N,degX,para): if gu<>FAIL then RETURN(gu): fi: for j from 1 to 1 do degX:=[seq([degX[i][1]-1,degX[i][2]+1],i=1..nops(degX))]: gu:=Z1Ld(F,POL,qQuad,kList,n,L,N,degX,para): if gu<>FAIL then RETURN(gu): fi: od: FAIL: end: Z1:=proc(F,POL,qQuad,kList,n,N,para) local L,gu: gu:=Z1L(F,POL,qQuad,kList,n,0,N,para): for L from 1 while gu=FAIL do gu:=Z1L(F,POL,qQuad,kList,n,L,N,para): od: gu: end: ####converting to the input needed by Z1### qFindFactorial:=proc(Prod1) local Prod2,i,lu,lu1,lu2: Prod2:=factor(Prod1): if type(Prod2,function) and op(0,Prod2)=qfac then RETURN(Prod1): fi: for i from 1 to nops(Prod2) do lu:=op(i,Prod2): if type(lu,function) and op(0,lu)=qfac then RETURN(lu): fi: if type(lu,`^`) then lu1:=op(1,lu): lu2:=op(2,lu): if type(lu1,function) and op(0,lu1)=qfac then if not type(lu2,integer) and lu2>0 then RETURN(FAIL): else RETURN(lu1): fi: fi: fi: od: FAIL: end: FindPower:=proc(Prod1,kList) local Prod2,i,lu: Prod2:=factor(Prod1): if Prod2=1 then RETURN(FAIL): fi: if type(Prod2,`^`) and degree(op(2,Prod1),{op(kList)})=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),{op(kList)})=1 then RETURN(lu): fi: od: FAIL: end: hafokh:=proc(F,kList,n) local F1,Num,Den,POL,z,F1t,F1b,lu,lu1,lu2,i: F1:=F: F1t:=factor(numer(F1)): F1b:=factor(denom(F1)): Num:=[]: lu:=qFindFactorial(F1t): while lu<>FAIL do lu1:=op(lu): if degree(lu1,{n,op(kList)})>1 or not type(coeff(lu1,n,1),integer) then RETURN(FAIL): fi: for i from 1 to nops(kList) do if not type(coeff(lu1,kList[i],1),integer) then RETURN(FAIL): fi: od: Num:=[op(Num),lu1]: F1t:=normal(F1t/lu): lu:=qFindFactorial(F1t): od: for i from 1 to nops(kList) do z[i]:=1: od: lu:=FindPower(F1t,kList): while lu<>FAIL do lu1:=op(1,lu): lu2:=op(2,lu): if degree(lu2,{op(kList)})>1 then RETURN(FAIL): else for i from 1 to nops(kList) do z[i]:=z[i]*lu1^coeff(lu2,kList[i],1): F1t:=normal(simplify(normal(F1t/lu1^(kList[i]*coeff(lu2,kList[i]))))): od: fi: lu:=FindPower(F1t,kList): od: POL:=F1t: Den:=[]: lu:=qFindFactorial(F1b): while lu<>FAIL do lu1:=op(lu): if degree(lu1,{n,op(kList)})>1 or not type(coeff(lu1,n,1),integer) then RETURN(FAIL): fi: for i from 1 to nops(kList) do if not type(coeff(lu1,kList[i],1),integer) then RETURN(FAIL): fi: od: Den:=[op(Den),lu1]: F1b:=normal(F1b/lu): lu:=qFindFactorial(F1b): od: lu:=FindPower(F1b,kList): while lu<>FAIL do lu1:=op(1,lu): lu2:=op(2,lu): if degree(lu2,{op(kList)} ) >1 then RETURN(FAIL): else for i from 1 to nops(kList) do z[i]:=z[i]/lu1^coeff(lu2,kList[i],1): F1b:=normal(simplify(normal(F1b/lu1^(kList[i]*coeff(lu2,kList[i],1))))): od: fi: lu:=FindPower(F1b,kList): 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 degree(POL,{n,op(kList)})=FAIL then RETURN(FAIL): fi: if POL<>1 then print(`Bad input, the first argument must be a pure hyp. term`): RETURN(FAIL): fi: [Num,Den,[seq(z[i],i=1..nops(kList))]]: end: ##End conversion part qMulZeil1:=proc(F,POL,qQuad,kList,n,N,para) local F1,lu: global t0: t0:=time(): F1:=F: F1:=hafokh(F1,kList,n): if F1=FAIL then print(`Bad input`): RETURN(FAIL): fi: lu:=Z1(F1,POL,qQuad,kList,n,N,para): if nops(lu)>1 then print(`There are more than one operator`): RETURN(lu): else RETURN(lu[1]): fi: end: #Linear solver LS:=proc(eq,var) local i,var1,ind,X,j: var1:=SolveTools[Linear](eq,{op(var)}): ind:={}: for i from 1 to nops(var1) do if op(1,var1[i])=op(2,var1[i]) then ind:=ind union {op(1,var1[i])}: fi: od: X:=[seq(expand(subs(var1,var[i])),i=1..nops(var))]: {seq([seq(coeff(X[i],ind[j],1),i=1..nops(X))],j=1..nops(ind))}: end: #End Linear solver #qMulZeil(F,POL,qQuad,kList,n,N,para):The operators plus the certificates qMulZeil:=proc(F,POL,qQuad,kList,n,N,para) local F1, lu,kadima,X,ope,i,L,bList: option remember: lu:=qMulZeil1(F,POL,qQuad,kList,n,N,para): if type(lu,set) then lu:=lu[1]: fi: ope:=lu[1]: L:=degree(ope,N): F1:=hafokh(F,kList,n): bList:=b(F1[1],F1[2],kList,n,L): kadima:=RatioH(F1[1],F1[2],n,L,0): X:=lu[2]: ope,[seq(X[i]/kadima*subs(kList[i]=kList[i]-1,bList[i]),i=1..nops(kList))]: end: #SumSet(n,kList,kLimits,n0): the set of summation given #by kLimits SumSet:=proc(n,kList,kLimits,n0) local i1,j1,gu,gu1,low1,high1,kLimits1,kList1: if nops(kList)<>nops(kLimits) or kList=[] then ERROR(`Bad input`): fi: low1:=subs(n=n0,kLimits[1][1]): high1:=subs(n=n0,kLimits[1][2]): if nops(kList)=1 then RETURN({seq([i1],i1=low1..high1)}): fi: gu:={}: for i1 from low1 to high1 do kLimits1:=subs(n=n0,[op(2..nops(kLimits),kLimits)]): kList1:=[op(2..nops(kList),kList)]: gu1:=SumSet(kList[1],kList1,kLimits1,i1): gu:=gu union {seq([i1,op(gu1[j1])],j1=1..nops(gu1))}: od: gu: end: #qSumF1(F,POL,qQuad,n,kList,kLimits,n0): The sum of F(n,[kList]) at n=n0 #given by the limits indicated by kLimits, in that order #For example, try: #qSumF1(qfac(i+j)*qfac(n-i)*qfac(n-j)/qfac(i)^2/qfac(j)^2/qfac(n-i-j)^2, #1,0,n,[i,j],[[0,n],[0,n-i]],5); qSumF1:=proc(F,POL,qQuad,n,kList,kLimits,n0) local Domain1,i1,j1,F1: F1:=subs({qfac=qfac1,qbin=qbin1},F): Domain1:=SumSet(n,kList,kLimits,n0): normal( add(eval( subs({n=n0,seq(kList[i1]=Domain1[j1][i1],i1=1..nops(kList))},F1*POL* q^qQuad)), j1=1..nops(Domain1))): end: #qSumF(F,POL,qQuad,n,kList,kLimits,N0): The first N0+1 terms #of sum of F(n,[kList]) at n=n0 #given by the limits indicated by kLimits, in that order #For example, try: #qSumF(qfac(i+j)*qfac(n-i)*qfac(n-j)/qfac(i)^2/qfac(j)^2/qfac(n-i-j)^2,1,0, #n,[i,j],[[0,n],[0,n-i]],5); qSumF:=proc(F,POL,qQuad,n,kList,kLimits,N0) local n0: [seq(qSumF1(F,POL,qQuad,n,kList,kLimits,n0),n0=1..N0)]: end: ###qfindrec (from SCHUTZENBERGER) qfindrec:=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]*q^(n*j)*N^i: var:=var union {a[i,j]}: od: od: eq:={}: for n0 from 1 to (1+DEGREE)*(1+ORDER)+4 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(`Too much slack`): fi: #for i from 1 to nops(kv) do # ope:=subs(op(i,kv)=1,ope): #od: if ope<>0 then RETURN(yafe(ope,N)[2]): else RETURN(0): fi: end: ###End qfindrec (from SCHUTZENBERGER) #ApplyOpe(ope,n,N,gu): applies the operator ope(n,N) to the sequence gu ApplyOpe:=proc(ope,n,N,gu) local i,i1: [seq(add(subs(n=i,coeff(ope,N,i1))*gu[i+i1],i1=0..degree(ope,N)), i=1..nops(gu)-degree(ope,N))]: end: #OpeToSeq(ope,n,N,Ini,n0,n1): Given an operator ope(n,N) and #Initial conditions [a[n0],a[n0+1], ... a[n0+L-1]] (L=degree(ope,N) #is the order of ope, outputs the sequence up to n=n0+n1-1 OpeToSeq:=proc(ope,n,N,Ini,n0,n1) local gu,i,ope1,L,j: ope1:=expand(normal((subs(n=n-ldegree(ope,N),ope)/N^ldegree(ope,N)))): L:=degree(ope1,N): if nops(Ini)<>L then ERROR(`The size of Ini should be`, L , ` the order of`, ope): fi: ope1:=normal(ope/coeff(ope,N,L)): ope1:=expand(normal(subs(n=n-L,ope1)/N^L)): gu:=Ini: for i from n0+L to n0+n1-1 do gu:=[op(gu),-add(subs(n=i,coeff(ope1,N,-j))*gu[nops(gu)-j+1],j=1..L)]: od: gu: end: #qMulZeilL(F,POL,qQuad,,kList,kLimits,n,N,para): first applies qMulZeil, then #tries to find a better operator (i.e. of smaller #order and complexity). #if successful, returns the original [ope,cert] followed by #the better operator. If it couldn't find a better operator #it returns, [ope,cert],0 #For example, try: #qMulZeilL(F,POL,qQuad,qfac(n)/qfac(i)/qfac(j)/ #qfac(n-i-j),[i,j],[[0,n],[0,n-i]],n,N): qMulZeilL:=proc(F,POL,qQuad,kList,kLimits,n,N,para) local lu,L,Ini,ORDER,DEGREE,Comp1,lu1,gu,ope,lu2,t,ope1: lu:=qMulZeil(F,POL,qQuad,kList,n,N,para): ope:=lu[1]: L:=degree(ope,N): ope1:=numer(normal(ope)): ope1:=expand(ope1): ope1:=subs(q^n=t,ope1): Comp1:=(L+1)*(degree(ope1,t)+1): Ini:=qSumF(F,POL,qQuad,n,kList,kLimits,L): gu:=OpeToSeq(ope,n,N,Ini,1,(L+2)*(degree(ope1,t)+2)+3): for ORDER from 0 to L-1 do for DEGREE from 0 to trunc(Comp1/(ORDER+1)) do lu1:=qfindrec(gu,DEGREE,ORDER,n,N): if lu1<>0 then lu2:=Khalek(ope,lu1,n,N): if lu2=FAIL then print(`Something is wrong`): RETURN(FAIL): else RETURN(lu1,lu[2],lu2): fi: fi: od: od: lu,1: end: #Khalek(ope1,ope2,n,N): given recurrence #operators ope1 and ope2, both monic in #finds the operators ope3 such that ope1=ope3*ope2 #or returns FAIL Khalek:=proc(ope1,ope2,n,N) local ope3,d,lu,mu: if degree(ope1,N)=degree(ope2,N) then if ope1=ope2 then RETURN(1): else RETURN(FAIL): fi: fi: d:=degree(ope1,N)-degree(ope2,N): ope3:=expand(ope1-N^d*subs(n=n+d,ope2)): lu:=yafe(ope3,N): mu:=Khalek(lu[2],ope2,n,N): if mu=FAIL then RETURN(FAIL): fi: N^d+lu[1]*mu: end: #qkTot(Pol,kList,t): Given a polynomial Pol in q^kList[i] #subsitutes q^k List[i] by t[i] qkTot:=proc(Pol,kList,t) local P1,i: P1:=expand(Pol): for i from 1 to nops(kList) do P1:=subs(q^kList[i]=t[i],P1): od: P1: end: FindDeg:=proc(aList,bList,c,t) local gu,r,T,i: r:=nops(aList): T:={seq(t[i],i=1..r)}: gu:=[]: for i from 1 to nops(aList) do gu:=[op(gu), [ldegree(expand(c),T)-min(ldegree(expand(aList[i]),T), ldegree(expand(bList[i]),T)), degree(expand(c),T)-max(degree(expand(aList[i]),T), degree(expand(bList[i]),T))]]: od: gu: end: qbin:=proc(n,k):qfac(n)/qfac(k)/qfac(n-k):end: qfac1:=proc(n) local i: convert([seq(1-q^i,i=1..n)],`*`):end: qbin1:=proc(n,k): expand(normal(qfac1(n)/qfac1(k)/qfac1(n-k))):end: qrf:=proc(k,n): if type(k, numeric) and k<=0 then 0: elif k=1 then qfac(n): else qfac(n+k-1)/qfac(k-1): fi: end: