#   -*-Maple-*-
########################## Ore polynomials #################################
##
##    Title:    OrePoly
##    Created:  Sept 1993
##    Author:   Bruno Salvy
##              <Bruno.Salvy@inria.fr>
##
## 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<da then [[0],f]
    else
        lc:=a[nops(a)];
        rem:=f;
        for i from df by -1 to da do
            if nops(rem)=i+1 then
                if right then
                    q[i-da]:=rem[i+1]/F[1](lc,i-da);
                    rem:=`OrePoly/+`(rem,
                        [seq(-j,j=`OrePoly/*`([0$(i-da),q[i-da]],a,F))])
                else
                    q[i-da]:=F[1](rem[i+1]/lc,-da);
                    rem:=`OrePoly/+`(rem,
                        [seq(-j,j=`OrePoly/*`(a,[0$(i-da),q[i-da]],F))])
                fi
            else q[i-da]:=0 fi
        od;
        [[seq(q[i],i=0..df-da)],rem]
    fi
end: # `OrePoly/rightleftquorem`

`OrePoly/rightquorem`:=proc(f,g,F) `OrePoly/rightleftquorem`(f,g,F,true) end:
`OrePoly/leftquorem`:=proc(f,g,F) `OrePoly/rightleftquorem`(f,g,F,false) end:

# OrePoly/rightquo, OrePoly/leftquo
#Input: both polynomials and the field
#Output: the right (resp. right) quotient
`OrePoly/rightquo`:=proc() op(1,`OrePoly/rightquorem`(args)) end:
`OrePoly/leftquo`:=proc() op(1,`OrePoly/leftquorem`(args)) end:

# OrePoly/rightrem, OrePoly/leftrem
#Input: both polynomials and the field
#Output: the right (resp. left) remainder
# Avoid storing the coefficients of the quotient,
# otherwise, same as rightquorem, leftquorem. 
# This difference justifies the duplication of code in very large computations.
`OrePoly/rightleftrem`:=proc (f,a,F,right)
local df, da, i, rem, j, lc;
    df:=`OrePoly/degree`(f); da:=`OrePoly/degree`(a);
    if df<da then f
    elif da=0 then 0
    else
        lc:=a[nops(a)];
        rem:=f;
        for i from df by -1 to da do if nops(rem)=i+1 then
            if right then
                rem:=`OrePoly/+`(rem,[seq(-j,j=`OrePoly/*`([0$(i-da),
                    rem[nops(rem)]/F[1](lc,i-da)],a,F))])
            else
                rem:=`OrePoly/+`(rem,[seq(-j,j=`OrePoly/*`(a,
                    [0$(i-da),F[1](rem[i+1]/lc,-da)],F))])
            fi fi od;
        rem
    fi
end: # `OrePoly/rightrem`

`OrePoly/rightrem`:=proc(f,a,F) `OrePoly/rightleftrem`(f,a,F,true) end:
`OrePoly/leftrem`:=proc(f,a,F) `OrePoly/rightleftrem`(f,a,F,false) end:

# OrePoly/leftgcd, OrePoly/rightgcd
#Input: both polynomials and the field
#Output: the left (right) gcd
`OrePoly/rightgcd`:=proc(f,g,F)
local da,db,a,b,q,dq;
    da:=`OrePoly/degree`(f); db:=`OrePoly/degree`(g);
    if da<db then a:=g; b:=f; q:=da; da:=db; db:=q 
    else a:=f; b:=g fi;
    do
        if `OrePoly/degree`(b)=0 then
            if b<>[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<db then       
        a:=g;b:=f;q:=da;da:=db;db:=q;u0:=[0];v0:=[1];u1:=[1];v1:=[0] 
    else a:=f;b:=g;u0:=[1];v0:=[0];u1:=[0];v1:=[1] fi;
    if b=[0] then u:=u1;v:=v1; RETURN(a) fi;
    do
        if db=0 then
            if b=[0] and not lcm then u:=u0; v:=v0; RETURN(a)
            elif b=[0] or not lcm then u:=u1; v:=v1; RETURN(b) fi fi;
        qr:=`OrePoly/rightleftquorem`(a,b,F,right);
        dq:=`OrePoly/degree`(qr[2]);
        da:=db;db:=dq;a:=b;b:=qr[2];
        t:=u0;u0:=u1;
        if right then u1:=`OrePoly/+`(t,[seq(-i,i=`OrePoly/*`(qr[1],u1,F))])
        else          u1:=`OrePoly/+`(t,[seq(-i,i=`OrePoly/*`(u1,qr[1],F))])fi;
        t:=v0;v0:=v1;
        if right then v1:=`OrePoly/+`(t,[seq(-i,i=`OrePoly/*`(qr[1],v1,F))])
        else          v1:=`OrePoly/+`(t,[seq(-i,i=`OrePoly/*`(v1,qr[1],F))]) fi
    od
end: # `OrePoly/rightleftgcdex`

`OrePoly/leftgcdex` :=proc() `OrePoly/rightleftgcdex`(args,false,false) end:
`OrePoly/rightgcdex`:=proc() `OrePoly/rightleftgcdex`(args,true, false) end:

# OrePoly/leftlcm, OrePoly/rightlcm
#Input: both polynomials and the field
#Output: their left (right) least common multiple
`OrePoly/leftlcm`:=proc(f,g,F)
local u,v;
    `OrePoly/rightleftgcdex`(f,g,'u','v',F,true,true);
    `OrePoly/*`(u,f,F)
end: # `OrePoly/leftlcm`

`OrePoly/rightlcm`:=proc(f,g,F)
local u,v;
    `OrePoly/rightleftgcdex`(f,g,'u','v',F,false,true);
    `OrePoly/*`(f,u,F)
end: # `OrePoly/leftlcm`

# OrePoly/rightannihilators, OrePoly/leftannihilators
#Input: two polynomials A and B, the field,
#       an optional 4th argument tells that the result should be normalized.
#Output: a pair of two polynomials U and V such that
#  AU + BV = 0 (UA + BV = 0)
`OrePoly/rightleftannihilators`:=proc(A,B,F,right,normalized)
local u,v,s,i;
        `OrePoly/rightleftgcdex`(A,B,u,v,F,true,not right);
        if nargs=5 then
            s:=lcm(op(map(denom,u)),op(map(denom,v)));
            u:=map(normal,[seq(i*s,i=u)]);
            v:=map(normal,[seq(i*s,i=v)])
        fi;
        [u,v]
end:
`OrePoly/rightannihilators`:=proc(A,B,F)
    `OrePoly/rightleftannihilators`(A,B,F,true,args[4..nargs]) end:
`OrePoly/leftannihilators`:=proc(A,B,F)
    `OrePoly/rightleftannihilators`(A,B,F,false,args[4..nargs]) end:

# OrePoly/rightpseudoquorem
#Input: both polynomials F and Q, the field plus a name C
#Output: the right pseudo-quotient and rem, such that
#  CF = QA+R
`OrePoly/rightpseudoquorem`:=proc(f,a,F,C)
local df,da,q,i,rem,j,lc,cof;
    df:=`OrePoly/degree`(f); da:=`OrePoly/degree`(a);
    if df<da then [[0],f]
    else
        lc[0]:=a[nops(a)];
        for i to df-da do lc[i]:=F[1](lc[i-1],i) od;
        cof:=lcm(seq(lc[i],i=0..df-da));
        rem:=`OrePoly/*`([cof],f,F);
        for i from df by -1 to da do
            if nops(rem)=i+1 then
                q[i-da]:=quo(rem[nops(rem)],lc[i-da],F[3]);
                rem:=`OrePoly/+`(rem,
                    [seq(-j,j=`OrePoly/*`([0$(i-da),q[i-da]],a,F))])
            else q[i-da]:=0 fi
        od;
        C:=cof;
        [[seq(q[i],i=0..df-da)],rem]
    fi
end: # `OrePoly/pseudorightquorem`

# OrePoly/apply
#Input: a polynomial, a constant c, an element of something on which 
# the polynomial operates, the field.
#Output: the result, knowing that the action of Delta is
#   Delta(f) = c sigma(f,1) + delta(f).
# This constant c should be 0 for differential and difference operators,
#                           1 for shift operators.
`OrePoly/apply`:=proc(f,c,a,F) local i, res;
    res[1]:=a;
    for i from 2 to nops(f) do res[i]:=c*F[1](res[i-1],1)+F[2](res[i-1]) od;
    convert([seq(f[i]*res[i],i=1..nops(f))],`+`)
end: # `OrePoly/apply`

############## HELP FILE GENERATED BY makehelp FROM OrePoly.help #######
`help/text/OrePoly` := TEXT(`HELP FOR: The OrePoly functions`,`   `,
`CALLING SEQUENCE:`,`   ``OrePoly/<function>``(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