--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--All rights reserved.
--
--Redistribution and use in source and binary forms, with or without
--modification, are permitted provided that the following conditions are
--met:
--
--    - Redistributions of source code must retain the above copyright
--      notice, this list of conditions and the following disclaimer.
--
--    - Redistributions in binary form must reproduce the above copyright
--      notice, this list of conditions and the following disclaimer in
--      the documentation and/or other materials provided with the
--      distribution.
--
--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
--      names of its contributors may be used to endorse or promote products
--      derived from this software without specific prior written permission.
--
--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 
)abbrev category OREPCAT UnivariateSkewPolynomialCategory
++ Author: Manuel Bronstein, Jean Della Dora, Stephen M. Watt
++ Date Created: 19 October 1993
++ Date Last Updated: 1 February 1994
++ Description:
++   This is the category of univariate skew polynomials over an Ore
++   coefficient ring.
++   The multiplication is given by \spad{x a = \sigma(a) x + \delta a}.
++   This category is an evolution of the types
++     MonogenicLinearOperator, OppositeMonogenicLinearOperator, and
++     NonCommutativeOperatorDivision
++   developped by Jean Della Dora and Stephen M. Watt.
UnivariateSkewPolynomialCategory(R:Ring):
  Category == Join(Ring, BiModule(R, R), FullyRetractableTo R) with
        degree: $ -> NonNegativeInteger
            ++ degree(l) is \spad{n} if
            ++   \spad{l = sum(monomial(a(i),i), i = 0..n)}.
        minimumDegree: $ -> NonNegativeInteger
            ++ minimumDegree(l) is the smallest \spad{k} such that
            ++ \spad{a(k) ~= 0} if
            ++   \spad{l = sum(monomial(a(i),i), i = 0..n)}.
        leadingCoefficient: $ -> R
            ++ leadingCoefficient(l) is \spad{a(n)} if
            ++   \spad{l = sum(monomial(a(i),i), i = 0..n)}.
        reductum: $ -> $
            ++ reductum(l) is \spad{l - monomial(a(n),n)} if
            ++   \spad{l = sum(monomial(a(i),i), i = 0..n)}.
        coefficient: ($, NonNegativeInteger) -> R
            ++ coefficient(l,k) is \spad{a(k)} if
            ++   \spad{l = sum(monomial(a(i),i), i = 0..n)}.
        monomial: (R, NonNegativeInteger) -> $
            ++ monomial(c,k) produces c times the k-th power of
            ++ the generating operator, \spad{monomial(1,1)}.
        coefficients: % -> List R
            ++ coefficients(l) returns the list of all the nonzero
            ++ coefficients of l.
        apply: (%, R, R) -> R
          ++ apply(p, c, m) returns \spad{p(m)} where the action is
          ++ given by \spad{x m = c sigma(m) + delta(m)}.
        if R has CommutativeRing then Algebra R
        if R has IntegralDomain then
          exquo: (%, R) -> Union(%, "failed")
            ++ exquo(l, a) returns the exact quotient of l by a, 
            ++ returning \axiom{"failed"} if this is not possible.
          monicLeftDivide:   (%, %) -> Record(quotient: %, remainder: %)
            ++ monicLeftDivide(a,b) returns the pair \spad{[q,r]} such that
            ++ \spad{a = b*q + r} and the degree of \spad{r} is
            ++ less than the degree of \spad{b}.
            ++ \spad{b} must be monic.
            ++ This process is called ``left division''.
          monicRightDivide:   (%, %) -> Record(quotient: %, remainder: %)
            ++ monicRightDivide(a,b) returns the pair \spad{[q,r]} such that
            ++ \spad{a = q*b + r} and the degree of \spad{r} is
            ++ less than the degree of \spad{b}.
            ++ \spad{b} must be monic.
            ++ This process is called ``right division''.
        if R has GcdDomain then
          content: % -> R
            ++ content(l) returns the gcd of all the coefficients of l.
          primitivePart: % -> %
            ++ primitivePart(l) returns l0 such that \spad{l = a * l0}
            ++ for some a in R, and \spad{content(l0) = 1}.
        if R has Field then
          leftDivide:   (%, %) -> Record(quotient: %, remainder: %)
            ++ leftDivide(a,b) returns the pair \spad{[q,r]} such that
            ++ \spad{a = b*q + r} and the degree of \spad{r} is
            ++ less than the degree of \spad{b}.
            ++ This process is called ``left division''.
          leftQuotient:  (%, %) -> %
            ++ leftQuotient(a,b) computes the pair \spad{[q,r]} such that
            ++ \spad{a = b*q + r} and the degree of \spad{r} is
            ++ less than the degree of \spad{b}.
            ++ The value \spad{q} is returned.
          leftRemainder:  (%, %) -> %
            ++ leftRemainder(a,b) computes the pair \spad{[q,r]} such that
            ++ \spad{a = b*q + r} and the degree of \spad{r} is
            ++ less than the degree of \spad{b}.
            ++ The value \spad{r} is returned.
          leftExactQuotient:(%, %) -> Union(%, "failed")
            ++ leftExactQuotient(a,b) computes the value \spad{q}, if it exists,
            ++  such that \spad{a = b*q}.
          leftGcd:   (%, %) -> %
            ++ leftGcd(a,b) computes the value \spad{g} of highest degree
            ++ such that
            ++    \spad{a = g*aa}
            ++    \spad{b = g*bb}
            ++ for some values \spad{aa} and \spad{bb}.
            ++ The value \spad{g} is computed using left-division.
          leftExtendedGcd:   (%, %) -> Record(coef1:%, coef2:%, generator:%)
            ++ leftExtendedGcd(a,b) returns \spad{[c,d]} such that
            ++ \spad{g = a * c + b * d = leftGcd(a, b)}.
          rightLcm:   (%, %) -> %
            ++ rightLcm(a,b) computes the value \spad{m} of lowest degree
            ++ such that \spad{m = a*aa = b*bb} for some values
            ++ \spad{aa} and \spad{bb}.  The value \spad{m} is
            ++ computed using left-division.
          rightDivide:   (%, %) -> Record(quotient: %, remainder: %)
            ++ rightDivide(a,b) returns the pair \spad{[q,r]} such that
            ++ \spad{a = q*b + r} and the degree of \spad{r} is
            ++ less than the degree of \spad{b}.
            ++ This process is called ``right division''.
          rightQuotient:  (%, %) -> %
            ++ rightQuotient(a,b) computes the pair \spad{[q,r]} such that
            ++ \spad{a = q*b + r} and the degree of \spad{r} is
            ++ less than the degree of \spad{b}.
            ++ The value \spad{q} is returned.
          rightRemainder:  (%, %) -> %
            ++ rightRemainder(a,b) computes the pair \spad{[q,r]} such that
            ++ \spad{a = q*b + r} and the degree of \spad{r} is
            ++ less than the degree of \spad{b}.
            ++ The value \spad{r} is returned.
          rightExactQuotient:(%, %) -> Union(%, "failed")
            ++ rightExactQuotient(a,b) computes the value \spad{q}, if it exists
            ++ such that \spad{a = q*b}.
          rightGcd:   (%, %) -> %
            ++ rightGcd(a,b) computes the value \spad{g} of highest degree
            ++ such that
            ++    \spad{a = aa*g}
            ++    \spad{b = bb*g}
            ++ for some values \spad{aa} and \spad{bb}.
            ++ The value \spad{g} is computed using right-division.
          rightExtendedGcd:   (%, %) -> Record(coef1:%, coef2:%, generator:%)
            ++ rightExtendedGcd(a,b) returns \spad{[c,d]} such that
            ++ \spad{g = c * a + d * b = rightGcd(a, b)}.
          leftLcm:   (%, %) -> %
            ++ leftLcm(a,b) computes the value \spad{m} of lowest degree
            ++ such that \spad{m = aa*a = bb*b} for some values
            ++ \spad{aa} and \spad{bb}.  The value \spad{m} is
            ++ computed using right-division.
 
   add
        coerce(x:R):% == monomial(x, 0)
 
        coefficients l ==
          ans:List(R) := empty()
          while l ~= 0 repeat
            ans := concat(leadingCoefficient l, ans)
            l   := reductum l
          ans
 
        a:R * y:% ==
          z:% := 0
          while y ~= 0 repeat
            z := z + monomial(a * leadingCoefficient y, degree y)
            y := reductum y
          z
 
        retractIfCan(x:%):Union(R, "failed") ==
          zero? x or zero? degree x => leadingCoefficient x
          "failed"
 
        if R has IntegralDomain then
          l exquo a ==
            ans:% := 0
            while l ~= 0 repeat
              (u := (leadingCoefficient(l) exquo a)) case "failed" =>
                 return "failed"
              ans := ans + monomial(u::R, degree l)
              l   := reductum l
            ans
 
        if R has GcdDomain then
          content l       == gcd coefficients l
          primitivePart l == (l exquo content l)::%
 
        if R has Field then
          leftEEA:  (%, %) -> Record(gcd:%, coef1:%, coef2:%, lcm:%)
          rightEEA: (%, %) -> Record(gcd:%, coef1:%, coef2:%, lcm:%)
          ncgcd:    (%, %, (%, %) -> %) -> %
          nclcm:  (%, %, (%, %) -> Record(gcd:%, coef1:%, coef2:%, lcm:%)) -> %
          exactQuotient: Record(quotient:%, remainder:%) -> Union(%, "failed")
          extended: (%, %, (%, %) -> Record(gcd:%, coef1:%, coef2:%, lcm:%)) ->
                                          Record(coef1:%, coef2:%, generator:%)
 
          leftQuotient(a, b)         == leftDivide(a,b).quotient
          leftRemainder(a, b)        == leftDivide(a,b).remainder
          leftExtendedGcd(a, b)      == extended(a, b, leftEEA)
          rightLcm(a, b)             == nclcm(a, b, leftEEA)
          rightQuotient(a, b)        == rightDivide(a,b).quotient
          rightRemainder(a, b)       == rightDivide(a,b).remainder
          rightExtendedGcd(a, b)     == extended(a, b, rightEEA)
          leftLcm(a, b)              == nclcm(a, b, rightEEA)
          leftExactQuotient(a, b)    == exactQuotient leftDivide(a, b)
          rightExactQuotient(a, b)   == exactQuotient rightDivide(a, b)
          rightGcd(a, b)             == ncgcd(a, b, rightRemainder)
          leftGcd(a, b)              == ncgcd(a, b, leftRemainder)
          exactQuotient qr  == (zero?(qr.remainder) => qr.quotient; "failed")
 
          -- returns [g = leftGcd(a, b), c, d, l = rightLcm(a, b)]
          -- such that g := a c + b d
          leftEEA(a, b) ==
            a0 := a
            u0:% := v:% := 1
            v0:% := u:% := 0
            while b ~= 0 repeat
              qr     := leftDivide(a, b)
              (a, b) := (b, qr.remainder)
              (u0, u):= (u, u0 - u * qr.quotient)
              (v0, v):= (v, v0 - v * qr.quotient)
            [a, u0, v0, a0 * u]
 
          ncgcd(a, b, ncrem) ==
            zero? a => b
            zero? b => a
            degree a < degree b => ncgcd(b, a, ncrem)
            while b ~= 0 repeat (a, b) := (b, ncrem(a, b))
            a
 
          extended(a, b, eea) ==
            zero? a => [0, 1, b]
            zero? b => [1, 0, a]
            degree a < degree b =>
              rec := eea(b, a)
              [rec.coef2, rec.coef1, rec.gcd]
            rec := eea(a, b)
            [rec.coef1, rec.coef2, rec.gcd]
 
          nclcm(a, b, eea) ==
            zero? a or zero? b => 0
            degree a < degree b => nclcm(b, a, eea)
            rec := eea(a, b)
            rec.lcm
 
          -- returns [g = rightGcd(a, b), c, d, l = leftLcm(a, b)]
          -- such that g := a c + b d
          rightEEA(a, b) ==
            a0 := a
            u0:% := v:% := 1
            v0:% := u:% := 0
            while b ~= 0 repeat
              qr     := rightDivide(a, b)
              (a, b) := (b, qr.remainder)
              (u0, u):= (u, u0 - qr.quotient * u)
              (v0, v):= (v, v0 - qr.quotient * v)
            [a, u0, v0, u * a0]

)abbrev package APPLYORE ApplyUnivariateSkewPolynomial
++ Author: Manuel Bronstein
++ Date Created: 7 December 1993
++ Date Last Updated: 1 February 1994
++ Description:
++   \spad{ApplyUnivariateSkewPolynomial} (internal) allows univariate
++   skew polynomials to be applied to appropriate modules.
ApplyUnivariateSkewPolynomial(R:Ring, M: LeftModule R,
    P: UnivariateSkewPolynomialCategory R): with
      apply: (P, M -> M, M) -> M
        ++ apply(p, f, m) returns \spad{p(m)} where the action is given
        ++ by \spad{x m = f(m)}.
        ++ \spad{f} must be an R-pseudo linear map on M.
   == add
      apply(p, f, m) ==
        w:M  := 0
        mn:M := m
        for i in 0..degree p repeat
          w  := w + coefficient(p, i) * mn
          mn := f mn
        w

import Integer
import NonNegativeInteger
)abbrev domain AUTOMOR Automorphism
++ Author: Manuel Bronstein
++ Date Created: 31 January 1994
++ Date Last Updated: 31 January 1994
++ References:
++ Description:
++       Automorphism R is the multiplicative group of automorphisms of R.
-- In fact, non-invertible endomorphism are allowed as partial functions.
-- This domain is noncanonical in that f*f^{-1} will be the identity
-- function but won't be equal to 1.
Automorphism(R:Ring): Join(Group, Eltable(R, R)) with
      morphism: (R -> R) -> %
        ++ morphism(f) returns the non-invertible morphism given by f.
      morphism: (R -> R, R -> R) -> %
        ++ morphism(f, g) returns the invertible morphism given by f, where
        ++ g is the inverse of f..
      morphism: ((R, Integer) -> R) -> %
        ++ morphism(f) returns the morphism given by \spad{f^n(x) = f(x,n)}.
   == add
      Rep == (R, Integer) -> R
      err:   R -> R
      ident: (R, Integer) -> R
      iter:  (R -> R, NonNegativeInteger, R) -> R
      iterat: (R -> R, R -> R, Integer, R) -> R
      apply: (Rep, R, Integer) -> R
  
      1                               == per ident
      err r                           == error "Morphism is not invertible"
      ident(r, n)                     == r
      f = g                           == %peq(f,g)$Foreign(Builtin)
      elt(f, r)                       == apply(rep f, r, 1)
      inv f                           == per apply(rep f, #1, - #2)
      (f: %) ** (n: Integer)          == per apply(rep f, #1, n * #2)
      coerce(f:%):OutputForm          == message("R -> R")
      morphism(f:(R, Integer) -> R):% == per f
      morphism(f:R -> R):%            == morphism(f, err)
      morphism(f, g)                  == per iterat(f, g, #2, #1)
      apply(f, r, n) == f(r, n)
 
      iterat(f, g, n, r) ==
          negative? n => iter(g, (-n)::NonNegativeInteger, r)
          iter(f, n::NonNegativeInteger, r)
 
      iter(f, n, r) ==
          for i in 1..n repeat r := f r
          r
 
      f * g ==
        f = g => f**2
        per iterat(f g #1, (inv g)(inv f) #1, #2, #1)

)abbrev package OREPCTO UnivariateSkewPolynomialCategoryOps
++ Author: Manuel Bronstein
++ Date Created: 1 February 1994
++ Date Last Updated: 1 February 1994
++ Description:
++   \spad{UnivariateSkewPolynomialCategoryOps} provides products and
++    divisions of univariate skew polynomials.
-- Putting those operations here rather than defaults in OREPCAT allows
-- OREPCAT to be defined independently of sigma and delta.
-- MB 2/94
UnivariateSkewPolynomialCategoryOps(R, C): Exports == Implementation where
    R: Ring
    C: UnivariateSkewPolynomialCategory R
 
    N   ==> NonNegativeInteger
    MOR ==> Automorphism R
    QUOREM ==> Record(quotient: C, remainder: C)
 
    Exports ==> with
        times: (C, C, MOR, R -> R) -> C
           ++ times(p, q, sigma, delta) returns \spad{p * q}.
           ++ \spad{\sigma} and \spad{\delta} are the maps to use.
        apply: (C, R, R, MOR, R -> R) -> R
          ++ apply(p, c, m, sigma, delta) returns \spad{p(m)} where the action
          ++ is given by \spad{x m = c sigma(m) + delta(m)}.
        if R has IntegralDomain then
            monicLeftDivide: (C, C, MOR) -> QUOREM
                ++ monicLeftDivide(a, b, sigma) returns the pair \spad{[q,r]}
                ++ such that \spad{a = b*q + r} and the degree of \spad{r} is
                ++ less than the degree of \spad{b}.
                ++ \spad{b} must be monic.
                ++ This process is called ``left division''.
                ++ \spad{\sigma} is the morphism to use.
            monicRightDivide: (C, C, MOR) -> QUOREM
                ++ monicRightDivide(a, b, sigma) returns the pair \spad{[q,r]}
                ++ such that \spad{a = q*b + r} and the degree of \spad{r} is
                ++ less than the degree of \spad{b}.
                ++ \spad{b} must be monic.
                ++ This process is called ``right division''.
                ++ \spad{\sigma} is the morphism to use.
        if R has Field then
            leftDivide: (C, C, MOR) -> QUOREM
                ++ leftDivide(a, b, sigma) returns the pair \spad{[q,r]} such
                ++ that \spad{a = b*q + r} and the degree of \spad{r} is
                ++ less than the degree of \spad{b}.
                ++ This process is called ``left division''.
                ++ \spad{\sigma} is the morphism to use.
            rightDivide: (C, C, MOR) -> QUOREM
                ++ rightDivide(a, b, sigma) returns the pair \spad{[q,r]} such
                ++ that \spad{a = q*b + r} and the degree of \spad{r} is
                ++ less than the degree of \spad{b}.
                ++ This process is called ``right division''.
                ++ \spad{\sigma} is the morphism to use.
 
    Implementation ==> add
        termPoly:         (R, N, C, MOR, R -> R) -> C
        localLeftDivide : (C, C, MOR, R) -> QUOREM
        localRightDivide: (C, C, MOR, R) -> QUOREM
 
        times(x, y, sigma, delta) ==
          zero? y => 0
          z:C := 0
          while x ~= 0 repeat
            z := z + termPoly(leadingCoefficient x, degree x, y, sigma, delta)
            x := reductum x
          z
 
        termPoly(a, n, y, sigma, delta) ==
          zero? y => 0
          (u := subtractIfCan(n, 1)) case "failed" => a * y
          n1 := u::N
          z:C := 0
          while y ~= 0 repeat
            m := degree y
            b := leadingCoefficient y
            z := z + termPoly(a, n1, monomial(sigma b, m + 1), sigma, delta)
                   + termPoly(a, n1, monomial(delta b, m), sigma, delta)
            y := reductum y
          z
 
        apply(p, c, x, sigma, delta) ==
          w:R  := 0
          xn:R := x
          for i in 0..degree p repeat
            w  := w + coefficient(p, i) * xn
            xn := c * sigma xn + delta xn
          w
 
        -- localLeftDivide(a, b) returns [q, r] such that a = q b + r
        -- b1 is the inverse of the leadingCoefficient of b
        localLeftDivide(a, b, sigma, b1) ==
            zero? b => error "leftDivide: division by 0"
            zero? a or
             (n := subtractIfCan(degree(a),(m := degree b))) case "failed" =>
                    [0,a]
            q  := monomial((sigma**(-m))(b1 * leadingCoefficient a), n::N)
            qr := localLeftDivide(a - b * q, b, sigma, b1)
            [q + qr.quotient, qr.remainder]
 
        -- localRightDivide(a, b) returns [q, r] such that a = q b + r
        -- b1 is the inverse of the leadingCoefficient of b
        localRightDivide(a, b, sigma, b1) ==
            zero? b => error "rightDivide: division by 0"
            zero? a or
              (n := subtractIfCan(degree(a),(m := degree b))) case "failed" =>
                    [0,a]
            q := monomial(leadingCoefficient(a) * (sigma**n) b1, n::N)
            qr := localRightDivide(a - q * b, b, sigma, b1)
            [q + qr.quotient, qr.remainder]
 
        if R has IntegralDomain then
            monicLeftDivide(a, b, sigma) ==
                unit?(u := leadingCoefficient b) =>
                    localLeftDivide(a, b, sigma, recip(u)::R)
                error "monicLeftDivide: divisor is not monic"
 
            monicRightDivide(a, b, sigma) ==
                unit?(u := leadingCoefficient b) =>
                    localRightDivide(a, b, sigma, recip(u)::R)
                error "monicRightDivide: divisor is not monic"
 
        if R has Field then
            leftDivide(a, b, sigma) ==
                localLeftDivide(a, b, sigma, inv leadingCoefficient b)
 
            rightDivide(a, b, sigma) ==
                localRightDivide(a, b, sigma, inv leadingCoefficient b)

)abbrev domain ORESUP SparseUnivariateSkewPolynomial
++ Author: Manuel Bronstein
++ Date Created: 19 October 1993
++ Date Last Updated: September, 2008
++ Description:
++   This is the domain of sparse univariate skew polynomials over an Ore
++   coefficient field.
++   The multiplication is given by \spad{x a = \sigma(a) x + \delta a}.
SparseUnivariateSkewPolynomial(R:Ring, sigma:Automorphism R, delta: R -> R):
 UnivariateSkewPolynomialCategory R with
      outputForm: (%, OutputForm) -> OutputForm
        ++ outputForm(p, x) returns the output form of p using x for the
        ++ otherwise anonymous variable.
   == SparseUnivariatePolynomial R add
      import UnivariateSkewPolynomialCategoryOps(R, %)
 
      x:% * y:%      == times(x, y, sigma, delta)
      apply(p, c, r) == apply(p, c, r, sigma, delta)
      x:% ** n:PositiveInteger == expt(x,n)$RepeatedSquaring(%)
      x:% ** n:NonNegativeInteger ==
        zero? n => 1
        expt(x,n::PositiveInteger)$RepeatedSquaring(%)
 
      if R has IntegralDomain then
          monicLeftDivide(a, b)  == monicLeftDivide(a, b, sigma)
          monicRightDivide(a, b) == monicRightDivide(a, b, sigma)
 
      if R has Field then
          leftDivide(a, b)  == leftDivide(a, b, sigma)
          rightDivide(a, b) == rightDivide(a, b, sigma)

)abbrev domain OREUP UnivariateSkewPolynomial
++ Author: Manuel Bronstein
++ Date Created: 19 October 1993
++ Date Last Updated: 1 February 1994
++ Description:
++   This is the domain of univariate skew polynomials over an Ore
++   coefficient field in a named variable.
++   The multiplication is given by \spad{x a = \sigma(a) x + \delta a}.
UnivariateSkewPolynomial(x:Symbol, R:Ring, sigma:Automorphism R, delta: R -> R):
 Join(UnivariateSkewPolynomialCategory R,CoercibleFrom Variable x)
  == SparseUnivariateSkewPolynomial(R, sigma, delta) add
     Rep := SparseUnivariateSkewPolynomial(R, sigma, delta)
     coerce(v:Variable(x)):% == monomial(1, 1)
     coerce(p:%):OutputForm  == outputForm(p, outputForm x)$Rep

