Ore_algebra, a Package for Skew Operators
by Frederic Chyzak, Algorithms Project, INRIA, France
Frederic.Chyzak@inria.fr
Several algorithms for integration and summation have a natural description in terms of linear differential and difference operators, which in turn are well described by skew (or Ore) polynomials. This was the starting point for the Ore_algebra package.
> with(Ore_algebra);
To work with an Ore algebra, we first have to declare it. The package creates a table which implements and remembers the
operations in this algebra. Here is the example of the algebra of linear differential operators in the differential operator
with (rational) coefficients in
.
> A:=diff_algebra([Dx,x]);
(This is the name of the table.)
Although the usual product in Maple is commutative, we use
to denote skew products with the convention that powers of
are on the right.
Multiplication of operators is obtained by the function skew_product.
> P1:=x*Dx^2-1: P2:=Dx-x: P1_P2:=skew_product(P1,P2,A);
Remember that in skew algebras of linear operators, factorizations are seldom unique.
> DEtools[DFactor](P1_P2,[Dx,x]);
(Univariate) Ore polynomial rings have a Euclidean algorithm. This is implemented in the following routines: skew_pdiv , skew_prem , skew_gcdex and skew_elim . An extension to the multivariate case is available.
As a first application, we show how the package can be used to derive a proof of the Reed-Dawson identity:
> U[n]=Sum(binomial(n,k)*binomial(2*k,k)*(-2)^(n-k),k=0..n);
is
> binomial(n,n/2);
for even
and
for odd
.
Let us introduce an algebra of linear recurrence operators in
and
.
and
denote the shift operators in
and
respectively.
> A:=shift_algebra([Sn,n],[Sk,k]):
We compute first order recurrences satisfied by the summand (which we denote by
).
> h:=binomial(n,k)*binomial(2*k,k)*(-2)^(n-k);
> Sn-normal(applyopr(Sn,h,A)/h,expanded);
> Sk-normal(applyopr(Sk,h,A)/h,expanded);
> G:=map(numer,{%%,%}):
We eliminate
between these operators by computing their greatest common right divisor.
> GCD:=skew_gcdex(op(G),k,A);
> series(GCD[1],Sk=1);
The factor
corresponds to taking a finite difference. Definite summation for
over all integers yields a telescoping series which collapses to 0. We thus get:
> Sum(applyopr(coeff(%,Sk-1,0),u[n,k],A),k=-infinity..infinity)=0;
Therefore, the sum
satisfies
> collect(primpart(coeff(%%,Sk-1,0),Sn),Sn,normal);
> applyopr(%,U[n],A);
Now, the announced result follows from initial conditions:
For
,
> U[1]:=add(subs(n=1,h),k=0..1);
so that
for odd
.
For
,
> U[0]:=add(subs(n=0,h),k=0..0);
so that
is given by
> 2*(p+1)*V(p+1)-4*(2*p+1)*V(p);
> LREtools[hypergeomsols](%,V(p),{V(0)=1},output=basis);
which is
.
As a more advanced application, we derive contiguity relations for Gauss hypergeometric function. This function is known to Maple as:
> f:=hypergeom([a,b],[c],z);
It is
for
> u:=pochhammer(a,n)*pochhammer(b,n)/pochhammer(c,n)/n!:
When considering the summand, we introduce the following algebra.
> A:=shift_algebra([Sn,n],[Sa,a],comm={b,c}):
The bivariate sequence
vanishes at both following operators:
> G:={(c+n)*(n+1)*Sn-(a+n)*(b+n),a*Sa-(a+n)}:
> normal(map(applyopr,G,u,A),expanded);
(There are general algorithms to find such operators.)
From the previous first order recurrences, we derive relations on
in the mixed differential difference algebra
> A:=skew_algebra(diff=[Dz,z],shift=[Sa,a],comm={b,c}):
The equations are obtained by multiplication of the recurrences by
, followed by summation over all non-negative
. Formally, this corresponds to changing
into
and
into
.
>
map(proc(p)
collect(add(add(coeff(coeff(p,Sn,i),n,j)*
skew_product(skew_power(z*Dz,j,A),1/z^i,A),
j=0..degree(coeff(p,Sn,i),n)),
i=0..degree(p,Sn)),Dz,factor)
end,G);
Therefore, we set:
> P:=z*(1-z)*Dz^2+(c-(a+b+1)*z)*Dz-a*b:
> H:=z*Dz/a+1:
> G:={P,numer(Sa-H)};
The linear differential operator
is called a step-up operator. It relates the forward shift of
to derivatives of
by the following equation
> hypergeom([a+1,b],[c],z)=z/a*diff(hypergeom([a,b],[c],z),z)+hypergeom([a,b],[c],z);
The elimination of
between this step-up operator and the differential equation yields a contiguity relation for
, i.e., a purely recurrence equation. It is obtained by the extended skew gcd algorithm:
> C:=collect(skew_elim(P,numer(Sa-H),Dz,A),Sa,factor);
In other words, Gauss hypergeometric function satisfies the following equation:
> applyopr(C,hypergeom([a,b],[c],z),A)=0;
More interestingly, the extended Euclidean algorithm yields a step-down operator for
, i.e., a relation between an inverse shift of
and its derivatives. This is obtained by computing an inverse of
modulo
.
> GCD:=skew_gcdex(P,numer(H),Dz,A);
From this result, we have
and
. In particular,
> B:=collect((a-1)*subs(a=a-1,GCD[3]/GCD[1]),Dz,factor);
is the step-down operator:
> hypergeom([a-1,b],[c],z)=z*(1-z)*diff(hypergeom([a,b],[c],z),z)/(a-c)+(c-a-z*b)*hypergeom([a,b],[bc],z)/(c-a);
There is also an algorithm to compute (one-sided) Groebner bases in Ore algebras. This will be demonstrated in the presentation of the Groebner package.