--Copyright (c) 1991-2002, The Numerical ALgorithms Group Ltd.
--All rights reserved.
--
--Redistribution and use in source and binary forms, with or without
--modification, are permitted provided that the following conditions are
--met:
--
--    - Redistributions of source code must retain the above copyright
--      notice, this list of conditions and the following disclaimer.
--
--    - Redistributions in binary form must reproduce the above copyright
--      notice, this list of conditions and the following disclaimer in
--      the documentation and/or other materials provided with the
--      distribution.
--
--    - Neither the name of The Numerical ALgorithms Group Ltd. nor the
--      names of its contributors may be used to endorse or promote products
--      derived from this software without specific prior written permission.
--
--THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
--IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
--TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
--PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
--OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
--EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
--PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
--PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
--LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
--NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
--SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

)abbrev domain INDE IndexedExponents
++ Author: James Davenport
++ Date Created:
++ Date Last Updated:
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++   IndexedExponents of an ordered set of variables gives a representation
++ for the degree of polynomials in commuting variables. It gives an ordered
++ pairing of non negative integer exponents with variables

IndexedExponents(Varset:OrderedSet): C == T where
  C == Join(OrderedAbelianMonoidSup,
            IndexedDirectProductCategory(NonNegativeInteger,Varset))
  T == IndexedDirectProductOrderedAbelianMonoidSup(NonNegativeInteger,Varset) add
      Term:=  Record(k:Varset,c:NonNegativeInteger)
      Rep:=  List Term
      coerceOF(t: Term):OutputForm ==     -- converts term to OutputForm
         t.c = 1 => (t.k)::OutputForm
         (t.k)::OutputForm ** (t.c)::OutputForm
      coerce(x: %):OutputForm == -- converts entire exponents to OutputForm
         null x => 1::Integer::OutputForm
         null rest x => coerceOF(first x)
         reduce("*",[coerceOF t for t in x])

)abbrev domain SMP SparseMultivariatePolynomial
++ Author: Dave Barton, Barry Trager
++ Date Created:
++ Date Last Updated: 30 November 1994
++ Fix History:
++ 30 Nov 94: added gcdPolynomial for float-type coefficients
++ Basic Functions: Ring, degree, eval, coefficient, monomial, differentiate,
++ resultant, gcd
++ Related Constructors: Polynomial, MultivariatePolynomial
++ Also See:
++ AMS Classifications:
++ Keywords: polynomial, multivariate
++ References:
++ Description:
++   This type is the basic representation of sparse recursive multivariate
++ polynomials. It is parameterized by the coefficient ring and the
++ variable set which may be infinite. The variable ordering is determined
++ by the variable set parameter. The coefficient ring may be non-commutative,
++ but the variables are assumed to commute.

SparseMultivariatePolynomial(R: Ring,VarSet: OrderedSet): C == T where
  pgcd ==> PolynomialGcdPackage(IndexedExponents VarSet,VarSet,R,%)
  C == PolynomialCategory(R,IndexedExponents(VarSet),VarSet)
  SUP ==> SparseUnivariatePolynomial
  T == add
    --constants
    --D := F(%) replaced by next line until compiler support completed

    --representations
      D := SparseUnivariatePolynomial(%)
      VPoly:=  Record(v:VarSet,ts:D)
      Rep:=  Union(R,VPoly)

    --local function


    --declarations
      fn: R -> R
      n: Integer
      k: NonNegativeInteger
      kp:PositiveInteger
      k1:NonNegativeInteger
      c: R
      mvar: VarSet
      val : R
      var:VarSet
      up: D
      p,p1,p2,pval: %
      Lval : List(R)
      Lpval : List(%)
      Lvar : List(VarSet)

    --define
      0  == 0$R::%
      1  == 1$R::%


      zero? p == p case R and zero?(p)$R
      one? p == p case R and one?(p)$R
    -- a local function
      red(p:%):% ==
         p case R => 0
         if ground?(reductum p.ts) then leadingCoefficient(reductum p.ts) else [p.v,reductum p.ts]$VPoly

      numberOfMonomials(p): NonNegativeInteger ==
        p case R => 
          zero?(p)$R => 0
          1
        +/[numberOfMonomials q for q in coefficients(p.ts)]

      coerce(mvar):% == [mvar,monomial(1,1)$D]$VPoly

      monomial? p ==
        p case R => true
        sup : D := p.ts
        1 ~= numberOfMonomials(sup) => false
        monomial? leadingCoefficient(sup)$D

--    local
      moreThanOneVariable?: % -> Boolean

      moreThanOneVariable? p == 
         p case R => false
         q:=p.ts
         any?(not ground? #1 ,coefficients q) => true
         false

      -- if we already know we use this (slighlty) faster function
      univariateKnown: % -> SparseUnivariatePolynomial R 

      univariateKnown p == 
        p case R => (leadingCoefficient p) :: SparseUnivariatePolynomial(R)
	monomial( leadingCoefficient p,degree p.ts)+ univariateKnown(red p)

      univariate p ==
        p case R =>(leadingCoefficient p) :: SparseUnivariatePolynomial(R)
        moreThanOneVariable?  p => error "not univariate"
        monomial( leadingCoefficient p,degree p.ts)+ univariate(red p)

      multivariate (u:SparseUnivariatePolynomial(R),var:VarSet) ==
        ground? u => (leadingCoefficient u) ::%
        [var,monomial(leadingCoefficient u,degree u)$D]$VPoly +
           multivariate(reductum u,var)

      univariate(p:%,mvar:VarSet):SparseUnivariatePolynomial(%) ==
        p case R or mvar>p.v  => monomial(p,0)$D
        pt:=p.ts
        mvar=p.v => pt
        monomial(1,p.v,degree pt)*univariate(leadingCoefficient pt,mvar)+
          univariate(red p,mvar)

--  a local functions, used in next definition
      unlikeUnivReconstruct(u:SparseUnivariatePolynomial(%),mvar:VarSet):% ==
        zero? (d:=degree u) => coefficient(u,0)
        monomial(leadingCoefficient u,mvar,d)+
            unlikeUnivReconstruct(reductum u,mvar)

      multivariate(u:SparseUnivariatePolynomial(%),mvar:VarSet):% ==
        ground? u => coefficient(u,0)
        uu:=u
        while not zero? uu repeat
          cc:=leadingCoefficient uu
          cc case R or mvar > cc.v => uu:=reductum uu
          return unlikeUnivReconstruct(u,mvar)
        [mvar,u]$VPoly

      ground?(p:%):Boolean ==
        p case R => true
        false

--      const p ==
--        p case R => p
--        error "the polynomial is not a constant"

      monomial(p,mvar,k1) ==
        zero? k1 or zero? p => p
        p case R or mvar>p.v => [mvar,monomial(p,k1)$D]$VPoly
        p*[mvar,monomial(1,k1)$D]$VPoly

      monomial(c:R,e:IndexedExponents(VarSet)):% ==
        zero? e => (c::%)
        monomial(1,leadingSupport e, leadingCoefficient e) *
            monomial(c,reductum e)

      coefficient(p:%, e:IndexedExponents(VarSet)) : R ==
        zero? e =>
          p case R  => p::R
          coefficient(coefficient(p.ts,0),e)
        p case R => 0
        ve := leadingSupport e
        vp := p.v
        ve < vp =>
          coefficient(coefficient(p.ts,0),e)
        ve > vp => 0
        coefficient(coefficient(p.ts,leadingCoefficient e),reductum e)

--    coerce(e:IndexedExponents(VarSet)) : % ==
--      e = 0 => 1
--      monomial(1,leadingSupport e, leadingCoefficient e) *
--          (reductum e)::%

--    retract(p:%):IndexedExponents(VarSet) ==
--      q:Union(IndexedExponents(VarSet),"failed"):=retractIfCan p
--      q :: IndexedExponents(VarSet)

--    retractIfCan(p:%):Union(IndexedExponents(VarSet),"failed") ==
--      p = 0 => degree p
--      reductum(p)=0 and leadingCoefficient(p)=1 => degree p
--      "failed"

      coerce(n) == n::R::%
      coerce(c) == c::%
      characteristic == characteristic$R

      recip(p) ==
        p case R => (uu:=recip(p::R);uu case "failed" => "failed"; uu::%)
        "failed"

      - p ==
          p case R => -$R p
          [p.v, - p.ts]$VPoly
      n * p  ==
          p case R => n * p::R
          mvar:=p.v
          up:=n*p.ts
          if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly
      c * p  ==
          c = 1 => p
          p case R => c * p::R
          mvar:=p.v
          up:=c*p.ts
          if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly
      p1 + p2  ==
         p1 case R and p2 case R => p1 +$R p2
         p1 case R => [p2.v, p1::D + p2.ts]$VPoly
	 p2 case R => [p1.v,  p1.ts + p2::D]$VPoly
	 p1.v = p2.v => 
              mvar:=p1.v
              up:=p1.ts+p2.ts
              if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly
	 p1.v < p2.v =>
              [p2.v, p1::D + p2.ts]$VPoly
         [p1.v, p1.ts + p2::D]$VPoly

      p1 - p2  ==
         p1 case R and p2 case R => p1 -$R p2
         p1 case R => [p2.v, p1::D - p2.ts]$VPoly
         p2 case R => [p1.v,  p1.ts - p2::D]$VPoly
         p1.v = p2.v =>
              mvar:=p1.v
              up:=p1.ts-p2.ts
              if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly
         p1.v < p2.v =>
              [p2.v, p1::D - p2.ts]$VPoly
         [p1.v, p1.ts - p2::D]$VPoly

      p1 = p2  ==
         p1 case R =>
             p2 case R => p1 =$R p2
             false
         p2 case R => false
         p1.v = p2.v => p1.ts = p2.ts
         false

      p1 * p2  ==
         p1 case R => p1::R * p2
         p2 case R => 
            mvar:=p1.v
            up:=p1.ts*p2
            if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly
	 p1.v = p2.v => 
            mvar:=p1.v
            up:=p1.ts*p2.ts
            if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly
         p1.v > p2.v => 
            mvar:=p1.v
            up:=p1.ts*p2
            if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly
            --- p1.v < p2.v 
         mvar:=p2.v
         up:=p1*p2.ts
         if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly

      p ** kp == p ** (kp pretend NonNegativeInteger )
      p ** k  ==
         p case R => p::R ** k
         -- univariate special case 
         not moreThanOneVariable? p => 
             multivariate( (univariateKnown p) ** k , p.v)
         mvar:=p.v
         up:=p.ts ** k
         if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly

      if R has IntegralDomain then
         UnitCorrAssoc ==> Record(unit:%,canonical:%,associate:%)
         unitNormal(p) ==
            u,c,a:R
            p case R =>
              (u,c,a):= unitNormal(p::R)$R
              [u::%,c::%,a::%]$UnitCorrAssoc
            (u,c,a):= unitNormal(leadingCoefficient(p))$R
            [u::%,(a*p)::%,a::%]$UnitCorrAssoc
         unitCanonical(p) ==
            p case R => unitCanonical(p::R)$R
            (u,c,a):= unitNormal(leadingCoefficient(p))$R
            a*p
         unit? p ==
            p case R => unit?(p::R)$R
            false
         associates?(p1,p2) ==
            p1 case R => p2 case R and associates?(p1,p2)$R
            p2 case VPoly and p1.v = p2.v and associates?(p1.ts,p2.ts)

         if R has approximate then
           p1  exquo  p2  ==
              p1 case R and p2 case R =>
                a:= (p1::R  exquo  p2::R)
                if a case "failed" then "failed" else a::%
              zero? p1 => p1
              one? p2 => p1
              p1 case R or p2 case VPoly and p1.v < p2.v => "failed"
              p2 case R or p1.v > p2.v =>
                 a:= (p1.ts  exquo  p2::D)
                 a case "failed" => "failed"
                 [p1.v,a]$VPoly::%
              -- The next test is useful in the case that R has inexact
              -- arithmetic (in particular when it is Interval(...)).
              -- In the case where the test succeeds, empirical evidence
              -- suggests that it can speed up the computation several times,
              -- but in other cases where there are a lot of variables
              -- and p1 and p2 differ only in the low order terms (e.g. p1=p2+1)
              -- it slows exquo down by about 15-20%.
              p1 = p2 => 1
              a:= p1.ts  exquo  p2.ts
              a case "failed" => "failed"
              mvar:=p1.v
              up:SUP %:=a
              if ground? (up) then leadingCoefficient(up) else [mvar,up]$VPoly::%
         else
           p1  exquo  p2  ==
              p1 case R and p2 case R =>
                a:= (p1::R  exquo  p2::R)
                if a case "failed" then "failed" else a::%
              zero? p1 => p1
              one? p2 => p1
              p1 case R or p2 case VPoly and p1.v < p2.v => "failed"
              p2 case R or p1.v > p2.v =>
                 a:= (p1.ts  exquo  p2::D)
                 a case "failed" => "failed"
                 [p1.v,a]$VPoly::%
              a:= p1.ts  exquo  p2.ts
              a case "failed" => "failed"
              mvar:=p1.v
              up:SUP %:=a
              if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly::%

      map(fn,p) ==
         p case R => fn(p)
         mvar:=p.v
         up:=map(map(fn,#1),p.ts)
         if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly

      if R has Field then
        (p : %) / (r : R) == inv(r) * p

      if R has GcdDomain then
        content(p) ==
           p case R => p
           c :R :=0
           up:=p.ts
           while not(zero? up) and not(one? c) repeat
               c:=gcd(c,content leadingCoefficient(up))
               up := reductum up
           c

      if R has EuclideanDomain and R has CharacteristicZero and not(R has FloatingPointSystem)  then
        content(p,mvar) ==
          p case R => p
          gcd(coefficients univariate(p,mvar))$pgcd

        gcd(p1,p2) ==  gcd(p1,p2)$pgcd

        gcd(lp:List %) ==  gcd(lp)$pgcd

        gcdPolynomial(a:SUP $,b:SUP $):SUP $ == gcd(a,b)$pgcd

      else if R has GcdDomain then
        content(p,mvar) ==
          p case R => p
          content univariate(p,mvar)

        gcd(p1,p2) ==
           p1 case R =>
              p2 case R => gcd(p1,p2)$R::%
              zero? p1 => p2
              gcd(p1, content(p2.ts))
           p2 case R =>
              zero? p2 => p1
              gcd(p2, content(p1.ts))
           p1.v < p2.v => gcd(p1, content(p2.ts))
           p1.v > p2.v => gcd(content(p1.ts), p2)
           mvar:=p1.v
           up:=gcd(p1.ts, p2.ts)
           if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly

        if R has FloatingPointSystem then
           -- eventually need a better notion of gcd's over floats
           -- this essentially computes the gcds of the monomial contents
           gcdPolynomial(a:SUP $,b:SUP $):SUP $ ==
	      ground? (a) =>
                  zero? a => b
		  gcd(leadingCoefficient a, content b)::SUP $
	      ground?(b) =>
                  zero? b => b
		  gcd(leadingCoefficient b, content a)::SUP $
	      conta := content a
	      mona:SUP $ := monomial(conta, minimumDegree a)
              if not one? mona then
		   a := (a exquo mona)::SUP $
	      contb := content b
	      monb:SUP $ := monomial(contb, minimumDegree b)
              if not one? monb then
		   b := (b exquo monb)::SUP $
	      mong:SUP $  := monomial(gcd(conta, contb),
                                      min(degree mona, degree monb))
              degree(a) >= degree b =>
		   not((a exquo b) case "failed") =>
			mong * b
		   mong
	      not((b exquo a) case "failed") => mong * a
	      mong

      coerce(p):OutputForm ==
        p case R => (p::R)::OutputForm
        outputForm(p.ts,p.v::OutputForm)

      coefficients p ==
        p case R => list(p :: R)$List(R)
        "append"/[coefficients(p1)$% for p1 in coefficients(p.ts)]

      retract(p:%):R ==
        p case R => p :: R
        error "cannot retract nonconstant polynomial"

      retractIfCan(p:%):Union(R, "failed") ==
        p case R => p::R
        "failed"

--      leadingCoefficientRecursive(p:%):% ==
--         p case R => p
--         leadingCoefficient p.ts

      mymerge:(List VarSet,List VarSet) ->List VarSet
      mymerge(l:List VarSet,m:List VarSet):List VarSet  ==
         empty? l => m
         empty? m => l
         first l = first m => 
            empty? rest l => 
                 setrest!(l,rest m)
                 l
            empty? rest m => l 
	    setrest!(l, mymerge(rest l, rest m))
            l
         first l > first m =>
            empty? rest l => 
                setrest!(l,m) 
                l
            setrest!(l, mymerge(rest l, m))
            l
         empty? rest m => 
             setrest!(m,l)
             m
         setrest!(m,mymerge(l,rest m))
         m
         
      variables p ==
         p case R => empty()
         lv:List VarSet:=empty()
         q := p.ts
         while not zero? q repeat
           lv:=mymerge(lv,variables leadingCoefficient q)
           q := reductum q
         cons(p.v,lv)

      mainVariable p ==
         p case R => "failed"
         p.v

      eval(p,mvar,pval) == univariate(p,mvar)(pval)
      eval(p,mvar,val) ==  univariate(p,mvar)(val)

      evalSortedVarlist(p,Lvar,Lpval):% ==
        p case R => p
        empty? Lvar or empty? Lpval => p
        mvar := Lvar.first
        mvar > p.v => evalSortedVarlist(p,Lvar.rest,Lpval.rest)
        pval := Lpval.first
        pts := map(evalSortedVarlist(#1,Lvar,Lpval),p.ts)
        mvar=p.v =>
             pval case R => pts (pval::R)
             pts pval
        multivariate(pts,p.v)

      eval(p,Lvar,Lpval) ==
	empty? rest Lvar => evalSortedVarlist(p,Lvar,Lpval)
	sorted?(#1 > #2, Lvar) => evalSortedVarlist(p,Lvar,Lpval)
        nlvar := sort(#1 > #2,Lvar)
        nlpval :=
           Lvar = nlvar => Lpval
           [Lpval.position(mvar,Lvar) for mvar in nlvar]
        evalSortedVarlist(p,nlvar,nlpval)

      eval(p,Lvar,Lval) ==
        eval(p,Lvar,[val::% for val in Lval]$(List %)) -- kill?

      degree(p,mvar) ==
        p case R => 0
        mvar= p.v => degree p.ts
        mvar > p.v => 0    -- might as well take advantage of the order
        max(degree(leadingCoefficient p.ts,mvar),degree(red p,mvar))

      degree(p,Lvar)  == [degree(p,mvar)  for mvar in Lvar]

      degree p ==
        p case R => 0
        degree(leadingCoefficient(p.ts)) + monomial(degree(p.ts), p.v)

      minimumDegree p ==
        p case R => 0
        md := minimumDegree p.ts
        minimumDegree(coefficient(p.ts,md)) + monomial(md, p.v)

      minimumDegree(p,mvar) ==
        p case R => 0
        mvar = p.v => minimumDegree p.ts
        md:=minimumDegree(leadingCoefficient p.ts,mvar)
        zero? (p1:=red p) => md
        min(md,minimumDegree(p1,mvar))

      minimumDegree(p,Lvar) ==
        [minimumDegree(p,mvar) for mvar in Lvar]

      totalDegree(p, Lvar) ==
        ground? p => 0
        null setIntersection(Lvar, variables p) => 0
        u := univariate(p, mv := mainVariable(p)::VarSet)
        weight:NonNegativeInteger := (member?(mv,Lvar) => 1; 0)
        tdeg:NonNegativeInteger := 0
        while u ~= 0 repeat
            termdeg:NonNegativeInteger := weight*degree u +
                           totalDegree(leadingCoefficient u, Lvar)
            tdeg := max(tdeg, termdeg)
            u := reductum u
        tdeg

      if R has CommutativeRing then
        differentiate(p,mvar) ==
          p case R => 0
          mvar=p.v =>  
             up:=differentiate p.ts
             if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly
          up:=map(differentiate(#1,mvar),p.ts)
          if ground? up then leadingCoefficient(up) else [p.v,up]$VPoly

      leadingCoefficient(p) ==
         p case R => p
         leadingCoefficient(leadingCoefficient(p.ts))

--      trailingCoef(p) ==
--        p case R => p
--        coef(p.ts,0) case R => coef(p.ts,0)
--        trailingCoef(red p)
--      TrailingCoef(p) == trailingCoef(p)

      leadingMonomial p ==
          p case R => p
          monomial(leadingMonomial leadingCoefficient(p.ts),
                   p.v, degree(p.ts))

      reductum(p) == 
          p case R => 0
          p - leadingMonomial p


--        if R is Integer then
--           pgcd := PolynomialGcdPackage(%,VarSet)
--           gcd(p1,p2) ==
--               gcd(p1,p2)$pgcd
--
--        else if R is RationalNumber then
--           gcd(p1,p2) ==
--               mrat:= MRationalFactorize(VarSet,%)
--               gcd(p1,p2)$mrat
--
--        else gcd(p1,p2) ==
--           p1 case R =>
--              p2 case R => gcd(p1,p2)$R::%
--              p1 = 0 => p2
--              gcd(p1, content(p2.ts))
--           p2 case R =>
--              p2 = 0 => p1
--              gcd(p2, content(p1.ts))
--           p1.v < p2.v => gcd(p1, content(p2.ts))
--           p1.v > p2.v => gcd(content(p1.ts), p2)
--           PSimp(p1.v, gcd(p1.ts, p2.ts))

)abbrev domain POLY Polynomial
++ Author: Dave Barton, Barry Trager
++ Date Created:
++ Date Last Updated:
++ Basic Functions: Ring, degree, eval, coefficient, monomial, differentiate,
++ resultant, gcd
++ Related Constructors: SparseMultivariatePolynomial, MultivariatePolynomial
++ Also See:
++ AMS Classifications:
++ Keywords: polynomial, multivariate
++ References:
++ Description:
++   This type is the basic representation of sparse recursive multivariate
++ polynomials whose variables are arbitrary symbols. The ordering
++ is alphabetic determined by the Symbol type.
++ The coefficient ring may be non commutative,
++ but the variables are assumed to commute.

Polynomial(R:Ring):
  PolynomialCategory(R, IndexedExponents Symbol, Symbol) with
   if R has Algebra Fraction Integer then
     integrate: (%, Symbol) -> %
       ++ integrate(p,x) computes the integral of \spad{p*dx}, i.e.
       ++ integrates the polynomial p with respect to the variable x.
 == SparseMultivariatePolynomial(R, Symbol) add

    import UserDefinedPartialOrdering(Symbol)

    coerce(p:%):OutputForm ==
      (r:= retractIfCan(p)@Union(R,"failed")) case R => r::R::OutputForm
      a :=
        userOrdered?() => largest variables p
        mainVariable(p)::Symbol
      outputForm(univariate(p, a), a::OutputForm)

    if R has Algebra Fraction Integer then
      integrate(p, x) == (integrate univariate(p, x)) (x::%)

)abbrev package POLY2 PolynomialFunctions2
++ Author:
++ Date Created:
++ Date Last Updated:
++ Basic Functions:
++ Related Constructors:
++ Also See:
++ AMS Classifications:
++ Keywords:
++ References:
++ Description:
++   This package takes a mapping between coefficient rings, and lifts
++ it to a mapping between polynomials over those rings.

PolynomialFunctions2(R:Ring, S:Ring): with
  map: (R -> S, Polynomial R) -> Polynomial S
    ++ map(f, p) produces a new polynomial as a result of applying
    ++ the function f to every coefficient of the polynomial p.
 == add
  map(f, p) == map(#1::Polynomial(S), f(#1)::Polynomial(S),
                   p)$PolynomialCategoryLifting(IndexedExponents Symbol,
                                   Symbol, R, Polynomial R, Polynomial S)


)abbrev domain MPOLY MultivariatePolynomial
++ Author: Dave Barton, Barry Trager
++ Date Created:
++ Date Last Updated:
++ Basic Functions: Ring, degree, eval, coefficient, monomial, differentiate,
++ resultant, gcd
++ Related Constructors: SparseMultivariatePolynomial, Polynomial
++ Also See:
++ AMS Classifications:
++ Keywords: polynomial, multivariate
++ References:
++ Description:
++   This type is the basic representation of sparse recursive multivariate
++ polynomials whose variables are from a user specified list of symbols.
++ The ordering is specified by the position of the variable in the list.
++ The coefficient ring may be non commutative,
++ but the variables are assumed to commute.

MultivariatePolynomial(vl:List Symbol, R:Ring)
   ==  SparseMultivariatePolynomial(--SparseUnivariatePolynomial,
           R, OrderedVariableList vl)

