# -*-Maple-*-
########################## Ore polynomials #################################
##
## Title: OrePoly
## Created: Sept 1993
## Author: Bruno Salvy
##
##
## Modified june 94. Fixed several bugs reported by Juerg Bolliger,
## op replaced by [].
## Modified july 94. Faster leftquorem, rightquorem, leftrem and rightrem.
## Modified oct 94. Fixed typos in these 4 functions.
## Addition becomes n-ary.
## Added OrePoly/pseudorightquorem.
## Added OrePoly/rightannihilators, OrePoly/leftannihilators.
## Many functions rewritten.
## Faster product in the shift and q-shift case.
##
## Description:
#
# Having a small ring of Ore polynomials is very fashionable.
# The ground field K(x) is passed around all the time. It is a list of three
# elements:
#
# sigma: endomorphism of the ring, which is
# Id in the differential case
# shift in the difference and shift cases
# x --> q*x in the q-difference and q-shift cases
# In practice, sigma takes two arguments:
# a rational function f(x) and an integer k,
# it yields sigma^{(k)}(f(x)), the kth iterate of sigma.
# The integer k may be negative, but sigma^{(-1)} need be
# defined only for operations on the left.
#
# delta: a derivation (but not necessarily)
# It will be d/dx in the differential case,
# f(x+1)-f(x) in the difference case,
# 0 in the shift and q-shift cases,
# f(qx)-f(x) in the q-difference case
# x: the variable.
# It would also be possible to have properties of the ground field handled
# this way, but we want this to look like Maple code for the moment.
#
# Ore polynomials are polynomials in K(x)[Delta] where Delta satisfies
# the following commutation rule:
#
# Delta.X = sigma(X) Delta + delta(X) 1.
#
# Operations on the right will usually work with very few constraints on
# sigma and delta, while operations on the left are more demanding.
#
# Ore polynomials will be represented internally as lists of
# rational functions coefficients.
# OrePoly/degree
#Input: the polynomial (as a list of coefficients)
#Output: the degree of the operator as a polynomial in delta
# In a better Maple, this would be a macro.
`OrePoly/degree`:=proc(f) nops(f)-1 end: # `OrePoly/degree`
# OrePoly/ldeg
#Input: the polynomial
#Output: its valuation
`OrePoly/ldeg`:=proc(f) local i;
for i to nops(f) do if f[i]<>0 then RETURN(i-1) fi od
end: # `OrePoly/ldeg`
# OrePoly/+
#Input: k polynomials
#Output: their sum
# less expensive than repetitive OrePoly/+
`OrePoly/+`:=proc()
local nr,nf,i,j,r;
r:=[0];
nr:=1;
for i to nargs do
nf:=nops(args[i]);
if nf<=nr then
r:=[seq(r[j]+args[i][j],j=1..nf),op(nf+1..nr,r)]
else r:=[seq(r[j]+args[i][j],j=1..nr),op(nr+1..nf,args[i])];
nr:=nf
fi
od;
r:=normal(r);
for i from nops(r) by -1 to 2 do
if r[i]<>0 then break
else r:=subsop(i=NULL,r) fi
od;
r
end: # `OrePoly/+`
# OrePoly/*
#Input: two operators in their list form, the field.
#Output: their composition in the same form
`OrePoly/*`:=proc(f,g,F)
local j,p,n;
if nops(f)=1 then # not really necessary, but this may give some speedup
[seq(f[1]*g[j],j=1..nops(g))]
elif F[2]=0 then # faster in shift case
map(convert,[seq([seq(f[p+1]*F[1](g[n-p+1],p),
p=max(0,n-nops(g)+1)..min(n,nops(f)-1))],n=0..nops(f)+nops(g)-2)],`+`)
else map(convert,[seq([seq(seq(
binomial(p,j)*f[p+1]*F[1]((F[2]@@(p-j))(g[n-j+1]),j),
j=max(0,n-nops(g)+1)..min(p,n)),p=0..nops(f)-1)],
n=0..nops(f)+nops(g)-2)],`+`)
fi
end: # `OrePoly/*`
# OrePoly/print
#Input: the polynomial,
# a name for the differential operator
#Output: A printed form of the polynomial.
# Since we are stuck with commutative Maple, the way it's done is by
# printing a series. The particular domains of gfun will redefine this by
# makerecoper and makediffoper.
`OrePoly/print`:=proc(f,DD) local i;
series(convert([seq(f[i]*DD^(i-1),i=1..nops(f))],`+`),DD,infinity)
end: # `OrePoly/print`
# OrePoly/rightquorem, OrePoly/leftquorem
#Input: both polynomials and the field
#Output: the left (resp. right) quotient and remainder
# F = QA+R (resp. F = AQ+R)
`OrePoly/rightleftquorem`:=proc(f,a,F,right)
local df,da,q,i,rem,j,lc;
df:=`OrePoly/degree`(f); da:=`OrePoly/degree`(a);
if df[0] then RETURN([1]) else RETURN(a) fi fi;
q:=`OrePoly/rightrem`(a,b,F);
dq:=`OrePoly/degree`(q);
da:=db;db:=dq;a:=b;b:=q
od
end: # `OrePoly/rightgcd`
`OrePoly/leftgcd`:=subs(`OrePoly/rightrem`=`OrePoly/leftrem`,
op(`OrePoly/rightgcd`)):
# OrePoly/leftgcdex, OrePoly/rightgcdex
#Input: two polynomials f and g, two names u and v, and the field
#Output: the left (right) gcd d, u and v having been assigned
# so that fu+gv=d (resp. uf+vg=d).
`OrePoly/rightleftgcdex`:=proc(f,g,u,v,F,right,lcm)
local da,db,a,b,q,dq,u0,u1,v0,v1,t,i,qr;
da:=`OrePoly/degree`(f); db:=`OrePoly/degree`(g);
if da``(args,field)`,` `,
`SYNOPSIS: `,
`- Ore polynomials are polynomials in two non-commutative variables, one of`,`\
which is a quasi-derivative. The ring of Ore polynomials provides a convenient`
,
`generalization of linear differential, difference and q-difference operators.`
,`In this ring, addition is the usual addition of operators, while product is`,
`composition of operators.`,` `,
`- All the functions of the package take a field as last argument. This field`,
`is given as a list of three elements:`,
` . sigma: sigma(.,1) is an endomorphism of the ring,`,
` sigma(.,k) is its kth iterate. `,
` . delta: a sigma-derivation, which means that `,
` delta(f g) = sigma(f) delta(g) + delta(f) g`,
` . x: the variable.`,` `,
`- In the differential case, sigma should be identity and delta d/dx.`,
`- In the difference case, sigma should be shift and delta(f)=shift(f)-f.`,
`- In the shift case, sigma should be shift and delta=0.`,
`- In the q-difference case, sigma should be Q:x->q*x and delta(f)=Q(f)-f.`,
` `,`- The functions available in this package are:`,` `,
` degree\tldeg\t\t+\t\t*\t\tprint`,
` rightquo\trightrem\trightgcd\trightgcdex\trightlcm`,
` leftquo\tleftrem\t\tleftgcd\t\tleftgcdex\tleftlcm`,
` apply\trightquorem\tleftquorem\trightpseudoquorem`,
` rightannihilators\t\tleftannihilators`,` `,
`- Most of these functions are the analogous of the usual Maple functions on`,
`commutative polynomials. However, non-commutativity makes it necessary to`,
`distinguish between left and right quotients, remainders, gcds and lcms.`,
`The function apply applies the operator to an element of the ground field.`,`\
The function print takes as argument an operator and a symbol for the derivati\
on`,
`variable. It displays the operator in "natural" form. The functions right and`
,`left annihilators compute two Ore polynomials U and V such that UA+VB=0 or`,
`AU+VB=0. `,` `,
`- Besides the Ore polynomial(s), these functions take the field as last `,
`argument. More detailed help on each of these functions has not (yet) been`,
`written. `,` `,
`- The Ore polynomials are represented internally as a list of coefficients`,
`of the differentiation. Currently there is no friendly user interface between`
,`this representation and (say) Maple recurrences or differential equations.`,
` `,`- If you have any comment or problem, please send me mail at`,
`Bruno.Salvy@inria.fr`,` `,`EXAMPLE: `,
`> DifferentialOperators:=[(f,k)->f,f->diff(f,x),x]:`,
`> A:=[x,2/x,1]:B:=[x,x^3,x^2+1]:C:=[x+3,x^2+x+1]:`,
`> F:=``OrePoly/*``(A,B,DifferentialOperators):`,
`> G:=``OrePoly/+``(F,C,DifferentialOperators):`,
`> ``OrePoly/print``(G,D,DifferentialOperators);`,` 2 4 \
2 2 2 2`,
` (x + 2/x + x + 3) + (x + 13 x + 5 + x ) D + (x (x + 1) + 6 + 8 x + x) D`
,` `,` / 2 \\`,
` | x + 1 3| 3 2 4`,
` + |2 ------ + 4 x + x | D + (x + 1) D`,
` \\ x /`,` `,
`> ``OrePoly/print``(``OrePoly/rightrem``(G,B,DifferentialOperators),D);`,
` 2`,
` (x + 3) + (x + x + 1) D`):
#save `OrePoly.m`;
#quit