# 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 paper.  This will
# be done in the future.
#
# Frederic Chyzak, July 8, 1997.

# Change this according to your configuration.
libname:=`/export/chyzak/ArchivesMgfun/Lib2.0-5.4-released`,libname:
with(Ore_algebra):
with(Groebner):

# Section 1.4.1. Jacobi polynomials
A:=skew_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)]:
skew_elim(G[1],G[2],Sn,A);

# Section 1.4.2. Gauss's hypergeometric function
A:=skew_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):
skew_elim(P,H,Dz,A);
A:=skew_algebra(comm=[a,b,c],diff=[Dz,z]):
GCD:=skew_gcdex(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:=numer(normal({Sn-subs(n=n+1,h)/h,Sk-subs(k=k+1,h)/h},expanded)):
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:=skew_algebra(diff=[Dz,z],shift=[Sn,n]):
P2:=collect(subs({seq(k^i=skew_power(z*Dz,i,A),
    i=1..degree(P,k))},P),Dz,factor);
Q2:=collect(subs({seq(k^i=skew_power(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(skew_elim(",M,Dz,A),Sn,factor);
# (This needs several minutes.)
# Improved method that returns an operator of order 3.
GCD:=subs(n=n+1,skew_gcdex(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(skew_product(U,P2,A),Dz,factor):
PDIV:=skew_pdiv(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]:=[skew_product(subs(n=n+1,"[1]),numer(R),A),subs(n=n+1,"[2])*denom(R)]:
T[Sn^3]:=[skew_product(subs(n=n+1,"[1]),numer(R),A),subs(n=n+1,"[2])*denom(R)]:
T[Sn^4]:=[skew_product(subs(n=n+1,"[1]),numer(R),A),subs(n=n+1,"[2])*denom(R)]:
skew_pdiv(T[Sn^2][1],M,Dz,A):
T[Sn^2]:=["[3],T[Sn^2][2]*"[1]]:
skew_pdiv(T[Sn^3][1],M,Dz,A):
T[Sn^3]:=["[3],T[Sn^3][2]*"[1]]:
skew_pdiv(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:=shift_algebra([Sn,n],[Sm,m],[comm,alpha],[comm,x],polynom=m):
(-1)^m*GAMMA(alpha+n-m)/m!/(n-2*m)!*(2*x)^(n-2*m):
numer(normal([Sn-subs(n=n+1,")/",Sm-subs(m=m+1,")/"],expanded)):
skew_elim(op("),m,A);
map(collect,series(",Sm=1),Sn,factor);

# Section 1.5. First example (Jacobi polynomials)
A:=skew_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
# The methods by rectangular systems or by the FGLM algorithm are not
# directly available in Mgfun2.0, so that they are not demonstrated.
# The method by intersection of ideals is as follows.
A:=diff_algebra([Dx,x],[Dy,y],[comm,mu],[comm,nu],[comm,t],polynom=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:=skew_algebra(comm=[a,b],shift=[Sn,n],diff=[Dx,x],diff=[Dy,y],polynom=n):
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,remove(has,GB,n)),[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));

# The functions hdefsum, hdefqsum and hdefint have not been
# reimplemented in version 2.0, and the examples of the following
# sections are not demonstrated here.
# Section 3.3. Extension of Takayama's algorithm for definite integral
# Section 3.4.1. An identity between Franel and Apery numbers
# Section 3.4.2. A Rogers-Ramanujan identity