--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.

import OrderedSet
import OrderedMonoid
import RetractableTo
)abbrev domain OFMONOID OrderedFreeMonoid
++ Author: Michel Petitot petitot@lifl.fr
++ Date Created: 91
++ Date Last Updated: 7 Juillet 92
++ Fix History: compilation v 2.1 le 13 dec 98
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++    The free monoid on a set \spad{S} is the monoid of finite products of
++ the form \spad{reduce(*,[si ** ni])} where the si's are in S, and the ni's
++ are non-negative integers. The multiplication is not commutative.
++ For two elements \spad{x} and \spad{y} the relation \spad{x < y}
++ holds if either \spad{length(x) < length(y)} holds or if these lengths
++ are equal and if \spad{x} is smaller than \spad{y} w.r.t. the lexicographical
++ ordering induced by \spad{S}.
++ This domain inherits implementation from \spadtype{FreeMonoid}.
++ Author: Michel Petitot (petitot@lifl.fr)

OrderedFreeMonoid(S: OrderedSet): OFMcategory == OFMdefinition where
    NNI ==> NonNegativeInteger
    REC ==> Record(gen:S, exp:NNI)
 
    OFMcategory == Join(FreeMonoidCategory S,OrderedMonoid) with
        first: % -> S
          ++ \spad{first(x)} returns the first letter of \spad{x}.
        rest:  % -> %
          ++ \spad{rest(x)} returns \spad{x} except the first letter.
        mirror: % -> %
          ++ \spad{mirror(x)} returns the reversed word of \spad{x}.
        lexico: (%,%) -> Boolean
          ++ \spad{lexico(x,y)} returns \spad{true} iff \spad{x} is smaller than \spad{y}
          ++ w.r.t. the pure lexicographical ordering induced by \spad{S}.
        lquo:   (%, S) -> Union(%, "failed")
          ++ \spad{lquo(x, s)} returns the exact left quotient of \spad{x} 
          ++ by \spad{s}. 
        rquo:   (%, S) -> Union(%, "failed")
          ++ \spad{rquo(x, s)} returns the exact right quotient 
          ++ of \spad{x} by \spad{s}.
        div:   (%, %) -> Union(Record(lm: %, rm: %), "failed")
          ++ \spad{x div y} returns the left and right exact quotients of
          ++ \spad{x} by \spad{y}, that is \spad{[l, r]} such that \spad{x = l * y * r}.
          ++ "failed" is returned iff \spad{x} is not of the form \spad{l * y * r}.
          ++ monomial of \spad{x}.
        length: % -> NNI
          ++ \spad{length(x)} returns the length of \spad{x}.
        varList: % -> List S
          ++ \spad{varList(x)} returns the list of variables of \spad{x}.

    OFMdefinition == FreeMonoid(S) add
        Rep := ListMonoidOps(S, NNI, 1)
        
      -- definitions
        lquo(w:%, l:S) == 
          x: List REC := listOfMonoms(w)$Rep
          null x        => "failed"
          fx: REC := first x
          fx.gen ~= l  => "failed"
          fx.exp = 1   => makeMulti rest(x)
          makeMulti [[fx.gen, (fx.exp - 1)::NNI ]$REC, :rest x]
       
        rquo(w:%, l:S) ==
          u:% := reverse w
          (r := lquo (u,l)) case "failed" => "failed"
          reverse! (r::%)

        length x == reduce("+" ,[f.exp for f in listOfMonoms x], 0)

        varList x ==
          le: List S := [t.gen for t in listOfMonoms x]
          sort! removeDuplicates(le)
 
        first w ==
          x: List REC := listOfMonoms w
          null x => error "empty word !!!"
          x.first.gen

        rest w ==
          x: List REC := listOfMonoms w
          null x => error "empty word !!!"
          fx: REC := first x
          fx.exp = 1 => makeMulti rest x
          makeMulti [[fx.gen , (fx.exp - 1)::NNI ]$REC , :rest x]

        lexico(a,b) ==         --  ordre lexicographique
            la := listOfMonoms a
            lb := listOfMonoms b
            while (not null la) and (not null lb) repeat
                la.first.gen > lb.first.gen => return false
                la.first.gen < lb.first.gen => return true
                if la.first.exp = lb.first.exp then
                    la:=rest la
                    lb:=rest lb
                else if la.first.exp > lb.first.exp then
                    la:=concat([la.first.gen,
                           (la.first.exp - lb.first.exp)::NNI], rest lb)
                    lb:=rest lb
                else
                    lb:=concat([lb.first.gen,
                             (lb.first.exp-la.first.exp)::NNI], rest la)
                    la:=rest la
            empty? la and not empty? lb


        a < b ==               --  ordre lexicographique par longueur
            la:NNI := length a; lb:NNI := length b
            la = lb =>  lexico(a,b)
            la < lb 

        mirror x == reverse(x)$Rep

import Ring
import SetCategory
import BiModule
import RetractableTo
)abbrev category FMCAT FreeModuleCat
++ Author: Michel Petitot petitot@lifl.fr
++ Date Created: 91
++ Date Last Updated: 7 Juillet 92
++ Fix History: compilation v 2.1 le 13 dec 98
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++   A domain of this category 
++   implements formal linear combinations
++   of elements from a domain \spad{Basis} with coefficients
++   in a domain \spad{R}. The domain \spad{Basis} needs only
++   to belong to the category \spadtype{SetCategory} and \spad{R}
++   to the category \spadtype{Ring}. Thus the coefficient ring
++   may be non-commutative.
++   See the \spadtype{XDistributedPolynomial} constructor
++   for examples of domains built with the \spadtype{FreeModuleCat}
++   category constructor.
++   Author: Michel Petitot (petitot@lifl.fr)

FreeModuleCat(R, Basis):Category == Exports where
   R: Ring
   Basis: SetCategory
   TERM ==> Record(k: Basis, c: R)
   
   Exports == Join(BiModule(R,R), RetractableTo Basis) with
        *                : (R, Basis) -> %
          ++ \spad{r*b} returns the product of \spad{r} by \spad{b}.
        coefficient        : (%, Basis) -> R
          ++ \spad{coefficient(x,b)} returns the coefficient 
          ++ of \spad{b} in \spad{x}.
        map                : (R -> R, %) -> %
          ++ \spad{map(fn,u)} maps function \spad{fn} onto the coefficients
          ++  of the non-zero monomials of \spad{u}.
        monom              : (Basis, R) -> %
          ++ \spad{monom(b,r)} returns the element with the single monomial
          ++  \spad{b} and coefficient \spad{r}.
        monomial?          : % -> Boolean
          ++ \spad{monomial?(x)} returns true if \spad{x} contains a single 
          ++ monomial.
        ListOfTerms        : % -> List TERM
          ++ \spad{ListOfTerms(x)} returns a list \spad{lt} of terms with type
          ++ \spad{Record(k: Basis, c: R)} such that \spad{x} equals
          ++ \spad{reduce(+, map(x +-> monom(x.k, x.c), lt))}.
        coefficients       : % -> List R           
          ++ \spad{coefficients(x)} returns the list of coefficients of \spad{x}.
        monomials          : % -> List %
          ++ \spad{monomials(x)} returns the list of \spad{r_i*b_i}
          ++ whose sum is \spad{x}.
        numberOfMonomials  : % -> NonNegativeInteger
          ++ \spad{numberOfMonomials(x)} returns the number of monomials of \spad{x}.
        leadingMonomial    : % -> Basis
          ++ \spad{leadingMonomial(x)} returns the first element from \spad{Basis}
          ++ which appears in \spad{ListOfTerms(x)}.
        leadingCoefficient : % -> R
          ++ \spad{leadingCoefficient(x)} returns the first coefficient
          ++ which appears in \spad{ListOfTerms(x)}.
        leadingTerm        : % -> TERM 
          ++ \spad{leadingTerm(x)} returns the first term which
          ++ appears in \spad{ListOfTerms(x)}.
        reductum           : % -> %
          ++ \spad{reductum(x)} returns \spad{x} minus its leading term.

      -- attributs
        if R has CommutativeRing then Module(R)

import Ring
import OrderedSet
)abbrev domain FM1 FreeModule1
++ Author: Michel Petitot petitot@lifl.fr
++ Date Created: 91
++ Date Last Updated: 7 Juillet 92
++ Fix History: compilation v 2.1 le 13 dec 98
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++   This domain implements linear combinations
++   of elements from the domain \spad{S} with coefficients
++   in the domain \spad{R} where \spad{S} is an ordered set
++   and \spad{R} is a ring (which may be non-commutative).
++   This domain is used by domains of non-commutative algebra such as:
++       \spadtype{XDistributedPolynomial},
++       \spadtype{XRecursivePolynomial}.
++   Author: Michel Petitot (petitot@lifl.fr)

FreeModule1(R:Ring,S:OrderedSet): FMcat == FMdef where
  EX ==> OutputForm
  TERM ==> Record(k:S,c:R)

  FMcat == FreeModuleCat(R,S) with
    *:(S,R) -> %
      ++ \spad{s*r} returns the product \spad{r*s}
      ++ used by \spadtype{XRecursivePolynomial} 
  FMdef == FreeModule(R,S) add
    -- representation
      Rep := List TERM  

    -- declarations
      lt: List TERM 
      x : %
      r : R
      s : S

    -- define
      numberOfMonomials p ==
         # (p::Rep)

      ListOfTerms(x) == x:List TERM 

      leadingTerm x == x.first
      leadingMonomial x == x.first.k
      coefficients x == [t.c for t in x]
      monomials x == [ monom (t.k, t.c) for t in x]

      retractIfCan x ==
         not one? numberOfMonomials(x) => "failed"
         x.first.c = 1 => x.first.k
         "failed"

      coerce(s:S):% == [[s,1$R]]
      retract x ==
         (rr := retractIfCan x) case "failed" => error "FM1.retract impossible"
         rr :: S

      if R has noZeroDivisors then
         r * x  ==
             r = 0 => 0
             [[u.k,r * u.c]$TERM for u in x]
         x * r  == 
             r = 0 => 0
             [[u.k,u.c * r]$TERM for u in x]
       else
         r * x  ==
             r = 0 => 0
             [[u.k,a] for u in x | not (a:=r*u.c)= 0$R]
         x * r  ==
             r = 0 => 0
             [[u.k,a] for u in x | not (a:=u.c*r)= 0$R]

      r * s ==
        r = 0 => 0
        [[s,r]$TERM]

      s * r ==
        r = 0 => 0
        [[s,r]$TERM]

      monom(b,r):% == [[b,r]$TERM] 

      outTerm(r:R, s:S):EX ==
            r=1  => s::EX
            r::EX * s::EX

      coerce(a:%):EX ==
            empty? a => (0$R)::EX
            reduce(_+, reverse! [outTerm(t.c, t.k) for t in a])$List(EX)

      coefficient(x,s) ==
         null x => 0$R
         x.first.k > s => coefficient(rest x,s)
         x.first.k = s => x.first.c
         0$R

import Ring
import BiModule
import CommutativeRing
import Algebra
)abbrev category XALG XAlgebra
++ Author: Michel Petitot petitot@lifl.fr
++ Date Created: 91
++ Date Last Updated: 7 Juillet 92
++ Fix History: compilation v 2.1 le 13 dec 98
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++   This is the category of algebras over non-commutative rings.
++   It is used by constructors of non-commutative algebras such as:
++       \spadtype{XPolynomialRing}.
++       \spadtype{XFreeAlgebra}
++   Author: Michel Petitot (petitot@lifl.fr)

XAlgebra(R: Ring): Category == 
  Join(Ring, BiModule(R,R),CoercibleFrom R) with
    -- attributs
      if R has CommutativeRing then Algebra(R)
      -- if R has CommutativeRing then Module(R) 
-- add
--  coerce(x:R):% == x * 1$%

import OrderedSet
import Ring
import XAlgebra
import RetractableTo
)abbrev category XFALG XFreeAlgebra
++ Author: Michel Petitot petitot@lifl.fr
++ Date Created: 91
++ Date Last Updated: 7 Juillet 92
++ Fix History: compilation v 2.1 le 13 dec 98
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++    This category specifies opeations for  polynomials
++    and formal series with non-commutative variables.
++ Author: Michel Petitot (petitot@lifl.fr)

XFreeAlgebra(vl:OrderedSet,R:Ring):Category == Catdef where
   WORD   ==> OrderedFreeMonoid(vl)          -- monoide libre
   NNI    ==> NonNegativeInteger
   I      ==> Integer
   TERM   ==> Record(k: WORD, c: R)

   Catdef == Join(Ring, XAlgebra(R), RetractableTo WORD) 
     with
       *: (vl,%) -> %
         ++ \spad{v * x} returns the product of a variable \spad{x} by \spad{x}.
       *: (%, R) -> %                 
         ++ \spad{x * r} returns the product of \spad{x} by \spad{r}.
         ++ Usefull if \spad{R} is a non-commutative Ring.
       mindeg: % -> WORD                
         ++ \spad{mindeg(x)} returns the little word which appears in \spad{x}.
         ++ Error if \spad{x=0}.
       mindegTerm: % -> TERM 
         ++ \spad{mindegTerm(x)} returns the term whose word is \spad{mindeg(x)}.
       coef  : (%,WORD) -> R            
         ++ \spad{coef(x,w)} returns the coefficient of the word \spad{w} in \spad{x}. 
       coef  : (%,%) -> R
         ++ \spad{coef(x,y)} returns scalar product of \spad{x} by \spad{y},
         ++ the set of words being regarded as an orthogonal basis.
       lquo  : (%,vl) -> %              
         ++ \spad{lquo(x,v)} returns the left simplification of \spad{x} by the variable \spad{v}.
       lquo  : (%,WORD) -> %            
         ++ \spad{lquo(x,w)} returns the left simplification of \spad{x} by the word \spad{w}.
       lquo  : (%,%) -> %
         ++ \spad{lquo(x,y)} returns the left simplification of \spad{x} by \spad{y}.
       rquo  : (%,vl) -> %
         ++ \spad{rquo(x,v)} returns the right simplification of \spad{x} by the variable \spad{v}.
       rquo  : (%,WORD) -> %
         ++ \spad{rquo(x,w)} returns the right simplification of \spad{x} by \spad{w}.
       rquo  : (%,%) -> %
         ++ \spad{rquo(x,y)} returns the right simplification of \spad{x} by \spad{y}.
       monom : (WORD , R) -> %
         ++ \spad{monom(w,r)} returns the product of the word \spad{w} by the coefficient \spad{r}.
       monomial? : % -> Boolean
         ++ \spad{monomial?(x)} returns true if \spad{x} is a monomial
       mirror: % -> %                   
         ++ \spad{mirror(x)} returns \spad{Sum(r_i mirror(w_i))} if \spad{x} writes \spad{Sum(r_i w_i)}. 
       coerce : vl -> %
         ++ \spad{coerce(v)} returns \spad{v}.
       constant?:% -> Boolean
         ++ \spad{constant?(x)} returns true if \spad{x} is constant.
       constant: % -> R   
         ++ \spad{constant(x)} returns the constant term of \spad{x}.
       quasiRegular? : % -> Boolean  
         ++ \spad{quasiRegular?(x)} return true if \spad{constant(x)} is zero. 
       quasiRegular : % -> %
         ++ \spad{quasiRegular(x)} return \spad{x} minus its constant term.
       if R has CommutativeRing then
          sh :(%,%) -> %
             ++ \spad{sh(x,y)} returns the shuffle-product of \spad{x} by \spad{y}.
             ++ This multiplication is associative and commutative.
          sh :(%,NNI) -> %
             ++ \spad{sh(x,n)} returns the shuffle power of \spad{x} to the \spad{n}.
       map   : (R -> R, %) -> %
         ++ \spad{map(fn,x)} returns \spad{Sum(fn(r_i) w_i)} if \spad{x} writes \spad{Sum(r_i w_i)}.
       varList: % -> List vl
         ++ \spad{varList(x)} returns the list of variables which appear in \spad{x}.

     -- Attributs
       if R has noZeroDivisors then noZeroDivisors

import OrderedSet
import XFreeAlgebra
)abbrev category XPOLYC XPolynomialsCat
++ Author: Michel Petitot petitot@lifl.fr
++ Date Created: 91
++ Date Last Updated: 7 Juillet 92
++ Fix History: compilation v 2.1 le 13 dec 98
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++   The Category of polynomial rings with non-commutative variables.
++   The coefficient ring may be non-commutative too. 
++   However coefficients commute with vaiables.
++ Author: Michel Petitot (petitot@lifl.fr)

XPolynomialsCat(vl:OrderedSet,R:Ring):Category == Export where
  WORD ==> OrderedFreeMonoid(vl)

  Export == XFreeAlgebra(vl,R) with
    maxdeg: % -> WORD 
      ++ \spad{maxdeg(p)} returns the greatest leading word in the support of \spad{p}.
    degree: % -> NonNegativeInteger 
      ++ \spad{degree(p)} returns the degree of \spad{p}. 
      ++  Note that the degree of a word is its length. 
    trunc : (% , NonNegativeInteger) -> %
      ++  \spad{trunc(p,n)} returns the polynomial \spad{p} truncated at order \spad{n}.

import Ring
import OrderedMonoid
import XAlgebra
import FreeMonoidCat
)abbrev domain XPR XPolynomialRing
++ Author: Michel Petitot petitot@lifl.fr
++ Date Created: 91
++ Date Last Updated: 7 Juillet 92
++ Fix History: compilation v 2.1 le 13 dec 98
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++ This domain represents generalized polynomials with coefficients
++ (from a not necessarily commutative ring), and words
++ belonging to an arbitrary \spadtype{OrderedMonoid}.
++ This type is used, for instance, by the \spadtype{XDistributedPolynomial} 
++ domain constructor where the Monoid is free.
++ Author: Michel Petitot (petitot@lifl.fr)

XPolynomialRing(R:Ring,E:OrderedMonoid): T == C where
  TERM   ==> Record(k: E, c: R)
  EX     ==> OutputForm
  NNI    ==> NonNegativeInteger

  T == Join(Ring, XAlgebra(R), FreeModuleCat(R,E),CoercibleFrom E) with
    --operations
      *: (%,R) -> %
        ++ \spad{p*r} returns the product of \spad{p} by \spad{r}.
      #: % -> NonNegativeInteger
        ++ \spad{# p} returns the number of terms in \spad{p}.
      maxdeg: % -> E
        ++ \spad{maxdeg(p)} returns the greatest word occurring in the polynomial \spad{p}
        ++ with a non-zero coefficient. An error is produced if  \spad{p} is zero.
      mindeg: % -> E
        ++ \spad{mindeg(p)} returns the smallest word occurring in the polynomial \spad{p}
        ++ with a non-zero coefficient. An error is produced if  \spad{p} is zero.
      reductum : % -> %   
        ++ \spad{reductum(p)} returns \spad{p} minus its leading term.
        ++ An error is produced if  \spad{p} is zero.
      coef  : (%,E) -> R
        ++ \spad{coef(p,e)} extracts the coefficient of the monomial \spad{e}.
        ++ Returns zero if \spad{e} is not present. 
      constant?:% -> Boolean
        ++ \spad{constant?(p)} tests whether the polynomial \spad{p} belongs to the
        ++ coefficient ring.
      constant: % -> R
        ++ \spad{constant(p)} return the constant term of \spad{p}.
      quasiRegular? : % -> Boolean
        ++ \spad{quasiRegular?(x)} return true if \spad{constant(p)} is zero.
      quasiRegular : % -> % 
        ++ \spad{quasiRegular(x)} return \spad{x} minus its constant term.
      map   : (R -> R, %) -> %
        ++ \spad{map(fn,x)} returns \spad{Sum(fn(r_i) w_i)} if \spad{x} writes \spad{Sum(r_i w_i)}.
      if R has Field then / : (%,R) -> %
        ++ \spad{p/r} returns \spad{p*(1/r)}.

    --assertions
      if R has noZeroDivisors then noZeroDivisors
      if R has unitsKnown then unitsKnown
      if R has canonicalUnitNormal then canonicalUnitNormal
          ++ canonicalUnitNormal guarantees that the function
          ++ unitCanonical returns the same representative for all
          ++ associates of any particular element.


  C == FreeModule1(R,E) add
    --representations
       Rep:=  List TERM
    --uses
       repeatMultExpt: (%,NonNegativeInteger) -> %
    --define
       1  == [[1$E,1$R]]
 
       characteristic  == characteristic$R
       #x == #$Rep x
       maxdeg p == if null p then  error " polynome nul !!"
                             else p.first.k
       mindeg p == if null p then  error " polynome nul !!" 
                             else (last p).k
       
       coef(p,e)  ==
          for tm in p repeat
            tm.k=e => return tm.c
            tm.k < e => return 0$R
          0$R

       constant? p == (p = 0) or (maxdeg(p) = 1$E)
       constant  p == coef(p,1$E)

       quasiRegular? p == (p=0) or not one?((last p).k)
       quasiRegular  p == 
          quasiRegular?(p) => p
          [t for t in p | not(t.k = 1$E)]

       recip(p) ==
           p=0 => "failed"
           p.first.k > 1$E => "failed"
           (u:=recip(p.first.c)) case "failed" => "failed"
           (u::R)::%
 
       coerce(r:R) == if r=0$R then 0$% else [[1$E,r]]
       coerce(n:Integer) == (n::R)::%
 
       if R has noZeroDivisors then
         p1:% * p2:%  ==
            null p1 => 0
            null p2 => 0
            p1.first.k = 1$E => p1.first.c * p2
            p2 = 1 => p1
--            +/[[[t1.k*t2.k,t1.c*t2.c]$TERM for t2 in p2]
--                   for t1 in reverse(p1)]
            +/[[[t1.k*t2.k,t1.c*t2.c]$TERM for t2 in p2]
                   for t1 in p1]
        else
         p1:% * p2:%  ==
            null p1 => 0
            null p2 => 0
            p1.first.k = 1$E => p1.first.c * p2
            p2 = 1 => p1
--            +/[[[t1.k*t2.k,r]$TERM for t2 in p2 | not (r:=t1.c*t2.c) =$R 0]
--                 for t1 in reverse(p1)]
            +/[[[t1.k*t2.k,r]$TERM for t2 in p2 | not (r:=t1.c*t2.c) =$R 0]
                   for t1 in p1]
       p:% ** nn:NNI  == repeatMultExpt(p,nn)
       repeatMultExpt(x,nn) ==
               nn = 0 => 1
               y:% := x
               for i in 2..nn repeat y:= x * y
               y
              
       outTerm(r:R, m:E):EX ==
            r=1 => m::EX
            m=1 => r::EX
            r::EX * m::EX

--       coerce(x:%) : EX ==
--         null x => (0$R) :: EX
--         le : List EX := nil
--         for rec in x repeat
--           rec.c = 1$R => le := cons(rec.k :: EX, le)
--           rec.k = 1$E => le := cons(rec.c :: EX, le)
--           le := cons(mkBinary("*"::EX,rec.c :: EX,
--             rec.k :: EX), le)
--         1 = #le => first le
--         mkNary("+" :: EX,le)

       coerce(a:%):EX ==
            empty? a => (0$R)::EX
            reduce(_+, reverse! [outTerm(t.c, t.k) for t in a])$List(EX)

 
       if R has Field then
          x/r == inv(r)*x

import OrderedSet
import Ring
import FreeModuleCat
import XPolynomialRing
import XPolynomialsCat
)abbrev domain XDPOLY XDistributedPolynomial
++ Author: Michel Petitot petitot@lifl.fr
++ Date Created: 91
++ Date Last Updated: 7 Juillet 92
++ Fix History: compilation v 2.1 le 13 dec 98
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++   This type supports distributed multivariate polynomials
++ whose variables do not commute.
++ The coefficient ring may be non-commutative too.
++ However, coefficients and variables commute.
++ Author: Michel Petitot (petitot@lifl.fr)

XDistributedPolynomial(vl:OrderedSet,R:Ring): XDPcat == XDPdef where

  WORD ==> OrderedFreeMonoid(vl)
  I    ==> Integer
  NNI  ==> NonNegativeInteger
  TERM ==> Record(k:WORD, c:R)

  XDPcat == Join(FreeModuleCat(R, WORD), XPolynomialsCat(vl,R))

  XDPdef == XPolynomialRing(R,WORD) add

       import( WORD, TERM)

    -- Representation
       Rep  :=  List TERM

    -- local functions
       shw: (WORD , WORD) -> %    -- shuffle de 2 mots

    -- definitions

       mindegTerm p == last(p)$Rep

       if R has CommutativeRing then
         sh(p:%, n:NNI):% ==
            n=0 => 1
            n=1 => p
            n1: NNI := (n-$I 1)::NNI
            sh(p, sh(p,n1))

      
         sh(p1:%, p2:%) ==
           p:% := 0 
           for t1 in p1 repeat
             for t2 in p2 repeat
                p := p + (t1.c * t2.c) * shw(t1.k,t2.k) 
           p

       coerce(v: vl):% == coerce(v::WORD)
       v:vl * p:% ==
         [[v * t.k , t.c]$TERM for t in p]

       mirror p == 
         null p => p
         monom(mirror$WORD leadingMonomial p, leadingCoefficient p) + _
               mirror reductum p

       degree(p) == length(maxdeg(p))$WORD

       trunc(p, n) ==
         p = 0 => p
         degree(p) > n => trunc( reductum p , n)
         p

       varList p ==
         constant? p => []
         le : List vl := "setUnion"/[varList(t.k) for t in p]
         sort!(le)

       rquo(p:% , w: WORD) == 
         [[r::WORD,t.c]$TERM for t in p | not (r:= rquo(t.k,w)) case "failed" ]
       lquo(p:% , w: WORD) ==
         [[r::WORD,t.c]$TERM for t in p | not (r:= lquo(t.k,w)) case "failed" ]
       rquo(p:% , v: vl) ==
         [[r::WORD,t.c]$TERM for t in p | not (r:= rquo(t.k,v)) case "failed" ]
       lquo(p:% , v: vl) ==
         [[r::WORD,t.c]$TERM for t in p | not (r:= lquo(t.k,v)) case "failed" ]

       shw(w1,w2) ==
         w1 = 1$WORD => w2::%
         w2 = 1$WORD => w1::%
         x: vl := first w1 ; y: vl := first w2
         x * shw(rest w1,w2) + y * shw(w1,rest w2)
 
       lquo(p:%,q:%):% ==
         +/  [r * t.c for t in q | (r := lquo(p,t.k)) ~= 0] 

       rquo(p:%,q:%):% ==
         +/  [r * t.c for t in q | (r := rquo(p,t.k)) ~= 0] 

       coef(p:%,q:%):R ==
         p = 0 => 0$R
         q = 0 => 0$R 
         p.first.k > q.first.k => coef(p.rest,q)
         p.first.k < q.first.k => coef(p,q.rest) 
         return p.first.c * q.first.c + coef(p.rest,q.rest)

import OrderedSet
import Ring
import XPolynomialsCat
import XDistributedPolynomial
)abbrev domain XRPOLY XRecursivePolynomial
++ Author: Michel Petitot petitot@lifl.fr
++ Date Created: 91
++ Date Last Updated: 7 Juillet 92
++ Fix History: compilation v 2.1 le 13 dec 98
++   extend renomme en expand
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++   This type supports multivariate polynomials
++ whose variables do not commute.
++ The representation is recursive.
++ The coefficient ring may be non-commutative.
++ Coefficients and variables commute.
++ Author: Michel Petitot (petitot@lifl.fr)

XRecursivePolynomial(VarSet:OrderedSet,R:Ring):  Xcat == Xdef where
  I      ==> Integer
  NNI    ==> NonNegativeInteger
  XDPOLY ==> XDistributedPolynomial(VarSet, R)
  EX     ==> OutputForm
  WORD   ==> OrderedFreeMonoid(VarSet)
  TERM   ==> Record(k:VarSet , c:%)
  LTERMS ==> List(TERM) 
  REGPOLY==> FreeModule1(%, VarSet) 
  VPOLY  ==> Record(c0:R, reg:REGPOLY)

  Xcat == XPolynomialsCat(VarSet,R) with
       expand: % -> XDPOLY
         ++ \spad{expand(p)} returns \spad{p} in distributed form.
       unexpand : XDPOLY -> %
         ++ \spad{unexpand(p)} returns \spad{p} in recursive form.
       RemainderList: % -> LTERMS
         ++ \spad{RemainderList(p)} returns the regular part of \spad{p}
         ++ as a list of terms.

  Xdef == add
       import(VPOLY)

    -- representation
       Rep     := Union(R,VPOLY)

    -- local functions
       construct: LTERMS -> REGPOLY
       simplifie: VPOLY -> %
       lquo1: (LTERMS,LTERMS) -> %        ++ a ajouter
       coef1: (LTERMS,LTERMS) -> R        ++ a ajouter
       outForm: REGPOLY -> EX

    --define
       construct(lt) == lt pretend REGPOLY
       p1:%  =  p2:%  ==
         p1 case R =>
             p2 case R => p1 =$R p2
             false
         p2 case R => false
         p1.c0 =$R p2.c0 and p1.reg =$REGPOLY p2.reg

       monom(w, r) == 
         r =0 => 0
         r * w::%

--       if R has Field then                  -- Bug non resolu !!!!!!!!
--         p:% / r: R == inv(r) * p
 
       rquo(p1:%, p2:%):% ==
         p2 case R => p1 * p2::R
         p1 case R => p1  * p2.c0
         x:REGPOLY := construct [[t.k, a]$TERM for t in ListOfTerms(p1.reg) _
                         | (a:= rquo(t.c,p2)) ~= 0$% ]$LTERMS
         simplifie [coef(p1,p2) , x]$VPOLY

       trunc(p,n) ==
         n = 0 or (p case R) => (constant p)::%
         n1: NNI := (n-1)::NNI
         lt: LTERMS := [[t.k, r]$TERM for t in ListOfTerms p.reg _
                        | (r := trunc(t.c, n1)) ~= 0]$LTERMS
         x: REGPOLY := construct lt
         simplifie [constant p, x]$VPOLY

       unexpand p ==
         constant? p => (constant p)::%
         vl: List VarSet := sort(#1 > #2, varList p)
         x : REGPOLY := _
           construct [[v, unexpand r]$TERM for v in vl| (r:=lquo(p,v)) ~= 0]
         [constant p, x]$VPOLY

       if R has CommutativeRing then
         sh(p:%, n:NNI):% ==
            n = 0 => 1
            p case R => (p::R)** n
            n1: NNI := (n-1)::NNI
            p1: % := n * sh(p, n1)  
            lt: LTERMS := [[t.k, sh(t.c, p1)]$TERM for t in ListOfTerms p.reg]
            [p.c0 ** n, construct lt]$VPOLY
 
         sh(p1:%, p2:%) ==
            p1 case R => p1::R * p2
            p2 case R => p1 * p2::R 
            lt1:LTERMS := ListOfTerms p1.reg ; lt2:LTERMS := ListOfTerms p2.reg
            x: REGPOLY := construct [[t.k,sh(t.c,p2)]$TERM for t in lt1]
            y: REGPOLY := construct [[t.k,sh(p1,t.c)]$TERM for t in lt2]
            [p1.c0*p2.c0,x + y]$VPOLY

       RemainderList p == 
           p case R => []
           ListOfTerms( p.reg)$REGPOLY
 
       lquo(p1:%,p2:%):% ==
         p2 case R => p1 * p2
         p1 case R => p1  *$R p2.c0
         p1 * p2.c0 +$% lquo1(ListOfTerms p1.reg, ListOfTerms p2.reg)

       lquo1(x:LTERMS,y:LTERMS):% ==
         null x => 0$%  
         null y => 0$%
         x.first.k < y.first.k => lquo1(x,y.rest)
         x.first.k = y.first.k => 
             lquo(x.first.c,y.first.c) + lquo1(x.rest,y.rest)
         return lquo1(x.rest,y)

       coef(p1:%, p2:%):R ==
         p1 case R => p1::R * constant p2
         p2 case R => p1.c0 * p2::R
         p1.c0 * p2.c0 +$R coef1(ListOfTerms p1.reg, ListOfTerms p2.reg)

       coef1(x:LTERMS,y:LTERMS):R ==
         null x => 0$R
         null y => 0$R
         x.first.k < y.first.k => coef1(x,y.rest)
         x.first.k = y.first.k =>
             coef(x.first.c,y.first.c) + coef1(x.rest,y.rest)
         return coef1(x.rest,y)

       --------------------------------------------------------------
       outForm(p:REGPOLY): EX ==
          le : List EX :=  [t.k::EX * t.c::EX for t in ListOfTerms p]
          reduce(_+, reverse! le)$List(EX)

       coerce(p:$): EX ==
          p case R => (p::R)::EX
          p.c0 = 0 => outForm p.reg
          p.c0::EX + outForm p.reg 

       0 == 0$R::%
       1 == 1$R::%
       constant? p ==  p case R
       constant p == 
          p case R => p
          p.c0

       simplifie p ==
         p.reg = 0$REGPOLY => (p.c0)::%
         p

       coerce (v:VarSet):% ==
         [0$R,coerce(v)$REGPOLY]$VPOLY

       coerce (r:R):% == r::%
       coerce (n:Integer) == n::R::%
       coerce (w:WORD) == 
         w = 1 => 1$R
         (first w) * coerce(rest w)
 
       expand p ==
         p case R => p::R::XDPOLY
         lt:LTERMS := ListOfTerms(p.reg)
         ep:XDPOLY := (p.c0)::XDPOLY
         for t in lt repeat
           ep:= ep + t.k * expand(t.c)
         ep
                
       - p:% ==
         p case R => -$R p
         [- p.c0, - p.reg]$VPOLY
 
       p1 + p2 ==
         p1 case R and p2 case R => p1 +$R p2
         p1 case R => [p1 + p2.c0 , p2.reg]$VPOLY
         p2 case R => [p2 + p1.c0 , p1.reg]$VPOLY 
         simplifie [p1.c0 + p2.c0 , p1.reg +$REGPOLY p2.reg]$VPOLY
 
       p1 - p2 ==
         p1 case R and p2 case R => p1 -$R p2
         p1 case R => [p1 - p2.c0 , -p2.reg]$VPOLY
         p2 case R => [p1.c0 - p2 , p1.reg]$VPOLY
         simplifie [p1.c0 - p2.c0 , p1.reg -$REGPOLY p2.reg]$VPOLY
 
       n:Integer * p:% ==
         n=0 => 0$%
         p case R => n *$R p
         -- [ n*p.c0,n*p.reg]$VPOLY
         simplifie [ n*p.c0,n*p.reg]$VPOLY

       r:R * p:% ==
         r=0 => 0$%
         p case R => r *$R p
         -- [ r*p.c0,r*p.reg]$VPOLY
         simplifie [ r*p.c0,r*p.reg]$VPOLY

       p:% * r:R ==
         r=0 => 0$%
         p case R => p *$R r
         -- [ p.c0 * r,p.reg * r]$VPOLY
         simplifie [ r*p.c0,r*p.reg]$VPOLY

       v:VarSet * p:% == 
          p = 0 => 0$%
          [0$R, v *$REGPOLY p]$VPOLY
 
       p1:% * p2:% ==
         p1 case R => p1::R * p2
         p2 case R => p1 * p2::R
         x:REGPOLY := p1.reg *$REGPOLY p2
         y:REGPOLY := (p1.c0)::% *$REGPOLY p2.reg  -- maladroit:(p1.c0)::% !!
         -- [ p1.c0 * p2.c0 , x+y ]$VPOLY
         simplifie [ p1.c0 * p2.c0 , x+y ]$VPOLY

       lquo(p:%, v:VarSet):% ==
         p case R => 0
         coefficient(p.reg,v)$REGPOLY

       lquo(p:%, w:WORD):% ==
         w = 1$WORD => p
         lquo(lquo(p,first w),rest w)

       rquo(p:%, v:VarSet):% ==
         p case R => 0
         x:REGPOLY := construct [[t.k, a]$TERM for t in ListOfTerms(p.reg)
                         | (a:= rquo(t.c,v)) ~= 0 ]
         simplifie [constant(coefficient(p.reg,v)) , x]$VPOLY 
        
       rquo(p:%, w:WORD):% ==
         w = 1$WORD => p
         rquo(rquo(p,rest w),first w)
 
       coef(p:%, w:WORD):R ==
         constant lquo(p,w)

       quasiRegular? p == 
         p case R => p = 0$R
         p.c0 = 0$R

       quasiRegular p ==
         p case R => 0$%
         [0$R,p.reg]$VPOLY

       characteristic == characteristic$R
       recip p ==
         p case R => recip(p::R)
         "failed"

       mindeg p ==
         p case R =>
           p = 0 => error "XRPOLY.mindeg: polynome nul !!"
           1$WORD
         p.c0 ~= 0 => 1$WORD
         "min"/[(t.k) *$WORD mindeg(t.c) for t in ListOfTerms p.reg] 

       maxdeg p ==
         p case R => 
            p = 0 => error "XRPOLY.maxdeg: polynome nul !!"
            1$WORD
         "max"/[(t.k) *$WORD maxdeg(t.c) for t in ListOfTerms p.reg] 

       degree p == 
          p = 0 => error "XRPOLY.degree: polynome nul !!"
          length(maxdeg p)

       map(fn,p) ==
         p case R => fn(p::R)
         x:REGPOLY := construct [[t.k,a]$TERM for t in ListOfTerms p.reg
                         |(a := map(fn,t.c)) ~= 0$R]
         simplifie [fn(p.c0),x]$VPOLY

       varList p ==
         p case R => []
         lv: List VarSet := "setUnion"/[varList(t.c) for t in ListOfTerms p.reg]
         lv:= setUnion(lv,[t.k for t in ListOfTerms p.reg])
         sort!(lv)

import XRecursivePolynomial
)abbrev domain XPOLY XPolynomial
++ Author: Michel Petitot petitot@lifl.fr
++ Date Created: 91
++ Date Last Updated: 7 Juillet 92
++ Fix History: compilation v 2.1 le 13 dec 98
++   extend renomme en expand
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++   This type supports multivariate polynomials
++ whose set of variables is \spadtype{Symbol}.
++ The representation is recursive.
++ The coefficient ring may be non-commutative and the variables 
++ do not commute.
++ However, coefficients and variables commute.
++ Author: Michel Petitot (petitot@lifl.fr)

XPolynomial(R:Ring) == XRecursivePolynomial(Symbol, R)

