\documentclass{article}% ===> this file was generated automatically by noweave --- better not edit it
\usepackage{axiom}
\begin{document}
\title{\$SPAD/src/algebra poly.spad}
\author{Dave Barton, James Davenport, Barry Trager, Patrizia Gianni, Marc Moreno Maza}
\date{June 2009}
\maketitle
\begin{abstract}
FreeModuleCategory

FreeModule

PolynomialRing

SparseUnivariatePolynomial

SparseUnivariatePolynomialFunctions2

UnivariatePolynomial

UnivariatePolynomialFunctions2

PolynomialToUnivariatePolynomial

UnivariatePolynomialSquareFree

PolynomialSquareFree

UnivariatePolynomialMultiplicationPackage

\end{abstract}
\eject
\tableofcontents
\eject
FreeModuleCategory

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{FreeModuleCategory}
category constructor.
Author: Michel Petitot (petitot@lifl.fr)Note (Franz Lehner, June 2009):
Since \spad{leadingTerm} makes no sense
for unordered base sets,
and at the time of this writing this domain was never used for such,
it was changed to \spad{OrderedSet}.
\spad{FreeModule} originally was not of FreeModuleCategory.
Some functions (like \spad{support}, \spad{coefficients},
\spad{monomials}, ...) from here could be moved to
\spad{IndexedDirectProductCategory}
but at the moment there is no need for this.

FreeModule

A bi-module is a free module
over a ring with generators indexed by an ordered set.
Each element can be expressed as a finite linear combination of
generators. Only non-zero terms are stored.

old domain FreeModule1 was merged to it in May 2009
The description of the latter:

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

PolynomialRing

This domain represents generalized polynomials with coefficients
(from a not necessarily commutative ring), and terms
indexed by their exponents (from an arbitrary ordered abelian monoid).
This type is used, for example,
by the \spadtype{DistributedMultivariatePolynomial} domain where
the exponent domain is a direct product of non negative integers.

SparseUnivariatePolynomial

This domain represents univariate polynomials over arbitrary
(not necessarily commutative) coefficient rings. The variable is
unspecified  so that the variable displays as \spad{?} on output.
If it is necessary to specify the variable name, use type \spadtype{UnivariatePolynomial}.
The representation is sparse
in the sense that only non-zero terms are represented.
Note: if the coefficient ring is a field, this domain forms a euclidean domain.

SparseUnivariatePolynomialFunctions2

This package lifts a mapping from coefficient rings R to S to
a mapping from sparse univariate polynomial over R to
a sparse univariate polynomial over S.
Note that the mapping is assumed
to send zero to zero, since it will only be applied to the non-zero
coefficients of the polynomial.

UnivariatePolynomial

This domain represents univariate polynomials in some symbol
over arbitrary (not necessarily commutative) coefficient rings.
The representation is sparse
in the sense that only non-zero terms are represented.
Note: if the coefficient ring is a field, then this domain forms a euclidean domain.

UnivariatePolynomialFunctions2

This package lifts a mapping from coefficient rings R to S to
a mapping from \spadtype{UnivariatePolynomial}(x,R) to
\spadtype{UnivariatePolynomial}(y,S). Note that the mapping is assumed
to send zero to zero, since it will only be applied to the non-zero
coefficients of the polynomial.

PolynomialToUnivariatePolynomial

This package is primarily to help the interpreter do coercions.
It allows you to view a polynomial as a
univariate polynomial in one of its variables with
coefficients which are again a polynomial in all the
other variables.

UnivariatePolynomialSquareFree

This package provides for square-free decomposition of
univariate polynomials over arbitrary rings, i.e.
a partial factorization such that each factor is a product
of irreducibles with multiplicity one and the factors are
pairwise relatively prime. If the ring
has characteristic zero, the result is guaranteed to satisfy
this condition. If the ring is an infinite ring of
finite characteristic, then it may not be possible to decide when
polynomials contain factors which are pth powers. In this
case, the flag associated with that polynomial is set to "nil"
(meaning that that polynomials are not guaranteed to be square-free).

PolynomialSquareFree

This package computes square-free decomposition of multivariate
polynomials over a coefficient ring which is an arbitrary gcd domain.
The requirement on the coefficient domain guarantees that the \spadfun{content} can be
removed so that factors will be primitive as well as square-free.
Over an infinite ring of finite characteristic,it may not be possible to
guarantee that the factors are square-free.

UnivariatePolynomialMultiplicationPackage

This package implements Karatsuba's trick for multiplying
(large) univariate polynomials. It could be improved with
a version doing the work on place and also with a special
case for squares. We've done this in Basicmath, but we
believe that this out of the scope of FriCAS.
Currently unused.

\section{category FMCAT FreeModuleCategory}
\nwfilename{PolySpad.pamphlet}\nwbegincode{1}\moddef{category FMCAT FreeModuleCategory}\endmoddef\nwstartdeflinemarkup\nwenddeflinemarkup
)abbrev category FMCAT FreeModuleCategory
++ 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\{FreeModuleCategory\}
++   category constructor.
++   Author: Michel Petitot (petitot@lifl.fr)
++
++   Note (Franz Lehner, June 2009):
++   Since \\spad\{leadingTerm\} makes no sense
++   for unordered base sets,
++   and at the time of this writing this domain was never used for such,
++   it was changed to \\spad\{OrderedSet\}.
++   \\spad\{FreeModule\} originally was not of FreeModuleCategory.
++   Some functions (like \\spad\{support\}, \\spad\{coefficients\},
++   \\spad\{monomials\}, ...) from here could be moved to
++   \\spad\{IndexedDirectProductCategory\}
++   but at the moment there is no need for this.
FreeModuleCategory(R, S):Category == Exports where
   R: Ring
   S: OrderedSet
   Term ==> Record(k: S, c: R)

   Exports == Join(BiModule(R,R), RetractableTo S, IndexedDirectProductCategory(R,S)) with
        "*"                : (R, S) -> %
          ++ \\spad\{r*b\} returns the product of \\spad\{r\} by \\spad\{b\}.
        "*":(S,R) -> %
          ++ \\spad\{s*r\} returns the product \\spad\{r*s\}
          ++ used by \\spadtype\{XRecursivePolynomial\}
        coefficients       : % -> List R
          ++ \\spad\{coefficients(x)\} returns the list of coefficients of \\spad\{x\}.
        support            : % -> List S
          ++ \\spad\{support(x)\} returns the list of basis elements with nonzero coefficients.
        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\}.
        monomial?: % -> Boolean
          ++ \\spad\{monomial?(x)\} returns true if \\spad\{x\} contains a single
          ++ monomial.
        leadingMonomial    : % -> S
          ++ \\spad\{leadingMonomial(x)\} returns the first element from \\spad\{S\}
          ++ which appears in \\spad\{listOfTerms(x)\}.
        leadingCoefficient : % -> R
          ++ \\spad\{leadingCoefficient(x)\} returns the first coefficient
          ++ which appears in \\spad\{listOfTerms(x)\}.
        monom : (S, R) -> %
          ++ \\spad\{monom(s,r)\} returns the product of the basis element \\spad\{s\} by the coefficient \\spad\{r\}.
        coefficient :(%,S) -> R
          ++ \\spad\{coefficient(x,s)\} returns the coefficient of the basis element s
        -- attributs
        if R has CommutativeRing then
             Module(R)
             linearExtend:(S->R,%)->R
          ++ \\spad\{linearExtend:(f,x)\} returns the linear extension
          ++ of a map defined on the basis applied to a linear combination
        if R has Comparable then Comparable
      add
        if R has Comparable then
          smaller?(p:%,q:%):Boolean ==
            if leadingMonomial(p)=leadingMonomial(q) then
              if leadingCoefficient(p)=leadingCoefficient(q) then
                smaller?(reductum p, reductum q)
              else
                smaller?(leadingCoefficient(p),leadingCoefficient(q))
            leadingMonomial(p)<leadingMonomial(q)

\nwendcode{}\nwbegindocs{2}\nwdocspar
\section{domain FM FreeModule}
\nwenddocs{}\nwbegincode{3}\moddef{domain FM FreeModule}\endmoddef\nwstartdeflinemarkup\nwenddeflinemarkup
)abbrev domain FM FreeModule
++ Author: Dave Barton, James Davenport, Barry Trager
++ Date Created:
++ Date Last Updated:
++ Basic Functions: BiModule(R,R)
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++ A bi-module is a free module
++ over a ring with generators indexed by an ordered set.
++ Each element can be expressed as a finite linear combination of
++ generators. Only non-zero terms are stored.

++ old domain FreeModule1 was merged to it in May 2009
++ The description of the latter:
++   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)


FreeModule(R:Ring,S:OrderedSet):
        Join(BiModule(R,R),FreeModuleCategory(R,S)) with
    if R has CommutativeRing then Module(R)
 == IndexedDirectProductAbelianGroup(R,S) add
    --representations
       Term ==>  Record(k:S,c:R)
       Rep :=  List Term
    --declarations
       x,y: %
       r: R
       n: Integer
       f: R -> R
       s: S
       lt: List Term
   --define
       if R has EntireRing then
         r * x  ==
             zero? r => 0
--             one? r => x
             (r = 1) => x
           --map(x+->r*x1,x)
             [[u.k,r*u.c] for u in x ]
       else
         r * x  ==
             zero? r => 0
--             one? r => x
             (r = 1) => x
           --map(x1+->r*x1,x)
             [[u.k,a] for u in x | (a:=r*u.c) ~= 0$R]
       if R has EntireRing then
         x * r  ==
             zero? r => 0
--             one? r => x
             (r = 1) => x
           --map(x1+->r*x1,x)
             [[u.k,u.c*r] for u in x ]
       else
         x * r  ==
             zero? r => 0
--             one? r => x
             (r = 1) => x
           --map(x1+->r*x1,x)
             [[u.k,a] for u in x | (a:=u.c*r) ~= 0$R]

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

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

       coerce(x) : OutputForm ==
         null x => (0$R) :: OutputForm
         le : List OutputForm := nil
         for rec in reverse x repeat
           rec.c = 1 => le := cons(rec.k :: OutputForm, le)
           le := cons(rec.c :: OutputForm *  rec.k :: OutputForm, le)
         reduce("+",le)

       leadingMonomial x == x.first.k

       support x == [t.k for t in x]

       coefficients x == [t.c for t in x]

       monomials x == [ monom (t.k, t.c) for t in x]

       retractIfCan x ==
         numberOfMonomials(x) ~= 1 => "failed"
         x.first.c = 1 => x.first.k
         "failed"

       retract x ==
         (rr := retractIfCan x) case "failed" => error "FM1.retract impossible"
         rr :: S

       coerce(s:S):% == [[s,1$R]]

-- the following is to be replaced by monomial(r,b) everywhere
       monom(b,r):% == [[b,r]$Term]

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

       monomial? x ==
         numberOfMonomials x = 1

       listOfTerms(x) ==
--          (x::Rep)
--          coerce(x)@Rep
         x pretend Rep

       numberOfMonomials x == # (listOfTerms x)

       if R has CommutativeRing then
           f:S->R
           x:%
           t:Term
           linearExtend(f,x) ==
               zero? x => 0
               res:R:= 0
               for t in listOfTerms x repeat
                   res := res + (t c)*f(t k)
               res


\nwendcode{}\nwbegindocs{4}\nwdocspar
\section{package FM2 FreeModuleFunctions2}
\nwenddocs{}\nwbegincode{5}\moddef{package FM2 FreeModuleFunctions2}\endmoddef\nwstartdeflinemarkup\nwenddeflinemarkup
)abbrev package FM2 FreeModuleFunctions2
++ Author: Franz Lehner lehner@finanz.math.tugraz.at
++ Date Created: 2009
++ Date Last Updated: 7 June 2009
++ Fix History:
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++   linear extensions of maps on the basis
FreeModuleFunctions2(R:CommutativeRing, S:OrderedSet,M1:FreeModuleCategory(R,S),M2:Module(R)): public == private where
    TERM1 ==> Record(k:S, c:R)

    public ==> with
        linearExtend:(S->M2,M1)->M2
        ++ \\spad\{linearExtend:(f,x)\} returns the linear extension
        ++ of a map defined on the basis of M2 applied to a linear combination

    private ==> add
        f:S->M2
        x:M1
        t:TERM1
        linearExtend(f,x) ==
            zero? x => 0
            res:M2:= 0
            for t in listOfTerms x repeat
                res := res + (t c)*f(t k)
            res


\nwendcode{}\nwbegindocs{6}\nwdocspar

\section{domain PR PolynomialRing}
\nwenddocs{}\nwbegincode{7}\moddef{domain PR PolynomialRing}\endmoddef\nwstartdeflinemarkup\nwenddeflinemarkup
)abbrev domain PR PolynomialRing
++ Author: Dave Barton, James Davenport, Barry Trager
++ Date Created:
++ Date Last Updated: 14.08.2000. Improved exponentiation [MMM/TTT]
++ Basic Functions: Ring, degree, coefficient, monomial, reductum
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++ This domain represents generalized polynomials with coefficients
++ (from a not necessarily commutative ring), and terms
++ indexed by their exponents (from an arbitrary ordered abelian monoid).
++ This type is used, for example,
++ by the \\spadtype\{DistributedMultivariatePolynomial\} domain where
++ the exponent domain is a direct product of non negative integers.

PolynomialRing(R:Ring,E:OrderedAbelianMonoid): T == C
 where
  T == FiniteAbelianMonoidRing(R,E) with
    --assertions
       if R has IntegralDomain and E has CancellationAbelianMonoid then
            fmecg: (%,E,R,%) -> %
            ++ fmecg(p1,e,r,p2) finds X : p1 - r * X^e * p2
       if R has canonicalUnitNormal then canonicalUnitNormal
          ++ canonicalUnitNormal guarantees that the function
          ++ unitCanonical returns the same representative for all
          ++ associates of any particular element.

  C == FreeModule(R,E) add
    --representations
       Term ==>  Record(k:E,c:R)
       Rep:=  List Term

       TermS ==>  Record(k : SingleInteger, c : R)
       RepS ==> List TermS

    --declarations
       x,y,p,p1,p2: %
       n: Integer
       nn: NonNegativeInteger
       np: PositiveInteger
       e: E
       r: R

       commutative := R has CommutativeRing

    --local operations
       1  == [[0$E,1$R]]
       characteristic  == characteristic$R
       numberOfMonomials x  == (# x)$Rep
       degree p == if null p then 0 else p.first.k
       minimumDegree p == if null p then 0 else (last p).k
       leadingCoefficient p == if null p then 0$R else p.first.c
       leadingMonomial p == if null p then 0 else [p.first]
       reductum p == if null p then p else p.rest
       retractIfCan(p:%):Union(R,"failed") ==
         null p => 0$R
         not null p.rest => "failed"
         zero?(p.first.k) => p.first.c
         "failed"

       if E is NonNegativeInteger then
           coefficient(p : %, e : E)  ==
               null(p) => 0$R
               degp := first(p).k pretend Integer
               degp < max()$SingleInteger =>
                   not(e pretend Integer < max()$SingleInteger) => 0$R
                   ps := p pretend RepS
                   es := e pretend SingleInteger
                   for tms in ps repeat
                       ks := tms.k pretend SingleInteger
                       ks = es => return tms.c
                       ks < es => return 0$R
                   return 0$R
               for tm in p repeat
                   tm.k=e => return tm.c
                   tm.k < e => return 0$R
               0$R
       else
           coefficient(p : %, e : E)  ==
               for tm in p repeat
                   tm.k=e => return tm.c
                   tm.k < e => return 0$R
               0$R

       recip(p) ==
           null p => "failed"
           p.first.k > 0$E => "failed"
           (u:=recip(p.first.c)) case "failed" => "failed"
           (u::R)::%

       coerce(r) == if zero? r then 0$% else [[0$E,r]]
       coerce(n) == (n::R)::%

       ground?(p): Boolean == empty? p or (empty? rest p and zero? degree p)

       -- qsetrest!: (Rep, Rep) -> Rep
       -- qsetrest!(l: Rep, e: Rep): Rep == RPLACD(l, e)$Lisp
       qsetrest! ==> RPLACD$Lisp

       times!: (R,    %) -> %
       times:  (R, E, %) -> %

       entireRing? := R has EntireRing

       times!(r: R, x: %): % ==
         res, endcell, newend, xx: Rep
         if entireRing? then
                for tx in x repeat tx.c := r*tx.c
                x
         else
                xx := x
                res := empty()
                while not empty? xx repeat
                        tx := first xx
                        tx.c := r * tx.c
                        if zero? tx.c then
                                xx := rest xx
                        else
                                newend := xx
                                xx := rest xx
                                if empty? res then
                                        res := newend
                                        endcell := res
                                else
                                        qsetrest!(endcell, newend)
                                        endcell := newend
                res;

        --- term * polynomial
       termTimes: (R, E, Term) -> Term
       termTimes(r: R, e: E, tx:Term): Term == [e+tx.k, r*tx.c]
       times(tco: R, tex: E, rx: %): % ==
        if entireRing? then
                map(x1+->termTimes(tco, tex, x1), rx::Rep)
        else
                [[tex + tx.k, r] for tx in rx::Rep | not zero? (r := tco * tx.c)]



       -- local addm!
       addm!: (Rep, R, E, Rep) -> Rep
        -- p1 + coef*x^E * p2
        -- `spare' (commented out) is for storage efficiency (not so good for
        -- performance though.
       ADDM_BODY ==>
                res     := empty()
                endcell := empty()
                --spare   := empty()
                while not empty? p1 and not empty? p2 repeat
                        tx := first p1
                        ty := first p2
                        exy := exp + ty.k
                        newcell := empty();
                        if tx.k = exy then
                                newcoef := tx.c + coef * ty.c
                                if not zero? newcoef then
                                        tx.c    := newcoef
                                        newcell := p1
                                --else
                                --      spare   := cons(p1, spare)
                                p1 := rest p1
                                p2 := rest p2
                        else if tx.k > exy then
                                newcell := p1
                                p1      := rest p1
                        else
                                newcoef := coef * ty.c
                                if not entireRing? and zero? newcoef then
                                        newcell := empty()
                                --else if empty? spare then
                                --      ttt := [exy, newcoef]
                                --      newcell := cons(ttt, empty())
                                --else
                                --      newcell := first spare
                                --      spare   := rest spare
                                --      ttt := first newcell
                                --      ttt.k := exy
                                --      ttt.c := newcoef
                                else
                                        ttt := [exy, newcoef]
                                        newcell := cons(ttt, empty())
                                p2 := rest p2
                        if not empty? newcell then
                                if empty? res then
                                        res := newcell
                                        endcell := res
                                else
                                        qsetrest!(endcell, newcell)
                                        endcell := newcell

       addm!(p1:Rep, coef:R, exp: E, p2:Rep): Rep ==
           res, newcell, endcell: Rep
           spare: List Rep
           ADDM_BODY
           if not empty? p1 then  -- then end is const * p1
               newcell := p1
           else  -- then end is (coef, exp) * p2
               newcell := times(coef, exp, p2)
           empty? res => newcell
           qsetrest!(endcell, newcell)
           res

       pomopo! (p1, r, e, p2) ==  addm!(p1, r, e, p2)

       if E is NonNegativeInteger then
           -- addms! works only for SingleInteger exponents
           -- but is much faster than general case
           addms!(pp1 : Rep, coef:R, exp1 : E, pp2 :Rep): Rep ==
               res, newcell, endcell: RepS
               spare: List RepS
               p1 : RepS := pp1 pretend RepS
               p2 : RepS := pp2 pretend RepS
               exp : SingleInteger := exp1 pretend SingleInteger
               ADDM_BODY
               if not empty? p1 then  -- then end is const * p1
                   newcell := p1
               else  -- then end is (coef, exp) * p2
                   newcell := times(coef, exp pretend NonNegativeInteger,
                                    p2 pretend Rep) pretend RepS
               empty? res => newcell pretend Rep
               qsetrest!(endcell, newcell)
               res pretend Rep

           p1 * p2 ==
                xx := p1::Rep
                empty? xx => p1
                yy := p2::Rep
                empty? yy => p2
                if commutative then
                    lx := # xx
                    ly := # yy
                    if lx > ly then
                        (xx, yy) := (yy, xx)
                        (p1, p2) := (p2, p1)
                degx := first(xx).k
                degy := first(yy).k
                zero? degx => first(xx).c * p2
                zero? degy => p1 * first(yy).c
                xx := reverse xx
                res : Rep := empty()
                if degx + degy < max()$SingleInteger then
                    for tx in xx repeat res := addms!(res, tx.c, tx.k, yy)
                else
                    for tx in xx repeat res := addm!(res, tx.c, tx.k, yy)
                res
       else
           p1 * p2 ==
                xx := p1::Rep
                empty? xx => p1
                yy := p2::Rep
                empty? yy => p2
                zero? first(xx).k => first(xx).c * p2
                zero? first(yy).k => p1 * first(yy).c
                --if #xx > #yy then
                --      (xx, yy) := (yy, xx)
                --      (p1, p2) := (p2, p1)
                xx := reverse xx
                res : Rep := empty()
                for tx in xx repeat res:=addm!(res,tx.c,tx.k,yy)
                res

--     if R has EntireRing then
--         p1 * p2  ==
--            null p1 => 0
--            null p2 => 0
--            zero?(p1.first.k) => p1.first.c * p2
--            one? p2 => p1
--            +/[[[t1.k+t2.k,t1.c*t2.c]$Term for t2 in p2]
--                   for t1 in reverse(p1)]
--                   -- This 'reverse' is an efficiency improvement:
--                   -- reduces both time and space [Abbott/Bradford/Davenport]
--        else
--         p1 * p2  ==
--            null p1 => 0
--            null p2 => 0
--            zero?(p1.first.k) => p1.first.c * p2
--            one? p2 => p1
--            +/[[[t1.k+t2.k,r]$Term for t2 in p2 | (r:=t1.c*t2.c) ~= 0]
--                 for t1 in reverse(p1)]
--                  -- This 'reverse' is an efficiency improvement:
--                  -- reduces both time and space [Abbott/Bradford/Davenport]
       if R has CommutativeRing  then
         p ^ np == p ^ (np pretend NonNegativeInteger)
         p ^ np  == p ^ (np pretend NonNegativeInteger)
         p ^ nn  == p ^ nn


         p ^ nn  ==
            null p => 0
            zero? nn => 1
--            one? nn => p
            (nn = 1) => p
            empty? p.rest =>
              zero?(cc:=p.first.c ^ nn) => 0
              [[nn * p.first.k, cc]]
            binomThmExpt([p.first], p.rest, nn)

       if R has Field then
         unitNormal(p) ==
            null p or (lcf:R:=p.first.c) = 1 => [1,p,1]
            a := inv lcf
            [lcf::%, [[p.first.k,1],:(a * p.rest)], a::%]
         unitCanonical(p) ==
            null p or (lcf:R:=p.first.c) = 1 => p
            a := inv lcf
            [[p.first.k,1],:(a * p.rest)]
       else if R has IntegralDomain then
         unitNormal(p) ==
            null p or p.first.c = 1 => [1,p,1]
            (u,cf,a):=unitNormal(p.first.c)
            [u::%, [[p.first.k,cf],:(a * p.rest)], a::%]
         unitCanonical(p) ==
            null p or p.first.c = 1 => p
            (u,cf,a):=unitNormal(p.first.c)
            [[p.first.k,cf],:(a * p.rest)]
       if R has IntegralDomain then
         associates?(p1,p2) ==
            null p1 => null p2
            null p2 => false
            p1.first.k = p2.first.k and
              associates?(p1.first.c,p2.first.c) and
               ((p2.first.c exquo p1.first.c)::R * p1.rest = p2.rest)
         p exquo r  ==
           [(if (a:= tm.c exquo r) case "failed"
               then return "failed" else [tm.k,a])
                  for tm in p] :: Union(%,"failed")
         if E has CancellationAbelianMonoid then
           fmecg(p1:%,e:E,r:R,p2:%):% ==       -- p1 - r * X^e * p2
              rout:%:= []
              r:= - r
              for tm in p2 repeat
                 e2:= e + tm.k
                 c2:= r * tm.c
                 c2 = 0 => "next term"
                 while not null p1 and p1.first.k > e2 repeat
                   (rout:=[p1.first,:rout]; p1:=p1.rest)  --use PUSH and POP?
                 null p1 or p1.first.k < e2 => rout:=[[e2,c2],:rout]
                 if (u:=p1.first.c + c2) ~= 0 then rout:=[[e2, u],:rout]
                 p1:=p1.rest
              NRECONC(rout,p1)$Lisp
           if R has Approximate then
             p1 exquo p2  ==
               null p2 => error "Division by 0"
               p2 = 1 => p1
               p1=p2 => 1
             --(p1.lastElt.c exquo p2.lastElt.c) case "failed" => "failed"
               rout:= []@List(Term)
               while not null p1 repeat
                  (a:= p1.first.c exquo p2.first.c)
                  a case "failed" => return "failed"
                  ee:= subtractIfCan(p1.first.k, p2.first.k)
                  ee case "failed" => return "failed"
                  p1:= fmecg(p1.rest, ee, a, p2.rest)
                  rout:= [[ee,a], :rout]
               null p1 => reverse(rout)::%    -- nreverse?
               "failed"
           else -- R not Approximate
             p1 exquo p2  ==
               null p2 => error "Division by 0"
               p2 = 1 => p1
             --(p1.lastElt.c exquo p2.lastElt.c) case "failed" => "failed"
               rout:= []@List(Term)
               while not null p1 repeat
                  (a:= p1.first.c exquo p2.first.c)
                  a case "failed" => return "failed"
                  ee:= subtractIfCan(p1.first.k, p2.first.k)
                  ee case "failed" => return "failed"
                  p1:= fmecg(p1.rest, ee, a, p2.rest)
                  rout:= [[ee,a], :rout]
               null p1 => reverse(rout)::%    -- nreverse?
               "failed"
       if R has Field then
          x/r == inv(r)*x

\nwendcode{}\nwbegindocs{8}\nwdocspar
\section{domain SUP SparseUnivariatePolynomial}
\nwenddocs{}\nwbegincode{9}\moddef{domain SUP SparseUnivariatePolynomial}\endmoddef\nwstartdeflinemarkup\nwenddeflinemarkup
)abbrev domain SUP SparseUnivariatePolynomial
++ Author: Dave Barton, Barry Trager
++ Date Created:
++ Date Last Updated:
++ Basic Functions: Ring, monomial, coefficient, reductum, differentiate,
++ elt, map, resultant, discriminant
++ Related Constructors: UnivariatePolynomial, Polynomial
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++ This domain represents univariate polynomials over arbitrary
++ (not necessarily commutative) coefficient rings. The variable is
++ unspecified  so that the variable displays as \\spad\{?\} on output.
++ If it is necessary to specify the variable name, use type \\spadtype\{UnivariatePolynomial\}.
++ The representation is sparse
++ in the sense that only non-zero terms are represented.
++ Note: if the coefficient ring is a field, this domain forms a euclidean domain.

SparseUnivariatePolynomial(R:Ring): UnivariatePolynomialCategory(R) with
     outputForm : (%,OutputForm) -> OutputForm
        ++ outputForm(p,var) converts the SparseUnivariatePolynomial p to
        ++ an output form (see \\spadtype\{OutputForm\}) printed as a polynomial in the
        ++ output form variable.
     fmecg: (%,NonNegativeInteger,R,%) -> %
        ++ fmecg(p1,e,r,p2) finds X : p1 - r * X^e * p2
    == PolynomialRing(R,NonNegativeInteger)
  add
   --representations
   Term ==> Record(k:NonNegativeInteger,c:R)
   Rep  := List Term
   p:%
   n:NonNegativeInteger
   np: PositiveInteger
   FP ==> SparseUnivariatePolynomial %
   pp,qq: FP
   lpp:List FP

   -- for karatsuba
   kBound: NonNegativeInteger := 63
   upmp := UnivariatePolynomialMultiplicationPackage(R,%)

   makeSUP(p : %) : SparseUnivariatePolynomial(R) ==
       p pretend SparseUnivariatePolynomial(R)

   unmakeSUP(sp : SparseUnivariatePolynomial(R)) : % ==
       sp pretend %

   if R has FieldOfPrimeCharacteristic  then
         p ^ np == p ^ (np pretend NonNegativeInteger)
         p ^ n  ==
            null p => 0
            zero? n => 1
--            one? n => p
            (n = 1) => p
            empty? p.rest =>
              zero?(cc:=p.first.c ^ n) => 0
              [[n * p.first.k, cc]]
            -- not worth doing special trick if characteristic is too small
            if characteristic()$R < 3 then return expt(p,n pretend PositiveInteger)$RepeatedSquaring(%)
            y:%:=1
            -- break up exponent in qn * characteristic + rn
            -- exponentiating by the characteristic is fast
            rec := divide(n, characteristic()$R)
            qn:= rec.quotient
            rn:= rec.remainder
            repeat
                if rn = 1 then y := y * p
                if rn > 1 then y:= y * binomThmExpt([p.first], p.rest, rn)
                zero? qn => return y
                -- raise to the characteristic power
                p:= [[t.k * characteristic()$R , primeFrobenius(t.c)$R ]$Term for t in p]
                rec := divide(qn, characteristic()$R)
                qn:= rec.quotient
                rn:= rec.remainder
            y



   zero?(p): Boolean == empty?(p)
--   one?(p):Boolean == not empty? p and (empty? rest p and zero? first(p).k and one? first(p).c)
   one?(p):Boolean == not empty? p and (empty? rest p and zero? first(p).k and (first(p).c = 1))
   ground?(p): Boolean == empty? p or (empty? rest p and zero? first(p).k)
   multiplyExponents(p,n) == [ [u.k*n,u.c] for u in p]
   divideExponents(p,n) ==
     null p => p
     m:= (p.first.k :: Integer exquo n::Integer)
     m case "failed" => "failed"
     u:= divideExponents(p.rest,n)
     u case "failed" => "failed"
     [[m::Integer::NonNegativeInteger,p.first.c],:u]
   karatsubaDivide(p, n)  ==
     zero? n => [p, 0]
     lowp: Rep := p
     highp: Rep := []
     repeat
       if empty? lowp then break
       t := first lowp
       if t.k < n then break
       lowp := rest lowp
       highp := cons([subtractIfCan(t.k,n)::NonNegativeInteger,t.c]$Term,highp)
     [ reverse highp,  lowp]
   shiftRight(p, n)  ==
      [[subtractIfCan(t.k,n)::NonNegativeInteger,t.c]$Term for t in p]
   shiftLeft(p, n)  ==
      [[t.k + n,t.c]$Term for t in p]
   pomopo!(p1,r,e,p2) ==
          rout:%:= []
          for tm in p2 repeat
             e2:= e + tm.k
             c2:= r * tm.c
             c2 = 0 => "next term"
             while not null p1 and p1.first.k > e2 repeat
               (rout:=[p1.first,:rout]; p1:=p1.rest)  --use PUSH and POP?
             null p1 or p1.first.k < e2 => rout:=[[e2,c2],:rout]
             if (u:=p1.first.c + c2) ~= 0 then rout:=[[e2, u],:rout]
             p1:=p1.rest
          NRECONC(rout,p1)$Lisp

-- implementation using karatsuba algorithm conditionally
--
--   p1 * p2 ==
--     xx := p1::Rep
--     empty? xx => p1
--     yy := p2::Rep
--     empty? yy => p2
--     zero? first(xx).k => first(xx).c * p2
--     zero? first(yy).k => p1 * first(yy).c
--     (first(xx).k > kBound) and (first(yy).k > kBound) and (#xx > kBound) and (#yy > kBound) =>
--       karatsubaOnce(p1,p2)$upmp
--     xx := reverse xx
--     res : Rep := empty()
--     for tx in xx repeat res:= rep pomopo!( res,tx.c,tx.k,p2)
--     res


   univariate(p:%) == p pretend SparseUnivariatePolynomial(R)
   multivariate(sup:SparseUnivariatePolynomial(R),v:SingletonAsOrderedSet) ==
      sup pretend %
   univariate(p:%,v:SingletonAsOrderedSet) ==
     zero? p => 0
     monomial(leadingCoefficient(p)::%,degree p) +
         univariate(reductum p,v)
   multivariate(supp:SparseUnivariatePolynomial(%),v:SingletonAsOrderedSet) ==
     zero? supp => 0
     lc:=leadingCoefficient supp
     degree lc > 0 => error "bad form polynomial"
     monomial(leadingCoefficient lc,degree supp) +
         multivariate(reductum supp,v)
   if R has FiniteFieldCategory and R has PolynomialFactorizationExplicit then
     RXY ==> SparseUnivariatePolynomial SparseUnivariatePolynomial R
     squareFreePolynomial pp ==
        squareFree(pp)$UnivariatePolynomialSquareFree(%,FP)
     factorPolynomial pp ==
          (generalTwoFactor(pp pretend RXY)$TwoFactorize(R))
                      pretend Factored SparseUnivariatePolynomial %
     factorSquareFreePolynomial pp ==
          (generalTwoFactor(pp pretend RXY)$TwoFactorize(R))
                      pretend Factored SparseUnivariatePolynomial %
     gcdPolynomial(pp,qq) == gcd(pp,qq)$FP
     factor p == factor(p)$DistinctDegreeFactorize(R,%)
     solveLinearPolynomialEquation(lpp,pp) ==
       solveLinearPolynomialEquation(lpp, pp)$FiniteFieldSolveLinearPolynomialEquation(R,%,FP)
   else if R has PolynomialFactorizationExplicit then
     import PolynomialFactorizationByRecursionUnivariate(R,%)
     solveLinearPolynomialEquation(lpp,pp)==
       solveLinearPolynomialEquationByRecursion(lpp,pp)
     factorPolynomial(pp) ==
       factorByRecursion(pp)
     factorSquareFreePolynomial(pp) ==
       factorSquareFreeByRecursion(pp)

   if R has IntegralDomain then
    if R has Approximate then
     p1 exquo p2  ==
        null p2 => error "Division by 0"
        p2 = 1 => p1
        p1=p2 => 1
      --(p1.lastElt.c exquo p2.lastElt.c) case "failed" => "failed"
        rout:= []@List(Term)
        while not null p1 repeat
           (a:= p1.first.c exquo p2.first.c)
           a case "failed" => return "failed"
           ee:= subtractIfCan(p1.first.k, p2.first.k)
           ee case "failed" => return "failed"
           p1:= fmecg(p1.rest, ee, a, p2.rest)
           rout:= [[ee,a], :rout]
        null p1 => reverse(rout)::%    -- nreverse?
        "failed"
    else -- R not Approximate
     p1 exquo p2  ==
        null p2 => error "Division by 0"
        p2 = 1 => p1
      --(p1.lastElt.c exquo p2.lastElt.c) case "failed" => "failed"
        rout:= []@List(Term)
        while not null p1 repeat
           (a:= p1.first.c exquo p2.first.c)
           a case "failed" => return "failed"
           ee:= subtractIfCan(p1.first.k, p2.first.k)
           ee case "failed" => return "failed"
           p1:= fmecg(p1.rest, ee, a, p2.rest)
           rout:= [[ee,a], :rout]
        null p1 => reverse(rout)::%    -- nreverse?
        "failed"
   fmecg(p1,e,r,p2) ==       -- p1 - r * X^e * p2
          rout:%:= []
          r:= - r
          for tm in p2 repeat
             e2:= e + tm.k
             c2:= r * tm.c
             c2 = 0 => "next term"
             while not null p1 and p1.first.k > e2 repeat
               (rout:=[p1.first,:rout]; p1:=p1.rest)  --use PUSH and POP?
             null p1 or p1.first.k < e2 => rout:=[[e2,c2],:rout]
             if (u:=p1.first.c + c2) ~= 0 then rout:=[[e2, u],:rout]
             p1:=p1.rest
          NRECONC(rout,p1)$Lisp
   pseudoRemainder(p1,p2) ==
     null p2 => error "PseudoDivision by Zero"
     null p1 => 0
     co:=p2.first.c;
     e:=p2.first.k;
     p2:=p2.rest;
     e1:=max((p1.first.k)::Integer-e+1,0)::NonNegativeInteger
     while not null p1 repeat
       if (u:=subtractIfCan(p1.first.k,e)) case "failed" then leave
       p1:=fmecg(co * p1.rest, u, p1.first.c, p2)
       e1:= (e1 - 1)::NonNegativeInteger
     e1 = 0 => p1
     co ^ e1 * p1
   toutput(t1:Term,v:OutputForm):OutputForm ==
     t1.k = 0 => t1.c :: OutputForm
     if t1.k = 1
       then mon:= v
       else mon := v ^ t1.k::OutputForm
     t1.c = 1 => mon
     t1.c = -1 and
          ((t1.c :: OutputForm) = (-1$Integer)::OutputForm)@Boolean => - mon
     t1.c::OutputForm * mon
   outputForm(p:%,v:OutputForm) ==
     l: List(OutputForm)
     l:=[toutput(t,v) for t in p]
     null l => (0$Integer)::OutputForm -- else FreeModule 0 problems
     reduce("+",l)

   coerce(p:%):OutputForm == outputForm(p, "?"::OutputForm)
   elt(p:%,val:R) ==
      null p => 0$R
      co:=p.first.c
      n:=p.first.k
      for tm in p.rest repeat
       co:= co * val ^ (n - (n:=tm.k))::NonNegativeInteger + tm.c
      n = 0 => co
      co * val ^ n
   elt(p:%,val:%) ==
      null p => 0$%
      coef:% := p.first.c :: %
      n:=p.first.k
      for tm in p.rest repeat
       coef:= coef * val ^ (n-(n:=tm.k))::NonNegativeInteger+(tm.c::%)
      n = 0 => coef
      coef * val ^ n

   monicDivide(p1:%,p2:%) ==
      null p2 => error "monicDivide: division by 0"
      leadingCoefficient p2 ~= 1 => error "Divisor Not Monic"
      p2 = 1 => [p1,0]
      null p1 => [0,0]
      degree p1 < (n:=degree p2) => [0,p1]
      rout:Rep := []
      p2 := p2.rest
      while not null p1 repeat
         (u:=subtractIfCan(p1.first.k, n)) case "failed" => leave
         rout:=[[u, p1.first.c], :rout]
         p1:=fmecg(p1.rest, rout.first.k, rout.first.c, p2)
      [reverse_!(rout),p1]

   if R has IntegralDomain then
       discriminant(p) == discriminant(p)$PseudoRemainderSequence(R,%)
--     discriminant(p) ==
--       null p or zero?(p.first.k) => error "cannot take discriminant of constants"
--       dp:=differentiate p
--       corr:=  p.first.c ^ ((degree p - 1 - degree dp)::NonNegativeInteger)
--       (-1)^((p.first.k*(p.first.k-1)) quo 2):NonNegativeInteger
--         * (corr * resultant(p,dp) exquo p.first.c)::R

       subResultantGcd(p1,p2) == subResultantGcd(p1,p2)$PseudoRemainderSequence(R,%)
--     subResultantGcd(p1,p2) ==    --args # 0, non-coef, prim, ans not prim
--       --see algorithm 1 (p. 4) of Brown's latest (unpublished) paper
--       if p1.first.k < p2.first.k then (p1,p2):=(p2,p1)
--       p:=pseudoRemainder(p1,p2)
--       co:=1$R;
--       e:= (p1.first.k - p2.first.k):NonNegativeInteger
--       while not null p and p.first.k ~= 0 repeat
--         p1:=p2; p2:=p; p:=pseudoRemainder(p1,p2)
--         null p or p.first.k = 0 => "enuf"
--         co:=(p1.first.c ^ e exquo co ^ max(0, (e-1))::NonNegativeInteger)::R
--         e:= (p1.first.k - p2.first.k):NonNegativeInteger;  c1:=co^e
--         p:=[[tm.k,((tm.c exquo p1.first.c)::R exquo c1)::R] for tm in p]
--       if null p then p2 else 1$%

       resultant(p1,p2) == resultant(p1,p2)$PseudoRemainderSequence(R,%)
--     resultant(p1,p2) ==      --SubResultant PRS Algorithm
--        null p1 or null p2 => 0$R
--        0 = degree(p1) => ((first p1).c)^degree(p2)
--        0 = degree(p2) => ((first p2).c)^degree(p1)
--        if p1.first.k < p2.first.k then
--          (if odd?(p1.first.k) then p1:=-p1;  (p1,p2):=(p2,p1))
--        p:=pseudoRemainder(p1,p2)
--        co:=1$R;  e:=(p1.first.k-p2.first.k):NonNegativeInteger
--        while not null p repeat
--           if not odd?(e) then p:=-p
--           p1:=p2;  p2:=p;  p:=pseudoRemainder(p1,p2)
--           co:=(p1.first.c ^ e exquo co ^ max(e:Integer-1,0):NonNegativeInteger)::R
--           e:= (p1.first.k - p2.first.k):NonNegativeInteger;  c1:=co^e
--           p:=(p exquo ((leadingCoefficient p1) * c1))::%
--        degree p2 > 0 => 0$R
--        (p2.first.c^e exquo co^((e-1)::NonNegativeInteger))::R
   if R has GcdDomain then
     content(p) == if null p then 0$R else "gcd"/[tm.c for tm in p]
        --make CONTENT more efficient?

     primitivePart(p) ==
        null p => p
        ct :=content(p)
        unitCanonical((p exquo ct)::%)
               -- exquo  present since % is now an IntegralDomain

     gcd(p1,p2) ==
          gcdPolynomial(p1 pretend SparseUnivariatePolynomial R,
                        p2 pretend SparseUnivariatePolynomial R) pretend %

   if R has Field then
     divide( p1, p2)  ==
       zero? p2 => error "Division by 0"
--       one? p2 => [p1,0]
       (p2 = 1) => [p1,0]
       ct:=inv(p2.first.c)
       n:=p2.first.k
       p2:=p2.rest
       rout:=empty()$List(Term)
       while p1 ~= 0 repeat
          (u:=subtractIfCan(p1.first.k, n)) case "failed" => leave
          rout:=[[u, ct * p1.first.c], :rout]
          p1:=fmecg(p1.rest, rout.first.k, rout.first.c, p2)
       [reverse_!(rout),p1]

     p / co == inv(co) * p

\nwendcode{}\nwbegindocs{10}\nwdocspar
\section{package SUP2 SparseUnivariatePolynomialFunctions2}
\nwenddocs{}\nwbegincode{11}\moddef{package SUP2 SparseUnivariatePolynomialFunctions2}\endmoddef\nwstartdeflinemarkup\nwenddeflinemarkup
)abbrev package SUP2 SparseUnivariatePolynomialFunctions2
++ Author:
++ Date Created:
++ Date Last Updated:
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++ This package lifts a mapping from coefficient rings R to S to
++ a mapping from sparse univariate polynomial over R to
++ a sparse univariate polynomial over S.
++ Note that the mapping is assumed
++ to send zero to zero, since it will only be applied to the non-zero
++ coefficients of the polynomial.

SparseUnivariatePolynomialFunctions2(R:Ring, S:Ring): with
  map:(R->S,SparseUnivariatePolynomial R) -> SparseUnivariatePolynomial S
    ++ map(func, poly) creates a new polynomial by applying func to
    ++ every non-zero coefficient of the polynomial poly.
 == add
  map(f, p) == map(f, p)$UnivariatePolynomialCategoryFunctions2(R,
           SparseUnivariatePolynomial R, S, SparseUnivariatePolynomial S)

\nwendcode{}\nwbegindocs{12}\nwdocspar
\section{domain UP UnivariatePolynomial}
\nwenddocs{}\nwbegincode{13}\moddef{domain UP UnivariatePolynomial}\endmoddef\nwstartdeflinemarkup\nwenddeflinemarkup
)abbrev domain UP UnivariatePolynomial
++ Author:
++ Date Created:
++ Date Last Updated:
++ Basic Functions: Ring, monomial, coefficient, reductum, differentiate,
++ elt, map, resultant, discriminant
++ Related Constructors: SparseUnivariatePolynomial, MultivariatePolynomial
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++ This domain represents univariate polynomials in some symbol
++ over arbitrary (not necessarily commutative) coefficient rings.
++ The representation is sparse
++ in the sense that only non-zero terms are represented.
++ Note: if the coefficient ring is a field, then this domain forms a euclidean domain.

UnivariatePolynomial(x:Symbol, R:Ring):
  UnivariatePolynomialCategory(R) with
    coerce: Variable(x) -> %
      ++ coerce(x) converts the variable x to a univariate polynomial.
    fmecg: (%,NonNegativeInteger,R,%) -> %
        ++ fmecg(p1,e,r,p2) finds X : p1 - r * X^e * p2
   == SparseUnivariatePolynomial(R)   add
    Rep:=SparseUnivariatePolynomial(R)
    coerce(p:%):OutputForm  == outputForm(p, outputForm x)
    coerce(v:Variable(x)):% == monomial(1, 1)

\nwendcode{}\nwbegindocs{14}\nwdocspar
\section{package UP2 UnivariatePolynomialFunctions2}
\nwenddocs{}\nwbegincode{15}\moddef{package UP2 UnivariatePolynomialFunctions2}\endmoddef\nwstartdeflinemarkup\nwenddeflinemarkup
)abbrev package UP2 UnivariatePolynomialFunctions2
++ Author:
++ Date Created:
++ Date Last Updated:
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++ This package lifts a mapping from coefficient rings R to S to
++ a mapping from \\spadtype\{UnivariatePolynomial\}(x,R) to
++ \\spadtype\{UnivariatePolynomial\}(y,S). Note that the mapping is assumed
++ to send zero to zero, since it will only be applied to the non-zero
++ coefficients of the polynomial.

UnivariatePolynomialFunctions2(x:Symbol, R:Ring, y:Symbol, S:Ring): with
  map: (R -> S, UnivariatePolynomial(x,R)) -> UnivariatePolynomial(y,S)
    ++ map(func, poly) creates a new polynomial by applying func to
    ++ every non-zero coefficient of the polynomial poly.
 == add
  map(f, p) == map(f, p)$UnivariatePolynomialCategoryFunctions2(R,
              UnivariatePolynomial(x, R), S, UnivariatePolynomial(y, S))

\nwendcode{}\nwbegindocs{16}\nwdocspar
\section{package POLY2UP PolynomialToUnivariatePolynomial}
\nwenddocs{}\nwbegincode{17}\moddef{package POLY2UP PolynomialToUnivariatePolynomial}\endmoddef\nwstartdeflinemarkup\nwenddeflinemarkup
)abbrev package POLY2UP PolynomialToUnivariatePolynomial
++ Author:
++ Date Created:
++ Date Last Updated:
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++ This package is primarily to help the interpreter do coercions.
++ It allows you to view a polynomial as a
++ univariate polynomial in one of its variables with
++ coefficients which are again a polynomial in all the
++ other variables.

PolynomialToUnivariatePolynomial(x:Symbol, R:Ring): with
  univariate: (Polynomial R, Variable x) ->
                                   UnivariatePolynomial(x, Polynomial R)
     ++ univariate(p, x) converts the polynomial p to a one of type
     ++ \\spad\{UnivariatePolynomial(x,Polynomial(R))\}, ie. as a member of \\spad\{R[...][x]\}.
 == add
  univariate(p, y) ==
    q:SparseUnivariatePolynomial(Polynomial R) := univariate(p, x)
    map(x1+->x1, q)$UnivariatePolynomialCategoryFunctions2(Polynomial R,
                  SparseUnivariatePolynomial Polynomial R, Polynomial R,
                      UnivariatePolynomial(x, Polynomial R))

\nwendcode{}\nwbegindocs{18}\nwdocspar
\section{package UPSQFREE UnivariatePolynomialSquareFree}
\nwenddocs{}\nwbegincode{19}\moddef{package UPSQFREE UnivariatePolynomialSquareFree}\endmoddef\nwstartdeflinemarkup\nwenddeflinemarkup
)abbrev package UPSQFREE UnivariatePolynomialSquareFree
++ Author: Dave Barton, Barry Trager
++ Date Created:
++ Date Last Updated:
++ Basic Functions: squareFree, squareFreePart
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++ This package provides for square-free decomposition of
++ univariate polynomials over arbitrary rings, i.e.
++ a partial factorization such that each factor is a product
++ of irreducibles with multiplicity one and the factors are
++ pairwise relatively prime. If the ring
++ has characteristic zero, the result is guaranteed to satisfy
++ this condition. If the ring is an infinite ring of
++ finite characteristic, then it may not be possible to decide when
++ polynomials contain factors which are pth powers. In this
++ case, the flag associated with that polynomial is set to "nil"
++ (meaning that that polynomials are not guaranteed to be square-free).

UnivariatePolynomialSquareFree(RC:IntegralDomain,P):C == T
  where
    fUnion ==> Union("nil", "sqfr", "irred", "prime")
    FF     ==> Record(flg:fUnion, fctr:P, xpnt:Integer)
    P:Join(UnivariatePolynomialCategory(RC),IntegralDomain) with
      gcd: (%,%) -> %
        ++ gcd(p,q) computes the greatest-common-divisor of p and q.

    C == with
      squareFree: P -> Factored(P)
        ++ squareFree(p) computes the square-free factorization of the
        ++ univariate polynomial p. Each factor has no repeated roots, and the
        ++ factors are pairwise relatively prime.
      squareFreePart: P -> P
        ++ squareFreePart(p) returns a polynomial which has the same
        ++ irreducible factors as the univariate polynomial p, but each
        ++ factor has multiplicity one.
      BumInSepFFE: FF -> FF
        ++ BumInSepFFE(f) is a local function, exported only because
        ++ it has multiple conditional definitions.

    T == add

      if RC has CharacteristicZero then
        squareFreePart(p:P) == (p exquo gcd(p, differentiate p))::P
      else
        squareFreePart(p:P) ==
          unit(s := squareFree(p)$%) * */[f.factor for f in factors s]

      if RC has FiniteFieldCategory then
        BumInSepFFE(ffe:FF) ==
           ["sqfr", map(charthRoot,ffe.fctr), characteristic$P*ffe.xpnt]
      else if RC has CharacteristicNonZero then
         BumInSepFFE(ffe:FF) ==
            np := multiplyExponents(ffe.fctr,characteristic$P::NonNegativeInteger)
            (nthrp := charthRoot(np)) case "failed" =>
               ["nil", np, ffe.xpnt]
            ["sqfr", nthrp, characteristic$P*ffe.xpnt]

      else
        BumInSepFFE(ffe:FF) ==
          ["nil",
           multiplyExponents(ffe.fctr,characteristic$P::NonNegativeInteger),
            ffe.xpnt]


      if RC has CharacteristicZero then
        squareFree(p:P) ==             --Yun's algorithm - see SYMSAC '76, p.27
           --Note ci primitive is, so GCD's don't need to %do contents.
           --Change gcd to return cofctrs also?
           ci:=p; di:=differentiate(p); pi:=gcd(ci,di)
           degree(pi)=0 =>
             (u,c,a):=unitNormal(p)
             makeFR(u,[["sqfr",c,1]])
           i:NonNegativeInteger:=0; lffe:List FF:=[]
           lcp := leadingCoefficient p
           while degree(ci) ~= 0 repeat
              ci:=(ci exquo pi)::P
              di:=(di exquo pi)::P - differentiate(ci)
              pi:=gcd(ci,di)
              i:=i+1
              degree(pi) > 0 =>
                 lcp:=(lcp exquo (leadingCoefficient(pi)^i))::RC
                 lffe:=[["sqfr",pi,i],:lffe]
           makeFR(lcp::P,lffe)

      else
        squareFree(p:P) ==           --Musser's algorithm - see SYMSAC '76, p.27
             --p MUST BE PRIMITIVE, Any characteristic.
             --Note ci primitive, so GCD's don't need to %do contents.
             --Change gcd to return cofctrs also?
           ci := gcd(p,differentiate(p))
           degree(ci)=0 =>
             (u,c,a):=unitNormal(p)
             makeFR(u,[["sqfr",c,1]])
           di := (p exquo ci)::P
           i:NonNegativeInteger:=0; lffe:List FF:=[]
           dunit : P := 1
           while degree(di) ~= 0 repeat
              diprev := di
              di := gcd(ci,di)
              ci:=(ci exquo di)::P
              i:=i+1
              degree(diprev) = degree(di) =>
                 lc := (leadingCoefficient(diprev) exquo leadingCoefficient(di))::RC
                 dunit := lc^i * dunit
              pi:=(diprev exquo di)::P
              lffe:=[["sqfr",pi,i],:lffe]
           dunit := dunit * di ^ (i+1)
           degree(ci)=0 => makeFR(dunit*ci,lffe)
           redSqfr:=squareFree(divideExponents(ci,characteristic$P)::P)
           lsnil:= [BumInSepFFE(ffe) for ffe in factorList redSqfr]
           lffe:=append(lsnil,lffe)
           makeFR(dunit*(unit redSqfr),lffe)

\nwendcode{}\nwbegindocs{20}\nwdocspar
\section{package PSQFR PolynomialSquareFree}
\nwenddocs{}\nwbegincode{21}\moddef{package PSQFR PolynomialSquareFree}\endmoddef\nwstartdeflinemarkup\nwenddeflinemarkup
)abbrev package PSQFR PolynomialSquareFree
++ Author:
++ Date Created:
++ Date Last Updated: November 1993, (P.Gianni)
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++ This package computes square-free decomposition of multivariate
++ polynomials over a coefficient ring which is an arbitrary gcd domain.
++ The requirement on the coefficient domain guarantees that the \\spadfun\{content\} can be
++ removed so that factors will be primitive as well as square-free.
++ Over an infinite ring of finite characteristic,it may not be possible to
++ guarantee that the factors are square-free.

PolynomialSquareFree(VarSet:OrderedSet,E,RC:GcdDomain,P):C == T where
  E:OrderedAbelianMonoidSup
  P:PolynomialCategory(RC,E,VarSet)

  C == with
    squareFree : P -> Factored P
      ++ squareFree(p) returns the square-free factorization of the
      ++ polynomial p.  Each factor has no repeated roots, and the
      ++ factors are pairwise relatively prime.

  T == add
    SUP    ==> SparseUnivariatePolynomial(P)
    NNI    ==> NonNegativeInteger
    fUnion ==> Union("nil", "sqfr", "irred", "prime")
    FF     ==> Record(flg:fUnion, fctr:P, xpnt:Integer)

    finSqFr : (P,List VarSet) -> Factored P
    pthPower : P -> Factored P
    pPolRoot : P -> P
    putPth   : P -> P

    chrc:=characteristic$RC

    if RC has CharacteristicNonZero then
    -- find the p-th root of a polynomial
      pPolRoot(f:P) : P ==
        lvar:=variables f
        empty? lvar => f
        mv:=first lvar
        uf:=univariate(f,mv)
        uf:=divideExponents(uf,chrc)::SUP
        uf:=map(pPolRoot,uf)
        multivariate(uf,mv)

    -- substitute variables with their p-th power
      putPth(f:P) : P ==
        lvar:=variables f
        empty? lvar => f
        mv:=first lvar
        uf:=univariate(f,mv)
        uf:=multiplyExponents(uf,chrc)::SUP
        uf:=map(putPth,uf)
        multivariate(uf,mv)

    -- the polynomial is a perfect power
      pthPower(f:P) : Factored P ==
        proot : P := 0
        isSq  : Boolean := false
        if (g:=charthRoot f) case "failed" then proot:=pPolRoot(f)
        else
          proot := g :: P
          isSq  := true
        psqfr:=finSqFr(proot,variables f)
        isSq  =>
          makeFR((unit psqfr)^chrc,[[u.flg,u.fctr,
           (u.xpnt)*chrc] for u in factorList psqfr])
        makeFR((unit psqfr),[["nil",putPth u.fctr,u.xpnt]
                             for u in factorList psqfr])

    -- compute the square free decomposition, finite characteristic case
      finSqFr(f:P,lvar:List VarSet) : Factored P ==
         empty? lvar => pthPower(f)
         mv:=first lvar
         lvar:=lvar.rest
         differentiate(f,mv)=0 => finSqFr(f,lvar)
         uf:=univariate(f,mv)
         cont := content uf
         cont1:P:=1
         uf := (uf exquo cont)::SUP
         squf := squareFree(uf)$UnivariatePolynomialSquareFree(P,SUP)
         pfaclist:List FF :=[]
         for u in factorList squf repeat
           uexp:NNI:=(u.xpnt)::NNI
           u.flg = "sqfr" =>  -- the square free factor is OK
             pfaclist:= cons([u.flg,multivariate(u.fctr,mv),uexp],
                              pfaclist)
           --listfin1:= finSqFr(multivariate(u.fctr,mv),lvar)
           listfin1:= squareFree multivariate(u.fctr,mv)
           flistfin1:=[[uu.flg,uu.fctr,uu.xpnt*uexp]
                        for uu in factorList listfin1]
           cont1:=cont1*((unit listfin1)^uexp)
           pfaclist:=append(flistfin1,pfaclist)
         cont:=cont*cont1
         cont ~= 1 =>
           sqp := squareFree cont
           pfaclist:= append (factorList sqp,pfaclist)
           makeFR(unit(sqp)*coefficient(unit squf,0),pfaclist)
         makeFR(coefficient(unit squf,0),pfaclist)

    squareFree(p:P) ==
       mv:=mainVariable p
       mv case "failed" => makeFR(p,[])$Factored(P)
       characteristic$RC ~= 0 => finSqFr(p,variables p)
       up:=univariate(p,mv)
       cont := content up
       up := (up exquo cont)::SUP
       squp := squareFree(up)$UnivariatePolynomialSquareFree(P,SUP)
       pfaclist:List FF :=
         [[u.flg,multivariate(u.fctr,mv),u.xpnt]
                                            for u in factorList squp]
       cont ~= 1 =>
         sqp := squareFree cont
         makeFR(unit(sqp)*coefficient(unit squp,0),
              append(factorList sqp, pfaclist))
       makeFR(coefficient(unit squp,0),pfaclist)

\nwendcode{}\nwbegindocs{22}\nwdocspar
\section{package UPMP UnivariatePolynomialMultiplicationPackage}
\nwenddocs{}\nwbegincode{23}\moddef{package UPMP UnivariatePolynomialMultiplicationPackage}\endmoddef\nwstartdeflinemarkup\nwenddeflinemarkup
)abbrev package UPMP UnivariatePolynomialMultiplicationPackage
++ Author: Marc Moreno Maza
++ Date Created: 14.08.2000
++ Description:
++ This package implements Karatsuba's trick for multiplying
++ (large) univariate polynomials. It could be improved with
++ a version doing the work on place and also with a special
++ case for squares. We've done this in Basicmath, but we
++ believe that this out of the scope of FriCAS.
-- Currently unused.

UnivariatePolynomialMultiplicationPackage(R: Ring, U: UnivariatePolynomialCategory(R)): C == T
  where
    HL ==> Record(quotient:U,remainder:U)
    C == with
      noKaratsuba: (U, U) -> U
        ++ \\spad\{noKaratsuba(a,b)\} returns \\spad\{a*b\} without
        ++ using Karatsuba's trick at all.
      karatsubaOnce: (U, U) -> U
        ++ \\spad\{karatsuba(a,b)\} returns \\spad\{a*b\} by applying
        ++ Karatsuba's trick once. The other multiplications
        ++ are performed by calling \\spad\{*\} from \\spad\{U\}.
      karatsuba: (U, U, NonNegativeInteger, NonNegativeInteger) -> U
        ++ \\spad\{karatsuba(a,b,l,k)\} returns \\spad\{a*b\} by applying
        ++ Karatsuba's trick provided that both \\spad\{a\} and \\spad\{b\}
        ++ have at least \\spad\{l\} terms and \\spad\{k > 0\} holds
        ++ and by calling \\spad\{noKaratsuba\} otherwise. The other
        ++ multiplications are performed by recursive calls with
        ++ the same third argument and \\spad\{k-1\} as fourth argument.

    T == add
      noKaratsuba(a,b) ==
        zero? a => a
        zero? b => b
        zero?(degree(a)) => leadingCoefficient(a) * b
        zero?(degree(b)) => a * leadingCoefficient(b)
        lu: List(U) := reverse monomials(a)
        res: U := 0;
        for u in lu repeat
          res := pomopo!(res, leadingCoefficient(u), degree(u), b)
        res
      karatsubaOnce(a:U,b:U): U ==
        da := minimumDegree(a)
        db := minimumDegree(b)
        if not zero? da then a := shiftRight(a,da)
        if not zero? db then b := shiftRight(b,db)
        d := da + db
        n: NonNegativeInteger := min(degree(a),degree(b)) quo 2
        rec: HL := karatsubaDivide(a, n)
        ha := rec.quotient
        la := rec.remainder
        rec := karatsubaDivide(b, n)
        hb := rec.quotient
        lb := rec.remainder
        w: U := (ha - la) * (lb - hb)
        u: U := la * lb
        v: U := ha * hb
        w := w + (u + v)
        w := shiftLeft(w,n) + u
        zero? d => shiftLeft(v,2*n) + w
        shiftLeft(v,2*n + d) + shiftLeft(w,d)
      karatsuba(a:U,b:U,l:NonNegativeInteger,k:NonNegativeInteger): U ==
        zero? k => noKaratsuba(a,b)
        degree(a) < l => noKaratsuba(a,b)
        degree(b) < l => noKaratsuba(a,b)
        numberOfMonomials(a) < l => noKaratsuba(a,b)
        numberOfMonomials(b) < l => noKaratsuba(a,b)
        da := minimumDegree(a)
        db := minimumDegree(b)
        if not zero? da then a := shiftRight(a,da)
        if not zero? db then b := shiftRight(b,db)
        d := da + db
        n: NonNegativeInteger := min(degree(a),degree(b)) quo 2
        k := subtractIfCan(k,1)::NonNegativeInteger
        rec: HL := karatsubaDivide(a, n)
        ha := rec.quotient
        la := rec.remainder
        rec := karatsubaDivide(b, n)
        hb := rec.quotient
        lb := rec.remainder
        w: U := karatsuba(ha - la, lb - hb, l, k)
        u: U := karatsuba(la, lb, l, k)
        v: U := karatsuba(ha, hb, l, k)
        w := w + (u + v)
        w := shiftLeft(w,n) + u
        zero? d => shiftLeft(v,2*n) + w
        shiftLeft(v,2*n + d) + shiftLeft(w,d)

\nwendcode{}\nwbegindocs{24}\nwdocspar
\section{License}
\nwenddocs{}\nwbegincode{25}\moddef{license}\endmoddef\nwstartdeflinemarkup\nwenddeflinemarkup
--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.
\nwendcode{}\nwbegindocs{26}\nwdocspar
\nwenddocs{}\nwbegincode{27}\moddef{*}\endmoddef\nwstartdeflinemarkup\nwenddeflinemarkup
\LA{}license\RA{}

\LA{}category FMCAT FreeModuleCategory\RA{}
\LA{}domain FM FreeModule\RA{}
\LA{}package FM2 FreeModuleFunctions2\RA{}
\LA{}domain PR PolynomialRing\RA{}
\LA{}package UPSQFREE UnivariatePolynomialSquareFree\RA{}
\LA{}package PSQFR PolynomialSquareFree\RA{}
\LA{}package UPMP UnivariatePolynomialMultiplicationPackage\RA{}
\LA{}domain SUP SparseUnivariatePolynomial\RA{}
\LA{}package SUP2 SparseUnivariatePolynomialFunctions2\RA{}
\LA{}domain UP UnivariatePolynomial\RA{}
\LA{}package UP2 UnivariatePolynomialFunctions2\RA{}
\LA{}package POLY2UP PolynomialToUnivariatePolynomial\RA{}
\nwendcode{}\nwbegindocs{28}\nwdocspar
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}
\nwenddocs{}
