# The following examples are presented as is. They accompany the # paper ``Non-commutative Elimination in Ore Algebras Proves # Multivariate Identities'', by Frederic Chyzak and Bruno Salvy. # # Some results of this session may differ from those of the paper due # to changes in the program. For the reason of changes and new bugs, # I also did not reproduce all the examples of the session. This will # be done in the future. # # Frederic Chyzak, May 14, 1997. # Modified (minor corrections), May 27, 1997. # Change this according to your configuration. libname:=`/home/musigny/algo/chyzak/ArchivesMgfun/Lib1.3-released`,libname: with(Mgfun): # Patch a bug in this version. `type/OreAlgebra`:=readlib(`type/Mgfun/OreAlgebra`): # Section 1.4.1. Jacobi polynomials A:=Ore_algebra(comm={a,b},shift=[Sn,n],diff=[Dx,x]): G:=[2*(n+2)*(n+a+b+2)*(2*n+a+b+2)*Sn^2 -((2*n+a+b+3)*(a^2-b^2)+(2*n+a+b+2) *(2*n+a+b+3)*(2*n+a+b+4)*x)*Sn +2*(n+a+1)*(n+b+1)*(2*n+a+b+4), (2*n+a+b+2)*(1-x^2)*Dx*Sn-(n+1) *(a-b-(2*n+a+b+2)*x)*Sn-2*(n+a+1)*(n+b+1)]: skewelim(G[1],G[2],Sn,A); # Section 1.4.2. Gauss's hypergeometric function A:=Ore_algebra(comm={b,c},diff=[Dz,z],shift=[Sa,a]): P:=z*(1-z)*Dz^2+(c-(a+b+1)*z)*Dz-a*b: H:=a*Sa-(z*Dz+a): skewelim(P,H,Dz,A); A:=Ore_algebra(comm={a,b,c},diff=[Dz,z]): GCD:=skewgcdex(P,z*Dz+a,Dz,A): B:=collect((a-1)*subs(a=a-1,GCD[3]/GCD[1]),Dz,factor); # Section 1.4.3. Partially hypergeometric series A:=shift_algebra([Sn,n],[Sk,k],[comm,z]): h:=(-1)^k*binomial(n,k)^2*binomial(n+k,k)^2; G:=hypergeomtoholon(h,A): G:=collect(G,{Sn,Sk},factor); L:=op(select(has,G,Sk)); # You can run this if you have gfun in the share library. # #with(share): #readshare(gfun,analysis): #with(gfun): #applyopr(L,u(k),A); #rectodiffeq({"},u(k),f(z)); #M:=collect(subs({seq((D@@i)(f)(z)=Dz^i,i=0..4)},"),Dz,factor); # # Otherwise, here is the result. M:=z^3*(z+1)*Dz^4+2*z^2*(4*z+3)*Dz^3-z*(-14*z+2*z*n+2*z*n^2-7)*Dz^2 +(-4*z*n^2+4*z-4*z*n+1)*Dz+(n+1)^2*n^2; H:=Sn; # Naive method that returns an operator of order 7. normal(applyopr(H,h,A)/h,expanded); P:=numer("): Q:=denom(""): A:=Ore_algebra(diff=[Dz,z],shift=[Sn,n]): P2:=collect(subs({seq(k^i=oppower(z*Dz,i,A),i=1..degree(P,k))},P),Dz,factor); Q2:=collect(subs({seq(k^i=oppower(z*Dz,i,A),i=1..degree(Q,k))},Q),Dz,factor); collect(Q2*H-P2,{Dz,Sn},distributed,factor); # The contiguity relation of order 7 is: C[7]:=collect(skewelim(",M,Dz,A),Sn,factor); # (This needs several minutes.) # Improved method that returns an operator of order 3. GCD:=subs(n=n+1,skewgcdex(subs(n=n-1,Q2),M,Dz,A)): U:=collect(GCD[2],Dz,factor); factor(GCD[1]); # We reduce U&*P/GCD[1] modulo M. PR:=collect(opprod(U,P2,A),Dz,factor): PDIV:=skewpdiv(PR,M,Dz,A): # Therefore, H=Sn is R:=collect(PDIV[3]/PDIV[1]/GCD[1],Dz,factor); # (Fraction-free) Gaussian elimination: T[1]:=[1,1]: T[Sn]:=[numer(R),denom(R)]: T[Sn^2]:=[opprod(subs(n=n+1,"[1]),numer(R),A),subs(n=n+1,"[2])*denom(R)]: T[Sn^3]:=[opprod(subs(n=n+1,"[1]),numer(R),A),subs(n=n+1,"[2])*denom(R)]: T[Sn^4]:=[opprod(subs(n=n+1,"[1]),numer(R),A),subs(n=n+1,"[2])*denom(R)]: skewpdiv(T[Sn^2][1],M,Dz,A): T[Sn^2]:=["[3],T[Sn^2][2]*"[1]]: skewpdiv(T[Sn^3][1],M,Dz,A): T[Sn^3]:=["[3],T[Sn^3][2]*"[1]]: skewpdiv(T[Sn^4][1],M,Dz,A): T[Sn^4]:=["[3],T[Sn^4][2]*"[1]]: Z:=collect(add(eta[i]*T[Sn^i][1]/T[Sn^i][2],i=0..4),Dz,numer@normal): SOL:=solve({coeffs(Z,Dz)},{seq(eta[i],i=0..4)}): # The contiguity relation of order 4 is: C[4]:=collect(primpart(subs(SOL,add(eta[i]*Sn^i,i=0..4)),Sn),Sn,factor); # Section 1.4.4 Sylvester's dialytic elimination A:=Ore_algebra(shift=[Sn,n],shift=[Sm,m],comm={alpha,x},polynom=m); skewelim(op(hypergeomtoholon( (-1)^m*GAMMA(alpha+n-m)/m!/(n-2*m)!*(2*x)^(n-2*m),A)),m,A); map(collect,series(",Sm=1),Sn,factor); # Section 1.5. First example (Jacobi polynomials) A:=Ore_algebra(comm={a,b},shift=[Sn,n],diff=[Dx,x]): G:=[2*(n+2)*(n+a+b+2)*(2*n+a+b+2)*Sn^2 -((2*n+a+b+3)*(a^2-b^2)+(2*n+a+b+2) *(2*n+a+b+3)*(2*n+a+b+4)*x)*Sn +2*(n+a+1)*(n+b+1)*(2*n+a+b+4), (2*n+a+b+2)*(1-x^2)*Dx*Sn-(n+1) *(a-b-(2*n+a+b+2)*x)*Sn-2*(n+a+1)*(n+b+1)]: T:=termorder(A,plex=[Sn,Dx]): remove(has,gbasis(G,T),Sn); # Section 2.2.2. Example of addition of partial finite functions # With rectangular systems. A:=Ore_algebra(diff=[Dx,x],diff=[Dy,y],comm={mu,nu}): T:=termorder(A,tdeg=[Dx,Dy]): GL[f]:=gbasis([Dx-mu,Dy-nu],T); GL[g]:=gbasis([x^2*Dx^2+x*Dx+x^2-mu^2,y^2*Dy^2+y*Dy+y^2-nu^2],T); # Patch a bug in this version. A[indet]:=A[all_indets]: collect(dependency(f(x,y)+g(x,y),x,3,GL,T),Dx,factor); collect(dependency(f(x,y)+g(x,y),y,3,GL,T),Dy,factor); # Using the FGLM algorithm: collect(hsum([[GL[f],T],[GL[g],T]],T),{Dx,Dy},distributed,factor); # By intersection of ideals. A:=Ore_algebra(diff=[Dx,x],diff=[Dy,y],comm={mu,nu,t}): T:=termorder(A,lexdeg=[[t],[Dx,Dy]]): G:=[t*(Dx-mu),t*(Dy-nu),(1-t)*(x^2*Dx^2+x*Dx+x^2-mu^2), (1-t)*(y^2*Dy^2+y*Dy+y^2-nu^2)]; GB:=gbasis(G,T): collect(remove(has,GB,t),{Dx,Dy},distributed,factor); # Section 3.2. Example of creative telescoping by Groebner bases A:=Ore_algebra(comm={a,b},shift=[Sn,n],diff=[Dx,x],diff=[Dy,y]): G:=[2*(n+2)*(n+a+b+2)*(2*n+a+b+2)*Sn^2 -y*(2*n+a+b+3)*(a^2-b^2+4*x*n^2+4*x*n*a+4*x*n*b +12*x*n+x*a^2+2*x*a*b+6*x*a+x*b^2+6*x*b+8*x)*Sn +2*(n+a+1)*(n+b+1)*(2*n+a+b+4)*y^2, -2*(n+a+1)*(n+b+1)*y+(n+1)*(-a+b+2*x*n+x*a+x*b+2*x)*Sn -(x-1)*(x+1)*(2*n+a+b+2)*Dx*Sn, n*(n+a+b+1)+(b-a-x*a-x*b-2*x)*Dx-(x-1)*(x+1)*Dx^2, y*Dy-n]: T:=termorder(A,lexdeg=[[n],[Sn,Dx,Dy]]): GB:=gbasis(G,T): CT:=collect(subs(Sn=1,GB[2..6]),[Dx,Dy],distributed,factor); R:=sqrt(1-2*x*y+y^2): P:=1/(R*(1-y+R)^a*(1+y+R)^b): map(simplify,map(applyopr,CT,P,A)); # Section 3.3. Extension of Takayama's algorithm for definite integral # to definite anti-partial hdefsum(G,n=vanishing,A): map(collect,",{Dx,Dy},distributed,factor); # The algebra defined by # A:=orealg(comm=[a,b,n],diff=[Dx,x],diff=[Dy,y]): # is used internally only. # Section 3.4.1. An identity between Franel and Apery numbers # Operator for the left-hand side. A:=shift_algebra([Sn,n],[Sk,k]): P1:=collect(op( hdefsum(hypergeomtoholon(binomial(n,k)^2*binomial(n+k,k)^2,A), k=vanishing,A)),Sn,factor); # Operator for the right-hand side. B:=shift_algebra([Sk,k],[Sj,j]): hdefsum(hypergeomtoholon(binomial(k,j)^3,B), j=vanishing,B); T:=termorder(A,tdeg=[Sn,Sk]): [gbasis([op(""),Sn-1],T),T]; [gbasis(hypergeomtoholon(binomial(n,k)*binomial(n+k,k),A),T),T]; P2:=collect(op(hdefsum(hprod(["","],T),k=vanishing,A)),Sn,factor); # Proof of the identity GCD:=skewgcdex(P1,P2,Sn,A); collect(opprod(GCD[4],P1,A),Sn,factor); seq(add(binomial(n,k)^2*binomial(n+k,k)^2,k=0..n),n=0..3); seq(add(binomial(n,k)*binomial(n+k,k)* add(binomial(k,j)^3,j=0..k),k=0..n),n=0..3); # Therefore both sides are equal and satisfy collect(GCD[1],Sn,factor); # Section 3.4.2. A Rogers-Ramanujan identity A:=Ore_algebra(comm=q,qshift=[Sn,qn,q],qshift=[Sk,qk,q]): h:=q^(k^2)/qfactorial(k)/qfactorial(n-k); G:=collect(numer(subs({q^n=qn,q^k=qk},normal( [Sn-subs(n=n+1,h)/h,Sk-subs(k=k+1,h)/h],expanded))),{Sn,Sk},factor); collect(op(hdefqsum(expand(G),qk=vanishing,A)),Sn,factor);