\documentclass{article}
\usepackage{open-axiom}
\begin{document}
\title{\$SPAD/src/algebra newpoly.spad}
\author{Marc Moreno Maza}
\maketitle
\begin{abstract}
\end{abstract}
\eject
\tableofcontents
\eject
Based on the {\bf PseudoRemainderSequence} package, the domain
constructor {\bf NewSparseUnivariatePolynomial} extends the
constructur {\bf SparseUnivariatePolynomial}. 
\section{domain NSUP NewSparseUnivariatePolynomial}
<<domain NSUP NewSparseUnivariatePolynomial>>=
)abbrev domain NSUP NewSparseUnivariatePolynomial
++ Author: Marc Moreno Maza
++ Date Created: 23/07/98
++ Date Last Updated: 14/12/98
++ Basic Operations: 
++ Related Domains:
++ Also See:
++ AMS Classifications:
++ Keywords: 
++ Examples:
++ References: 
++ Description: A post-facto extension for \axiomType{SUP} in order
++ to speed up operations related to pseudo-division and gcd for
++ both \axiomType{SUP} and, consequently, \axiomType{NSMP}.

NewSparseUnivariatePolynomial(R): Exports == Implementation where

  R:Ring
  NNI ==> NonNegativeInteger
  SUPR ==> SparseUnivariatePolynomial R

  Exports == Join(UnivariatePolynomialCategory(R), 
   CoercibleTo(SUPR),RetractableTo(SUPR)) with
     fmecg : (%,NNI,R,%) -> %
        ++ \axiom{fmecg(p1,e,r,p2)} returns \axiom{p1 - r * X**e * p2}
        ++ where \axiom{X} is \axiom{monomial(1,1)}
     monicModulo : ($, $) -> $ 
        ++ \axiom{monicModulo(a,b)} returns \axiom{r} such that \axiom{r} is 
        ++ reduced w.r.t. \axiom{b} and \axiom{b} divides \axiom{a -r}
        ++ where \axiom{b} is monic.
     lazyResidueClass: ($,$) -> Record(polnum:$, polden:R, power:NNI)
        ++ \axiom{lazyResidueClass(a,b)} returns \axiom{[r,c,n]} such that 
        ++ \axiom{r} is reduced w.r.t. \axiom{b} and \axiom{b} divides 
        ++ \axiom{c^n * a - r} where \axiom{c} is \axiom{leadingCoefficient(b)} 
        ++ and \axiom{n} is as small as possible with the previous properties.
     lazyPseudoRemainder: ($,$) -> $
        ++ \axiom{lazyPseudoRemainder(a,b)} returns \axiom{r} if \axiom{lazyResidueClass(a,b)}
        ++ returns \axiom{[r,c,n]}. This lazy pseudo-remainder is computed by 
        ++ means of the \axiomOpFrom{fmecg}{NewSparseUnivariatePolynomial} operation.
     lazyPseudoDivide: ($,$) -> Record(coef:R, gap:NNI, quotient:$, remainder: $)
        ++ \axiom{lazyPseudoDivide(a,b)} returns \axiom{[c,g,q,r]} such that 
        ++ \axiom{c^n * a = q*b +r} and \axiom{lazyResidueClass(a,b)} returns \axiom{[r,c,n]}
        ++ where \axiom{n + g = max(0, degree(b) - degree(a) + 1)}.
     lazyPseudoQuotient: ($,$) -> $
        ++ \axiom{lazyPseudoQuotient(a,b)} returns \axiom{q} if \axiom{lazyPseudoDivide(a,b)}
        ++ returns \axiom{[c,g,q,r]} 
     if R has IntegralDomain
     then 
       subResultantsChain: ($, $) -> List $
         ++ \axiom{subResultantsChain(a,b)} returns the list of the non-zero
         ++ sub-resultants of \axiom{a} and \axiom{b} sorted by increasing 
         ++ degree.
       lastSubResultant: ($, $) -> $
         ++ \axiom{lastSubResultant(a,b)} returns \axiom{resultant(a,b)}
         ++ if \axiom{a} and \axiom{b} has no non-trivial gcd in \axiom{R^(-1) P}
         ++ otherwise the non-zero sub-resultant with smallest index.
       extendedSubResultantGcd: ($, $) -> Record(gcd: $, coef1: $, coef2: $)
         ++ \axiom{extendedSubResultantGcd(a,b)} returns \axiom{[g,ca, cb]} such
         ++ that \axiom{g} is a gcd of \axiom{a} and \axiom{b} in \axiom{R^(-1) P}
         ++ and \axiom{g = ca * a + cb * b}
       halfExtendedSubResultantGcd1: ($, $) -> Record(gcd: $, coef1: $)
         ++ \axiom{halfExtendedSubResultantGcd1(a,b)} returns \axiom{[g,ca]} such that
         ++ \axiom{extendedSubResultantGcd(a,b)} returns \axiom{[g,ca, cb]}
       halfExtendedSubResultantGcd2: ($, $) -> Record(gcd: $, coef2: $)
         ++ \axiom{halfExtendedSubResultantGcd2(a,b)} returns \axiom{[g,cb]} such that
         ++ \axiom{extendedSubResultantGcd(a,b)} returns \axiom{[g,ca, cb]}
       extendedResultant: ($, $) -> Record(resultant: R, coef1: $, coef2: $)
         ++ \axiom{extendedResultant(a,b)} returns  \axiom{[r,ca,cb]} such that 
         ++ \axiom{r} is the resultant of \axiom{a} and \axiom{b} and
         ++ \axiom{r = ca * a + cb * b}
       halfExtendedResultant1: ($, $) -> Record(resultant: R, coef1: $)
         ++ \axiom{halfExtendedResultant1(a,b)} returns \axiom{[r,ca]} such that 
         ++ \axiom{extendedResultant(a,b)}  returns \axiom{[r,ca, cb]} 
       halfExtendedResultant2: ($, $) -> Record(resultant: R, coef2: $)
         ++ \axiom{halfExtendedResultant2(a,b)} returns \axiom{[r,ca]} such that 
         ++ \axiom{extendedResultant(a,b)} returns \axiom{[r,ca, cb]} 

  Implementation == SparseUnivariatePolynomial(R) add
   
     Term == Record(k:NonNegativeInteger,c:R)
     Rep ==> List Term

     rep(s:$):Rep == s pretend Rep
     per(l:Rep):$ == l pretend $

     coerce (p:$):SUPR == 
       p pretend SUPR

     coerce (p:SUPR):$ == 
       p pretend $

     retractIfCan (p:$) : Union(SUPR,"failed") == 
       (p pretend SUPR)::Union(SUPR,"failed")

     monicModulo(x,y) ==
		zero? y => 
		   error "in monicModulo$NSUP: division by 0"
		ground? y =>
		   error "in monicModulo$NSUP: ground? #2"
                yy := rep y
                not one? (yy.first.c) => 
		   error "in monicModulo$NSUP: not monic #2"
                xx := rep x; empty? xx => x
                e := yy.first.k; y := per(yy.rest)                
                -- while (not empty? xx) repeat
                repeat
                  if (u:=subtractIfCan(xx.first.k,e)) case nothing then break
                  xx:= rep fmecg(per rest(xx), u, xx.first.c, y)
                  if empty? xx then break
                per xx

     lazyResidueClass(x,y) ==
		zero? y => 
		   error "in lazyResidueClass$NSUP: division by 0"
		ground? y =>
		   error "in lazyResidueClass$NSUP: ground? #2"
                yy := rep y; co := yy.first.c; xx: Rep := rep x
                empty? xx => [x, co, 0]
                pow: NNI := 0; e := yy.first.k; y := per(yy.rest); 
                repeat
                  if (u:=subtractIfCan(xx.first.k,e)) case nothing then break
                  xx:= rep fmecg(co * per rest(xx), u, xx.first.c, y)
                  pow := pow + 1
                  if empty? xx then break
                [per xx, co, pow]

     lazyPseudoRemainder(x,y) ==
		zero? y => 
		   error "in lazyPseudoRemainder$NSUP: division by 0"
		ground? y =>
		   error "in lazyPseudoRemainder$NSUP: ground? #2"
		ground? x => x
                yy := rep y; co := yy.first.c
                one? co => monicModulo(x,y)
                (co = -1) => - monicModulo(-x,-y)
                xx:= rep x; e := yy.first.k; y := per(yy.rest)
                repeat
                  if (u:=subtractIfCan(xx.first.k,e)) case nothing then break
                  xx:= rep fmecg(co * per rest(xx), u, xx.first.c, y)
                  if empty? xx then break
                per xx

     lazyPseudoDivide(x,y) ==
		zero? y => 
		   error "in lazyPseudoDivide$NSUP: division by 0"
		ground? y =>
		   error "in lazyPseudoDivide$NSUP: ground? #2"
                yy := rep y; e := yy.first.k; 
                xx: Rep := rep x; co := yy.first.c
		(empty? xx) or (xx.first.k < e) => [co,0,0,x]
                pow: NNI := subtractIfCan(xx.first.k,e)::NNI + 1
                qq: Rep := []; y := per(yy.rest)
                repeat
                  if (u:=subtractIfCan(xx.first.k,e)) case nothing then break
                  qq := cons([u::NNI, xx.first.c]$Term, rep (co * per qq))
                  xx := rep fmecg(co * per rest(xx), u, xx.first.c, y)
                  pow := subtractIfCan(pow,1)::NNI
                  if empty? xx then break
                [co, pow, per reverse qq, per xx]

     lazyPseudoQuotient(x,y) ==
		zero? y => 
		   error "in lazyPseudoQuotient$NSUP: division by 0"
		ground? y =>
		   error "in lazyPseudoQuotient$NSUP: ground? #2"
                yy := rep y; e := yy.first.k; xx: Rep := rep x
		(empty? xx) or (xx.first.k < e) => 0
                qq: Rep := []; co := yy.first.c; y := per(yy.rest)
                repeat
                  if (u:=subtractIfCan(xx.first.k,e)) case nothing then break
                  qq := cons([u::NNI, xx.first.c]$Term, rep (co * per qq))
                  xx := rep fmecg(co * per rest(xx), u, xx.first.c, y)
                  if empty? xx then break
                per reverse qq

     if R has IntegralDomain
     then 
       pack ==> PseudoRemainderSequence(R, %)

       subResultantGcd(p1,p2) == subResultantGcd(p1,p2)$pack

       subResultantsChain(p1,p2) == chainSubResultants(p1,p2)$pack

       lastSubResultant(p1,p2) == lastSubResultant(p1,p2)$pack

       resultant(p1,p2) == resultant(p1,p2)$pack

       extendedResultant(p1,p2) == 
          re: Record(coef1: $, coef2: $, resultant: R) := resultantEuclidean(p1,p2)$pack
          [re.resultant, re.coef1, re.coef2]

       halfExtendedResultant1(p1:$, p2: $): Record(resultant: R, coef1: $) ==
          re: Record(coef1: $, resultant: R) := semiResultantEuclidean1(p1, p2)$pack
          [re.resultant, re.coef1]

       halfExtendedResultant2(p1:$, p2: $): Record(resultant: R, coef2: $) ==
          re: Record(coef2: $, resultant: R) := semiResultantEuclidean2(p1, p2)$pack
          [re.resultant, re.coef2]

       extendedSubResultantGcd(p1,p2) == 
          re: Record(coef1: $, coef2: $, gcd: $) := subResultantGcdEuclidean(p1,p2)$pack
          [re.gcd, re.coef1, re.coef2]

       halfExtendedSubResultantGcd1(p1:$, p2: $): Record(gcd: $, coef1: $) ==
          re: Record(coef1: $, gcd: $) := semiSubResultantGcdEuclidean1(p1, p2)$pack
          [re.gcd, re.coef1]

       halfExtendedSubResultantGcd2(p1:$, p2: $): Record(gcd: $, coef2: $) ==
          re: Record(coef2: $, gcd: $) := semiSubResultantGcdEuclidean2(p1, p2)$pack
          [re.gcd, re.coef2]

       pseudoDivide(x,y) ==
		zero? y => 
		   error "in pseudoDivide$NSUP: division by 0"
		ground? y =>
		   error "in pseudoDivide$NSUP: ground? #2"
                yy := rep y; e := yy.first.k
                xx: Rep := rep x; co := yy.first.c
		(empty? xx) or (xx.first.k < e) => [co,0,x]
                pow: NNI := subtractIfCan(xx.first.k,e)::NNI + 1
                qq: Rep := []; y := per(yy.rest)
                repeat
                  if (u:=subtractIfCan(xx.first.k,e)) case nothing then break
                  qq := cons([u::NNI, xx.first.c]$Term, rep (co * per qq))
                  xx := rep fmecg(co * per rest(xx), u, xx.first.c, y)
                  pow := subtractIfCan(pow,1)::NNI
                  if empty? xx then break
                zero? pow => [co, per reverse qq, per xx]
                default: R := co ** pow
                q := default * (per reverse qq)
                x := default * (per xx)
                [co, q, x]

       pseudoQuotient(x,y) ==
		zero? y => 
		   error "in pseudoDivide$NSUP: division by 0"
		ground? y =>
		   error "in pseudoDivide$NSUP: ground? #2"
                yy := rep y; e := yy.first.k; xx: Rep := rep x
		(empty? xx) or (xx.first.k < e) => 0
                pow: NNI := subtractIfCan(xx.first.k,e)::NNI + 1
                qq: Rep := []; co := yy.first.c; y := per(yy.rest)
                repeat
                  if (u:=subtractIfCan(xx.first.k,e)) case nothing then break
                  qq := cons([u::NNI, xx.first.c]$Term, rep (co * per qq))
                  xx := rep fmecg(co * per rest(xx), u, xx.first.c, y)
                  pow := subtractIfCan(pow,1)::NNI
                  if empty? xx then break
                zero? pow => per reverse qq
                (co ** pow) * (per reverse qq)

@
\section{package NSUP2 NewSparseUnivariatePolynomialFunctions2}
<<package NSUP2 NewSparseUnivariatePolynomialFunctions2>>=
)abbrev package NSUP2 NewSparseUnivariatePolynomialFunctions2
++ 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.

NewSparseUnivariatePolynomialFunctions2(R:Ring, S:Ring): with
  map:(R->S,NewSparseUnivariatePolynomial R) -> NewSparseUnivariatePolynomial S
    ++ \axiom{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,
           NewSparseUnivariatePolynomial R, S, NewSparseUnivariatePolynomial S)

@
\section{category RPOLCAT RecursivePolynomialCategory}
<<category RPOLCAT RecursivePolynomialCategory>>=
)abbrev category RPOLCAT RecursivePolynomialCategory
++ Author: Marc Moreno Maza
++ Date Created: 04/22/1994
++ Date Last Updated: 14/12/1998
++ Basic Functions: mvar, mdeg, init, head, tail, prem, lazyPrem
++ Related Constructors:
++ Also See: 
++ AMS Classifications:
++ Keywords: polynomial, multivariate, ordered variables set
++ References:
++ Description:
++ A category for general multi-variate polynomials with coefficients 
++ in a ring, variables in an ordered set, and exponents from an 
++ ordered abelian monoid, with a \axiomOp{sup} operation.
++ When not constant, such a polynomial is viewed as a univariate polynomial in its 
++ main variable w. r. t. to the total ordering on the elements in the ordered set, so that some
++ operations usually defined for univariate polynomials make sense here.

RecursivePolynomialCategory(R:Ring, E:OrderedAbelianMonoidSup, V:OrderedSet): Category == 
  PolynomialCategory(R, E, V) with
     mvar : $ -> V
         ++ \axiom{mvar(p)} returns an error if \axiom{p} belongs to \axiom{R}, 
         ++ otherwise returns its main variable w. r. t. to the total ordering 
         ++ on the elements in \axiom{V}.
     mdeg  : $ -> NonNegativeInteger 
         ++ \axiom{mdeg(p)} returns an error if \axiom{p} is \axiom{0}, 
         ++ otherwise, if \axiom{p} belongs to \axiom{R} returns \axiom{0}, 
         ++ otherwise, returns the degree of \axiom{p} in its main variable.
     init : $ -> $
         ++ \axiom{init(p)} returns an error if \axiom{p} belongs to \axiom{R}, 
         ++ otherwise returns its leading coefficient, where \axiom{p} is viewed 
         ++ as a univariate polynomial in its main variable.
     head  : $ -> $
         ++ \axiom{head(p)} returns \axiom{p} if \axiom{p} belongs to \axiom{R}, 
         ++ otherwise returns its leading term (monomial in the AXIOM sense), 
         ++ where \axiom{p} is viewed as a univariate polynomial in its main variable.
     tail  : $ -> $
         ++ \axiom{tail(p)} returns its reductum, where \axiom{p} is viewed as a univariate
         ++ polynomial in its main variable.
     deepestTail : $ -> $
         ++ \axiom{deepestTail(p)} returns \axiom{0} if \axiom{p} belongs to \axiom{R}, 
         ++ otherwise returns tail(p), if \axiom{tail(p)} belongs to  \axiom{R}
         ++ or \axiom{mvar(tail(p)) < mvar(p)}, otherwise returns \axiom{deepestTail(tail(p))}.
     iteratedInitials : $ -> List $ 
         ++ \axiom{iteratedInitials(p)} returns \axiom{[]} if \axiom{p} belongs to \axiom{R}, 
         ++ otherwise returns the list of the iterated initials of \axiom{p}.
     deepestInitial : $ -> $ 
         ++ \axiom{deepestInitial(p)} returns an error if \axiom{p} belongs to \axiom{R}, 
         ++ otherwise returns the last term of \axiom{iteratedInitials(p)}.
     leadingCoefficient : ($,V) -> $
         ++ \axiom{leadingCoefficient(p,v)} returns the leading coefficient of \axiom{p}, 
         ++ where \axiom{p} is viewed as A univariate polynomial in \axiom{v}.
     reductum  : ($,V) -> $
         ++ \axiom{reductum(p,v)} returns the reductum of \axiom{p}, where \axiom{p} is viewed as 
         ++ a univariate polynomial in \axiom{v}. 
     monic? : $ -> Boolean
         ++ \axiom{monic?(p)} returns false if \axiom{p} belongs to \axiom{R}, 
         ++ otherwise returns true iff \axiom{p} is monic as a univariate polynomial 
         ++ in its main variable.
     quasiMonic? : $ -> Boolean
         ++ \axiom{quasiMonic?(p)} returns false if \axiom{p} belongs to \axiom{R}, 
         ++ otherwise returns true iff the initial of \axiom{p} lies in the base ring \axiom{R}.
     mainMonomial : $ -> $ 
         ++ \axiom{mainMonomial(p)} returns an error if \axiom{p} is \axiom{O}, 
         ++ otherwise, if \axiom{p} belongs to \axiom{R} returns \axiom{1},
         ++ otherwise, \axiom{mvar(p)} raised to the power \axiom{mdeg(p)}.
     leastMonomial : $ -> $ 
         ++ \axiom{leastMonomial(p)} returns an error if \axiom{p} is \axiom{O}, 
         ++ otherwise, if \axiom{p} belongs to \axiom{R} returns \axiom{1},
         ++ otherwise, the monomial of \axiom{p} with lowest degree, 
         ++ where \axiom{p} is viewed as a univariate polynomial in its main variable.
     mainCoefficients : $ -> List $ 
         ++ \axiom{mainCoefficients(p)} returns an error if \axiom{p} is \axiom{O}, 
         ++ otherwise, if \axiom{p} belongs to \axiom{R} returns [p], 
         ++ otherwise returns the list of the coefficients of \axiom{p}, 
         ++ where \axiom{p} is viewed as a univariate polynomial in its main variable.
     mainMonomials : $ -> List $ 
         ++ \axiom{mainMonomials(p)} returns an error if \axiom{p} is \axiom{O}, 
         ++ otherwise, if \axiom{p} belongs to \axiom{R} returns [1], 
         ++ otherwise returns the list of the monomials of \axiom{p}, 
         ++ where \axiom{p} is viewed as a univariate polynomial in its main variable.
     RittWuCompare : ($, $) -> Union(Boolean,"failed")
         ++ \axiom{RittWuCompare(a,b)} returns \axiom{"failed"} if \axiom{a} and \axiom{b} have same rank w.r.t. 
         ++ Ritt and Wu Wen Tsun ordering using the refinement of Lazard, 
         ++ otherwise returns \axiom{infRittWu?(a,b)}.
     infRittWu?  : ($, $) -> Boolean
         ++ \axiom{infRittWu?(a,b)} returns true if \axiom{a} is less than \axiom{b}
         ++ w.r.t. the Ritt and Wu Wen Tsun ordering using the refinement of Lazard.
     supRittWu? : ($, $) -> Boolean
         ++ \axiom{supRittWu?(a,b)} returns true if \axiom{a} is greater than \axiom{b}
         ++ w.r.t. the Ritt and Wu Wen Tsun ordering using the refinement of Lazard.
     reduced? : ($,$) -> Boolean
         ++ \axiom{reduced?(a,b)} returns true iff \axiom{degree(a,mvar(b)) < mdeg(b)}.
     reduced? : ($,List($)) -> Boolean
         ++ \axiom{reduced?(q,lp)} returns true iff \axiom{reduced?(q,p)} holds 
         ++ for every \axiom{p} in \axiom{lp}.
     headReduced? : ($,$) -> Boolean
         ++ \axiom{headReduced?(a,b)} returns true iff \axiom{degree(head(a),mvar(b)) < mdeg(b)}.
     headReduced? : ($,List($)) -> Boolean
         ++ \axiom{headReduced?(q,lp)} returns true iff \axiom{headReduced?(q,p)} holds 
         ++ for every \axiom{p} in \axiom{lp}.
     initiallyReduced? : ($,$) -> Boolean
         ++ \axiom{initiallyReduced?(a,b)} returns false iff there exists an iterated initial 
         ++ of \axiom{a} which is not reduced w.r.t \axiom{b}.
     initiallyReduced? :  ($,List($)) -> Boolean
         ++ \axiom{initiallyReduced?(q,lp)} returns true iff \axiom{initiallyReduced?(q,p)} holds
         ++ for every \axiom{p} in \axiom{lp}.
     normalized? : ($,$) -> Boolean
         ++ \axiom{normalized?(a,b)} returns true iff \axiom{a} and its iterated initials have
         ++ degree zero w.r.t. the main variable of \axiom{b}
     normalized? : ($,List($)) -> Boolean
         ++ \axiom{normalized?(q,lp)} returns true iff \axiom{normalized?(q,p)} holds 
         ++ for every \axiom{p} in \axiom{lp}.
     prem : ($, $) -> $
         ++ \axiom{prem(a,b)} computes the pseudo-remainder of \axiom{a} by \axiom{b}, 
         ++ both viewed as univariate polynomials in the main variable of \axiom{b}.
     pquo : ($, $) -> $
         ++ \axiom{pquo(a,b)} computes the pseudo-quotient of \axiom{a} by \axiom{b}, 
         ++ both viewed as univariate polynomials in the main variable of \axiom{b}.
     prem : ($, $, V) -> $
         ++ \axiom{prem(a,b,v)} computes the pseudo-remainder of \axiom{a} by \axiom{b}, 
         ++ both viewed as univariate polynomials in \axiom{v}.
     pquo : ($, $, V) -> $
         ++ \axiom{pquo(a,b,v)} computes the pseudo-quotient of \axiom{a} by \axiom{b}, 
         ++ both viewed as univariate polynomials in \axiom{v}.
     lazyPrem : ($, $) ->  $
         ++ \axiom{lazyPrem(a,b)} returns the polynomial \axiom{r} reduced w.r.t. \axiom{b} 
         ++ and such that \axiom{b} divides \axiom{init(b)^e a - r} where \axiom{e} 
         ++ is the number of steps of this pseudo-division.
     lazyPquo : ($, $) ->  $
         ++ \axiom{lazyPquo(a,b)} returns the polynomial \axiom{q} such that 
         ++ \axiom{lazyPseudoDivide(a,b)} returns \axiom{[c,g,q,r]}.
     lazyPrem : ($, $, V) -> $
         ++ \axiom{lazyPrem(a,b,v)} returns the polynomial \axiom{r} 
         ++ reduced w.r.t. \axiom{b} viewed as univariate polynomials in the variable 
         ++ \axiom{v} such that \axiom{b} divides \axiom{init(b)^e a - r}
         ++ where \axiom{e} is the number of steps of this pseudo-division.
     lazyPquo : ($, $, V) ->  $
         ++ \axiom{lazyPquo(a,b,v)} returns the polynomial \axiom{q} such that 
         ++ \axiom{lazyPseudoDivide(a,b,v)} returns \axiom{[c,g,q,r]}.
     lazyPremWithDefault : ($, $) ->  Record (coef : $, gap : NonNegativeInteger, remainder : $)
         ++ \axiom{lazyPremWithDefault(a,b)} returns \axiom{[c,g,r]}
         ++ such that \axiom{r = lazyPrem(a,b)} and \axiom{(c**g)*r = prem(a,b)}.
     lazyPremWithDefault : ($, $, V) ->  Record (coef : $, gap : NonNegativeInteger, remainder : $)
         ++ \axiom{lazyPremWithDefault(a,b,v)} returns \axiom{[c,g,r]} 
         ++ such that \axiom{r = lazyPrem(a,b,v)} and \axiom{(c**g)*r = prem(a,b,v)}.
     lazyPseudoDivide : ($,$) ->  Record(coef:$, gap: NonNegativeInteger,quotient:$, remainder:$)
         ++ \axiom{lazyPseudoDivide(a,b)} returns \axiom{[c,g,q,r]} 
         ++ such that \axiom{[c,g,r] = lazyPremWithDefault(a,b)} and
         ++ \axiom{q} is the pseudo-quotient computed in this lazy pseudo-division.
     lazyPseudoDivide : ($,$,V) ->  Record(coef:$, gap:NonNegativeInteger, quotient:$, remainder: $)
         ++ \axiom{lazyPseudoDivide(a,b,v)} returns \axiom{[c,g,q,r]} such that  
         ++ \axiom{r = lazyPrem(a,b,v)}, \axiom{(c**g)*r = prem(a,b,v)} and \axiom{q}
         ++ is the pseudo-quotient computed in this lazy pseudo-division.
     pseudoDivide : ($, $) -> Record (quotient : $, remainder : $)
         ++ \axiom{pseudoDivide(a,b)} computes \axiom{[pquo(a,b),prem(a,b)]}, both
         ++ polynomials viewed as univariate polynomials in the main variable of \axiom{b}, 
         ++ if \axiom{b} is not a constant polynomial.
     monicModulo : ($, $) -> $ 
         ++ \axiom{monicModulo(a,b)} computes \axiom{a mod b}, if \axiom{b} is 
         ++ monic as univariate polynomial in its main variable.
     lazyResidueClass : ($,$) -> Record(polnum:$, polden:$, power:NonNegativeInteger)
         ++ \axiom{lazyResidueClass(a,b)} returns \axiom{[p,q,n]} where \axiom{p / q**n} 
         ++ represents the residue class of \axiom{a} modulo \axiom{b} 
         ++ and \axiom{p} is reduced w.r.t. \axiom{b} and \axiom{q} is \axiom{init(b)}.
     headReduce: ($, $) ->  $
         ++ \axiom{headReduce(a,b)} returns a polynomial \axiom{r} such that 
         ++ \axiom{headReduced?(r,b)} holds and there exists an integer \axiom{e}
         ++ such that \axiom{init(b)^e a - r} is zero modulo \axiom{b}.
     initiallyReduce: ($, $) ->  $
         ++ \axiom{initiallyReduce(a,b)} returns a polynomial \axiom{r} such that 
         ++ \axiom{initiallyReduced?(r,b)} holds and there exists an integer \axiom{e}
         ++ such that \axiom{init(b)^e a - r} is zero modulo \axiom{b}.

     if (V has ConvertibleTo(Symbol))
     then 
       CoercibleTo(Polynomial R)
       ConvertibleTo(Polynomial R)
       if R has Algebra Fraction Integer
         then 
           retractIfCan : Polynomial Fraction Integer -> Union($,"failed")
               ++ \axiom{retractIfCan(p)} returns \axiom{p} as an element of the current domain
               ++ if all its variables belong to \axiom{V}.
           retract : Polynomial Fraction Integer -> $
               ++ \axiom{retract(p)} returns \axiom{p} as an element of the current domain
               ++ if \axiom{retractIfCan(p)} does not return "failed", otherwise an error 
               ++ is produced.
           convert : Polynomial Fraction Integer -> $
               ++ \axiom{convert(p)} returns the same as \axiom{retract(p)}.
           retractIfCan : Polynomial Integer -> Union($,"failed")
               ++ \axiom{retractIfCan(p)} returns \axiom{p} as an element of the current domain
               ++ if all its variables belong to \axiom{V}.
           retract : Polynomial Integer -> $
               ++ \axiom{retract(p)} returns \axiom{p} as an element of the current domain
               ++ if \axiom{retractIfCan(p)} does not return "failed", otherwise an error 
               ++ is produced.
           convert : Polynomial Integer -> $
               ++ \axiom{convert(p)} returns the same as \axiom{retract(p)}
           if not (R has QuotientFieldCategory(Integer))
             then
               retractIfCan : Polynomial R -> Union($,"failed")
                 ++ \axiom{retractIfCan(p)} returns \axiom{p} as an element of the current domain
                 ++ if all its variables belong to \axiom{V}.
               retract : Polynomial R -> $
                 ++ \axiom{retract(p)} returns \axiom{p} as an element of the current domain
                 ++ if \axiom{retractIfCan(p)} does not return "failed", otherwise an error 
                 ++ is produced.
       if (R has Algebra Integer) and not(R has Algebra Fraction Integer)
         then 
           retractIfCan : Polynomial Integer -> Union($,"failed")
               ++ \axiom{retractIfCan(p)} returns \axiom{p} as an element of the current domain
               ++ if all its variables belong to \axiom{V}.
           retract : Polynomial Integer -> $
               ++ \axiom{retract(p)} returns \axiom{p} as an element of the current domain
               ++ if \axiom{retractIfCan(p)} does not return "failed", otherwise an error 
               ++ is produced.
           convert : Polynomial Integer -> $
               ++ \axiom{convert(p)} returns the same as \axiom{retract(p)}.
           if not (R has IntegerNumberSystem)
             then
               retractIfCan : Polynomial R -> Union($,"failed")
                 ++ \axiom{retractIfCan(p)} returns \axiom{p} as an element of the current domain
                 ++ if all its variables belong to \axiom{V}.
               retract : Polynomial R -> $
                 ++ \axiom{retract(p)} returns \axiom{p} as an element of the current domain
                 ++ if \axiom{retractIfCan(p)} does not return "failed", otherwise an error 
                 ++ is produced.
       if not(R has Algebra Integer) and not(R has Algebra Fraction Integer)
         then 
           retractIfCan : Polynomial R -> Union($,"failed")
             ++ \axiom{retractIfCan(p)} returns \axiom{p} as an element of the current domain
             ++ if all its variables belong to \axiom{V}.
           retract : Polynomial R -> $
             ++ \axiom{retract(p)} returns \axiom{p} as an element of the current domain
             ++ if \axiom{retractIfCan(p)} does not return "failed", otherwise an error 
             ++ is produced.
       convert : Polynomial R -> $
         ++ \axiom{convert(p)} returns \axiom{p} as an element of the current domain if all 
         ++ its variables belong to \axiom{V}, otherwise an error is produced.

       if R has RetractableTo(Integer)
       then
         ConvertibleTo(String)

     if R has IntegralDomain
     then
       primPartElseUnitCanonical : $ -> $
           ++ \axiom{primPartElseUnitCanonical(p)} returns \axiom{primitivePart(p)}
           ++ if \axiom{R} is a gcd-domain, otherwise \axiom{unitCanonical(p)}.
       primPartElseUnitCanonical! : $ -> $
           ++ \axiom{primPartElseUnitCanonical!(p)} replaces  \axiom{p} 
           ++ by \axiom{primPartElseUnitCanonical(p)}.
       exactQuotient : ($,R) -> $
           ++ \axiom{exactQuotient(p,r)} computes the exact quotient of \axiom{p} 
           ++ by \axiom{r}, which is assumed to be a divisor of \axiom{p}. 
           ++ No error is returned if this exact quotient fails!
       exactQuotient! : ($,R) -> $
           ++ \axiom{exactQuotient!(p,r)} replaces \axiom{p} by \axiom{exactQuotient(p,r)}.
       exactQuotient : ($,$) -> $
           ++ \axiom{exactQuotient(a,b)} computes the exact quotient of \axiom{a} 
           ++ by \axiom{b}, which is assumed to be a divisor of \axiom{a}. 
           ++ No error is returned if this exact quotient fails!
       exactQuotient! : ($,$) -> $
           ++ \axiom{exactQuotient!(a,b)} replaces \axiom{a} by \axiom{exactQuotient(a,b)}
       subResultantGcd : ($, $) -> $ 
           ++ \axiom{subResultantGcd(a,b)} computes a gcd of \axiom{a} and \axiom{b} 
           ++ where \axiom{a} and \axiom{b} are assumed to have the same main variable \axiom{v}
           ++ and are viewed as univariate polynomials in \axiom{v} with coefficients
           ++ in the fraction field of the polynomial ring generated by their other variables 
           ++ over \axiom{R}.
       extendedSubResultantGcd : ($, $) -> Record (gcd : $, coef1 : $, coef2 : $)
           ++ \axiom{extendedSubResultantGcd(a,b)} returns \axiom{[ca,cb,r]} 
           ++ such that \axiom{r} is \axiom{subResultantGcd(a,b)} and we have
           ++ \axiom{ca * a + cb * cb = r} .
       halfExtendedSubResultantGcd1: ($, $) -> Record (gcd : $, coef1 : $)
           ++ \axiom{halfExtendedSubResultantGcd1(a,b)} returns \axiom{[g,ca]}
           ++ if \axiom{extendedSubResultantGcd(a,b)} returns \axiom{[g,ca,cb]}
           ++ otherwise produces an error.
       halfExtendedSubResultantGcd2: ($, $) -> Record (gcd : $, coef2 : $)
           ++ \axiom{halfExtendedSubResultantGcd2(a,b)} returns \axiom{[g,cb]}
           ++ if \axiom{extendedSubResultantGcd(a,b)} returns \axiom{[g,ca,cb]}
           ++ otherwise produces an error.
       resultant  : ($, $) -> $ 
           ++ \axiom{resultant(a,b)} computes the resultant of \axiom{a} and \axiom{b} 
           ++ where \axiom{a} and \axiom{b} are assumed to have the same main variable \axiom{v}
           ++ and are viewed as univariate polynomials in \axiom{v}.
       subResultantChain : ($, $) -> List $
           ++ \axiom{subResultantChain(a,b)}, where \axiom{a} and \axiom{b} are not 
           ++ contant polynomials with the same
           ++ main variable, returns the subresultant chain of \axiom{a} and \axiom{b}.
       lastSubResultant: ($, $) -> $
           ++ \axiom{lastSubResultant(a,b)} returns the last non-zero subresultant 
           ++ of \axiom{a} and \axiom{b} where \axiom{a} and \axiom{b} are assumed to have 
           ++ the same main variable \axiom{v} and are viewed as univariate polynomials in \axiom{v}.
       LazardQuotient: ($, $, NonNegativeInteger) -> $
           ++ \axiom{LazardQuotient(a,b,n)} returns \axiom{a**n exquo b**(n-1)}
           ++ assuming that this quotient does not fail.
       LazardQuotient2: ($, $, $, NonNegativeInteger) -> $
           ++ \axiom{LazardQuotient2(p,a,b,n)} returns 
           ++ \axiom{(a**(n-1) * p) exquo b**(n-1)}
           ++ assuming that this quotient does not fail.
       next_subResultant2: ($, $, $, $) -> $
           ++ \axiom{nextsubResultant2(p,q,z,s)} is the multivariate version
           ++ of the operation \axiomOpFrom{next_sousResultant2}{PseudoRemainderSequence} from
           ++ the \axiomType{PseudoRemainderSequence} constructor.

     if R has GcdDomain
     then
       gcd : (R,$) -> R
           ++ \axiom{gcd(r,p)} returns the gcd of \axiom{r} and the content of \axiom{p}.
       primitivePart! : $ -> $
           ++ \axiom{primitivePart!(p)} replaces \axiom{p}  by its primitive part.
       mainContent : $ -> $
           ++ \axiom{mainContent(p)} returns the content of \axiom{p} viewed as a univariate 
           ++ polynomial in its main variable and with coefficients in the 
           ++ polynomial ring generated by its other variables over \axiom{R}.
       mainPrimitivePart : $ -> $
           ++ \axiom{mainPrimitivePart(p)} returns the primitive part of \axiom{p} viewed as a
           ++ univariate polynomial in its main variable and with coefficients 
           ++ in the polynomial ring generated by its other variables over \axiom{R}.
       mainSquareFreePart : $ -> $
           ++ \axiom{mainSquareFreePart(p)} returns the square free part of \axiom{p} viewed as a 
           ++ univariate polynomial in its main variable and with coefficients 
           ++ in the polynomial ring generated by its other variables over \axiom{R}.

 add
     O ==> OutputForm
     NNI ==> NonNegativeInteger
     INT ==> Integer

     exactQuo : (R,R) -> R

     coerce(p:$):O ==
       ground? (p) => (ground(p))::O
       if one?((ip := init(p)))
         then
           if zero?((tp := tail(p)))
             then
               if one?((dp := mdeg(p)))
                 then
                   return((mvar(p))::O)
                 else
                   return(((mvar(p))::O **$O (dp::O)))
             else
               if one?((dp := mdeg(p)))
                 then
                   return((mvar(p))::O +$O (tp::O))
                 else
                   return(((mvar(p))::O **$O (dp::O)) +$O (tp::O))
          else
           if zero?((tp := tail(p)))
             then
               if one?((dp := mdeg(p)))
                 then
                   return((ip::O) *$O  (mvar(p))::O)
                 else
                   return((ip::O) *$O ((mvar(p))::O **$O (dp::O)))
             else
               if one?(mdeg(p))
                 then
                   return(((ip::O) *$O  (mvar(p))::O) +$O (tp::O))
       ((ip)::O *$O ((mvar(p))::O **$O ((mdeg(p)::O))) +$O (tail(p)::O))

     mvar p ==
       ground?(p) => error"Error in mvar from RPOLCAT : #1 is constant."
       mainVariable(p)::V

     mdeg p == 
       ground?(p) => 0$NNI
       degree(p,mainVariable(p)::V)

     init p ==
       ground?(p) => error"Error in mvar from RPOLCAT : #1 is constant."
       v := mainVariable(p)::V
       coefficient(p,v,degree(p,v))

     leadingCoefficient (p,v) ==
       zero? (d := degree(p,v)) => p
       coefficient(p,v,d)

     head p ==
       ground? p => p
       v := mainVariable(p)::V
       d := degree(p,v)
       monomial(coefficient(p,v,d),v,d)

     reductum(p,v) ==
       zero? (d := degree(p,v)) => 0$$
       p - monomial(coefficient(p,v,d),v,d)

     tail p ==
       ground? p => 0$$
       p - head(p)

     deepestTail p ==
       ground? p => 0$$
       ground? tail(p) => tail(p)
       mvar(p) > mvar(tail(p)) => tail(p)
       deepestTail(tail(p))

     iteratedInitials p == 
       ground? p => []
       p := init(p)
       cons(p,iteratedInitials(p))

     localDeepestInitial (p : $) : $ == 
       ground? p => p
       localDeepestInitial init p

     deepestInitial p == 
       ground? p => error"Error in deepestInitial from RPOLCAT : #1 is constant."
       localDeepestInitial init p

     monic? p ==
       ground? p => false
       (recip(init(p))$$ case $)@Boolean

     quasiMonic?  p ==
       ground? p => false
       ground?(init(p))

     mainMonomial p == 
       zero? p => error"Error in mainMonomial from RPOLCAT : #1 is zero"
       ground? p => 1$$
       v := mainVariable(p)::V
       monomial(1$$,v,degree(p,v))

     leastMonomial p == 
       zero? p => error"Error in leastMonomial from RPOLCAT : #1 is zero"
       ground? p => 1$$
       v := mainVariable(p)::V
       monomial(1$$,v,minimumDegree(p,v))

     mainCoefficients p == 
       zero? p => error"Error in mainCoefficients from RPOLCAT : #1 is zero"
       ground? p => [p]
       v := mainVariable(p)::V
       coefficients(univariate(p,v)@SparseUnivariatePolynomial($))

     mainMonomials p == 
       zero? p => error"Error in mainMonomials from RPOLCAT : #1 is zero"
       ground? p => [1$$]
       v := mainVariable(p)::V
       lm := monomials(univariate(p,v)@SparseUnivariatePolynomial($))
       [monomial(1$$,v,degree(m)) for m in lm]

     RittWuCompare (a,b) ==
       (ground? b and  ground? a) => "failed"::Union(Boolean,"failed")
       ground? b => false::Union(Boolean,"failed")
       ground? a => true::Union(Boolean,"failed")
       mvar(a) < mvar(b) => true::Union(Boolean,"failed")
       mvar(a) > mvar(b) => false::Union(Boolean,"failed")
       mdeg(a) < mdeg(b) => true::Union(Boolean,"failed")
       mdeg(a) > mdeg(b) => false::Union(Boolean,"failed")
       lc := RittWuCompare(init(a),init(b))
       lc case Boolean => lc
       RittWuCompare(tail(a),tail(b))

     infRittWu? (a,b) ==
       lc : Union(Boolean,"failed") := RittWuCompare(a,b)
       lc case Boolean => lc::Boolean
       false
       
     supRittWu? (a,b) ==
       infRittWu? (b,a)

     prem (a:$, b:$)  : $ == 
       cP := lazyPremWithDefault (a,b)
       ((cP.coef) ** (cP.gap)) * cP.remainder

     pquo (a:$, b:$)  : $ == 
       cPS := lazyPseudoDivide (a,b)
       c := (cPS.coef) ** (cPS.gap)
       c * cPS.quotient

     prem (a:$, b:$, v:V) : $ ==
       cP := lazyPremWithDefault (a,b,v)
       ((cP.coef) ** (cP.gap)) * cP.remainder  

     pquo (a:$, b:$, v:V)  : $ == 
       cPS := lazyPseudoDivide (a,b,v)
       c := (cPS.coef) ** (cPS.gap)
       c * cPS.quotient     

     lazyPrem (a:$, b:$) : $ ==
       (not ground?(b)) and (monic?(b)) => monicModulo(a,b)
       (lazyPremWithDefault (a,b)).remainder
       
     lazyPquo (a:$, b:$) : $ ==
       (lazyPseudoDivide (a,b)).quotient

     lazyPrem (a:$, b:$, v:V) : $ ==
       zero? b =>  error"Error in lazyPrem : ($,$,V) -> $ from RPOLCAT : #2 is zero"
       ground?(b) => 0$$
       (v = mvar(b)) => lazyPrem(a,b)
       dbv : NNI := degree(b,v)
       zero? dbv => 0$$
       dav : NNI  := degree(a,v)
       zero? dav => a
       test : INT := dav::INT - dbv 
       lcbv : $ := leadingCoefficient(b,v)
       while not zero?(a) and not negative?(test) repeat
         lcav := leadingCoefficient(a,v)
         term := monomial(lcav,v,test::NNI)
         a := lcbv * a - term * b
         test := degree(a,v)::INT - dbv 
       a
         
     lazyPquo (a:$, b:$, v:V) : $ ==
       (lazyPseudoDivide (a,b,v)).quotient

     headReduce (a:$,b:$) == 
       ground? b => error"Error in headReduce : ($,$) -> Boolean from TSETCAT : #2 is constant"
       ground? a => a
       mvar(a) = mvar(b) => lazyPrem(a,b)
       while not reduced?((ha := head a),b) repeat
         lrc := lazyResidueClass(ha,b)
         if zero? tail(a)
           then
             a := lrc.polnum
           else
             a := lrc.polnum +  (lrc.polden)**(lrc.power) * tail(a)
       a

     initiallyReduce(a:$,b:$) ==
       ground? b => error"Error in initiallyReduce : ($,$) -> Boolean from TSETCAT : #2 is constant"
       ground? a => a
       v := mvar(b)
       mvar(a) = v => lazyPrem(a,b)
       ia := a
       ma := 1$$
       ta := 0$$
       while (not ground?(ia)) and (mvar(ia) >= mvar(b)) repeat
         if (mvar(ia) = mvar(b)) and (mdeg(ia) >= mdeg(b))
           then
             iamodb := lazyResidueClass(ia,b)
             ia := iamodb.polnum
             if not zero? ta
               then
                 ta :=  (iamodb.polden)**(iamodb.power) * ta
         if zero? ia 
           then 
             ia := ta
             ma := 1$$
             ta := 0$$
           else
             if not ground?(ia)
               then
                 ta := tail(ia) * ma + ta
                 ma := mainMonomial(ia) * ma
                 ia := init(ia)
       ia * ma + ta

     lazyPremWithDefault (a,b) == 
       ground?(b) => error"Error in lazyPremWithDefault from RPOLCAT : #2 is constant"
       ground?(a) => [1$$,0$NNI,a]
       xa := mvar a
       xb := mvar b
       xa < xb => [1$$,0$NNI,a]
       lcb : $ := init b 
       db : NNI := mdeg b
       test : INT := degree(a,xb)::INT - db
       delta : INT := max(test + 1$INT, 0$INT) 
       if xa = xb 
         then
           b := tail b
           while not zero?(a) and not negative?(test) repeat 
             term := monomial(init(a),xb,test::NNI)
             a := lcb * tail(a) - term * b 
             delta := delta - 1$INT 
             test := degree(a,xb)::INT - db
         else 
           while not zero?(a) and not negative?(test) repeat 
             term := monomial(leadingCoefficient(a,xb),xb,test::NNI)
             a := lcb * a - term * b
             delta := delta - 1$INT 
             test := degree(a,xb)::INT - db
       [lcb, (delta::NNI), a]

     lazyPremWithDefault (a,b,v) == 
       zero? b =>  error"Error in lazyPremWithDefault : ($,$,V) -> $ from RPOLCAT : #2 is zero"
       ground?(b) => [b,1$NNI,0$$]
       (v = mvar(b)) => lazyPremWithDefault(a,b)
       dbv : NNI := degree(b,v)
       zero? dbv => [b,1$NNI,0$$]
       dav : NNI  := degree(a,v)
       zero? dav => [1$$,0$NNI,a]
       test : INT := dav::INT - dbv 
       delta : INT := max(test + 1$INT, 0$INT) 
       lcbv : $ := leadingCoefficient(b,v)
       while not zero?(a) and not negative?(test) repeat
         lcav := leadingCoefficient(a,v)
         term := monomial(lcav,v,test::NNI)
         a := lcbv * a - term * b
         delta := delta - 1$INT 
         test := degree(a,v)::INT - dbv 
       [lcbv, (delta::NNI), a]

     pseudoDivide (a,b) == 
       cPS := lazyPseudoDivide (a,b)
       c := (cPS.coef) ** (cPS.gap)
       [c * cPS.quotient, c * cPS.remainder]

     lazyPseudoDivide (a,b) == 
       ground?(b) => error"Error in lazyPseudoDivide from RPOLCAT : #2 is constant"
       ground?(a) => [1$$,0$NNI,0$$,a]
       xa := mvar a 
       xb := mvar b
       xa < xb => [1$$,0$NNI,0$$, a]
       lcb : $ := init b 
       db : NNI := mdeg b
       q := 0$$
       test : INT := degree(a,xb)::INT - db
       delta : INT := max(test + 1$INT, 0$INT) 
       if xa = xb 
         then
           b := tail b
           while not zero?(a) and not negative?(test) repeat 
             term := monomial(init(a),xb,test::NNI)
             a := lcb * tail(a) - term * b 
             q := lcb * q + term
             delta := delta - 1$INT 
             test := degree(a,xb)::INT - db
         else 
           while not zero?(a) and not negative?(test) repeat 
             term := monomial(leadingCoefficient(a,xb),xb,test::NNI)
             a := lcb * a - term * b
             q := lcb * q + term
             delta := delta - 1$INT 
             test := degree(a,xb)::INT - db
       [lcb, (delta::NNI), q, a]

     lazyPseudoDivide (a,b,v) == 
       zero? b =>  error"Error in lazyPseudoDivide : ($,$,V) -> $ from RPOLCAT : #2 is zero"
       ground?(b) => [b,1$NNI,a,0$$]
       (v = mvar(b)) => lazyPseudoDivide(a,b)
       dbv : NNI := degree(b,v)
       zero? dbv => [b,1$NNI,a,0$$]
       dav : NNI  := degree(a,v)
       zero? dav => [1$$,0$NNI,0$$, a]
       test : INT := dav::INT - dbv 
       delta : INT := max(test + 1$INT, 0$INT) 
       lcbv : $ := leadingCoefficient(b,v)
       q := 0$$
       while not zero?(a) and not negative?(test) repeat
         lcav := leadingCoefficient(a,v)
         term := monomial(lcav,v,test::NNI)
         a := lcbv * a - term * b
         q := lcbv * q + term
         delta := delta - 1$INT 
         test := degree(a,v)::INT - dbv 
       [lcbv, (delta::NNI), q, a]

     monicModulo (a,b) == 
       ground?(b) => error"Error in monicModulo from RPOLCAT : #2 is constant"
       rec : Union($,"failed") 
       rec := recip((ib := init(b)))$$
       (rec case "failed")@Boolean => error"Error in monicModulo from RPOLCAT : #2 is not monic"
       ground? a => a
       ib * ((lazyPremWithDefault ((rec::$) * a,(rec::$) * b)).remainder)

     lazyResidueClass(a,b) ==
       zero? b => [a,1$$,0$NNI]
       ground? b => [0$$,1$$,0$NNI]
       ground? a => [a,1$$,0$NNI]
       xa := mvar a
       xb := mvar b
       xa < xb => [a,1$$,0$NNI]
       monic?(b) => [monicModulo(a,b),1$$,0$NNI]
       lcb : $ := init b 
       db : NNI := mdeg b
       test : INT := degree(a,xb)::INT - db
       pow : NNI := 0
       if xa = xb 
         then
           b := tail b
           while not zero?(a) and not negative?(test) repeat 
             term := monomial(init(a),xb,test::NNI)
             a := lcb * tail(a) - term * b 
             pow := pow + 1$NNI
             test := degree(a,xb)::INT - db
         else 
           while not zero?(a) and not negative?(test) repeat 
             term := monomial(leadingCoefficient(a,xb),xb,test::NNI)
             a := lcb * a - term * b
             pow := pow + 1$NNI
             test := degree(a,xb)::INT - db
       [a,lcb,pow]

     reduced? (a:$,b:$) : Boolean ==
       degree(a,mvar(b)) < mdeg(b)

     reduced? (p:$, lq : List($)) : Boolean ==
       ground? p => true
       while (not empty? lq) and (reduced?(p, first lq)) repeat
         lq := rest lq
       empty? lq

     headReduced? (a:$,b:$) : Boolean ==
       reduced?(head(a),b)

     headReduced? (p:$, lq : List($)) : Boolean ==
       reduced?(head(p),lq)

     initiallyReduced? (a:$,b:$) : Boolean ==
       ground? b => error"Error in initiallyReduced? : ($,$) -> Bool. from RPOLCAT : #2 is constant"
       ground?(a) => true
       mvar(a) < mvar(b) => true
       (mvar(a) = mvar(b)) => reduced?(a,b)
       initiallyReduced?(init(a),b)

     initiallyReduced? (p:$, lq : List($)) : Boolean ==
       ground? p => true
       while (not empty? lq) and (initiallyReduced?(p, first lq)) repeat
         lq := rest lq
       empty? lq

     normalized?(a:$,b:$) : Boolean ==
       ground? b => error"Error in  normalized? : ($,$) -> Boolean from TSETCAT : #2 is constant"
       ground? a => true
       mvar(a) < mvar(b) => true
       (mvar(a) = mvar(b)) => false
       normalized?(init(a),b)

     normalized? (p:$, lq : List($)) : Boolean ==
       while (not empty? lq) and (normalized?(p, first lq)) repeat
         lq := rest lq
       empty? lq       

     if R has IntegralDomain
     then

       exactQuo(r:R,s:R):R ==
         R has EuclideanDomain => r quo$R s
         (r exquo$R s)::R

       exactQuotient (p:$,r:R) ==
         (p exquo$$ r)::$

       exactQuotient (a:$,b:$) ==
         ground? b => exactQuotient(a,ground(b))
         (a exquo$$ b)::$

       exactQuotient! (a:$,b:$) ==
         ground? b => exactQuotient!(a,ground(b))
         a := (a exquo$$ b)::$

       if (R has GcdDomain) and not(R has Field)
       then

         primPartElseUnitCanonical p ==
           primitivePart p

         primitivePart! p ==
           zero? p => p
           if one?(cp := content(p))
             then
               p := unitCanonical p
             else
               p := unitCanonical exactQuotient!(p,cp) 
           p

         primPartElseUnitCanonical! p ==
           primitivePart! p

       else
         primPartElseUnitCanonical p ==
           unitCanonical p

         primPartElseUnitCanonical! p ==
           p := unitCanonical p


     if R has GcdDomain
     then

       gcd(r:R,p:$):R ==
         one? r => r
         zero? p => r
         ground? p => gcd(r,ground(p))$R
         gcd(gcd(r,init(p)),tail(p))

       mainContent p ==
         zero? p => p
         "gcd"/mainCoefficients(p)

       mainPrimitivePart p ==
         zero? p => p
         (unitNormal((p exquo$$ mainContent(p))::$)).canonical

       mainSquareFreePart p ==
         ground? p => p
         v := mainVariable(p)::V
         sfp : SparseUnivariatePolynomial($)
         sfp := squareFreePart(univariate(p,v)@SparseUnivariatePolynomial($))
         multivariate(sfp,v)

     if (V has ConvertibleTo(Symbol))
       then

         PR ==> Polynomial R
         PQ ==> Polynomial Fraction Integer
         PZ ==> Polynomial Integer
         IES ==> IndexedExponents(Symbol)
         Q ==> Fraction Integer
         Z ==> Integer

         convert(p:$) : PR ==
           ground? p => (ground(p)$$)::PR
           v : V := mvar(p)
           d : NNI := mdeg(p)
           convert(init(p))@PR *$PR ((convert(v)@Symbol)::PR)**d +$PR convert(tail(p))@PR

         coerce(p:$) : PR ==
           convert(p)@PR

         localRetract : PR -> $
         localRetractPQ : PQ -> $
         localRetractPZ : PZ -> $
         localRetractIfCan : PR -> Union($,"failed")
         localRetractIfCanPQ : PQ -> Union($,"failed")
         localRetractIfCanPZ : PZ -> Union($,"failed")

         if V has Finite
           then 

             sizeV : NNI := size()$V
             lv : List Symbol
             lv := [convert(index(i::PositiveInteger)$V)@Symbol for i in 1..sizeV]

             localRetract(p : PR) : $ ==
               ground? p => (ground(p)$PR)::$
               mvp : Symbol := (mainVariable(p)$PR)::Symbol
               d : NNI
               imvp : PositiveInteger := (position(mvp,lv)$(List Symbol))::PositiveInteger 
               vimvp : V := index(imvp)$V
               xvimvp,c : $ 
               newp := 0$$
               while (not zero? (d := degree(p,mvp))) repeat
                 c := localRetract(coefficient(p,mvp,d)$PR)
                 xvimvp := monomial(c,vimvp,d)$$
                 newp := newp +$$ xvimvp
                 p := p -$PR monomial(coefficient(p,mvp,d)$PR,mvp,d)$PR
               newp +$$ localRetract(p)

             if R has Algebra Fraction Integer
               then 
                 localRetractPQ(pq:PQ):$ ==
                   ground? pq => ((ground(pq)$PQ)::R)::$
                   mvp : Symbol := (mainVariable(pq)$PQ)::Symbol
                   d : NNI
                   imvp : PositiveInteger := (position(mvp,lv)$(List Symbol))::PositiveInteger 
                   vimvp : V := index(imvp)$V
                   xvimvp,c : $ 
                   newp := 0$$
                   while (not zero? (d := degree(pq,mvp))) repeat
                     c := localRetractPQ(coefficient(pq,mvp,d)$PQ)
                     xvimvp := monomial(c,vimvp,d)$$
                     newp := newp +$$ xvimvp
                     pq := pq -$PQ monomial(coefficient(pq,mvp,d)$PQ,mvp,d)$PQ
                   newp +$$ localRetractPQ(pq)

             if R has Algebra Integer
               then 
                 localRetractPZ(pz:PZ):$ ==
                   ground? pz => ((ground(pz)$PZ)::R)::$
                   mvp : Symbol := (mainVariable(pz)$PZ)::Symbol
                   d : NNI
                   imvp : PositiveInteger := (position(mvp,lv)$(List Symbol))::PositiveInteger 
                   vimvp : V := index(imvp)$V
                   xvimvp,c : $ 
                   newp := 0$$
                   while (not zero? (d := degree(pz,mvp))) repeat
                     c := localRetractPZ(coefficient(pz,mvp,d)$PZ)
                     xvimvp := monomial(c,vimvp,d)$$
                     newp := newp +$$ xvimvp
                     pz := pz -$PZ monomial(coefficient(pz,mvp,d)$PZ,mvp,d)$PZ
                   newp +$$ localRetractPZ(pz)

             retractable?(p:PR):Boolean ==
               lvp := variables(p)$PR
               while not empty? lvp and member?(first lvp,lv) repeat
                 lvp := rest lvp
               empty? lvp   
                     
             retractablePQ?(p:PQ):Boolean ==
               lvp := variables(p)$PQ
               while not empty? lvp and member?(first lvp,lv) repeat
                 lvp := rest lvp
               empty? lvp       
                 
             retractablePZ?(p:PZ):Boolean ==
               lvp := variables(p)$PZ
               while not empty? lvp and member?(first lvp,lv) repeat
                 lvp := rest lvp
               empty? lvp                        

             localRetractIfCan(p : PR): Union($,"failed") ==
               not retractable?(p) => "failed"::Union($,"failed")
               localRetract(p)::Union($,"failed")

             localRetractIfCanPQ(p : PQ): Union($,"failed") ==
               not retractablePQ?(p) => "failed"::Union($,"failed")
               localRetractPQ(p)::Union($,"failed")

             localRetractIfCanPZ(p : PZ): Union($,"failed") ==
               not retractablePZ?(p) => "failed"::Union($,"failed")
               localRetractPZ(p)::Union($,"failed")

         if R has Algebra Fraction Integer
           then 

             mpc2Z := MPolyCatFunctions2(Symbol,IES,IES,Z,R,PZ,PR)
             mpc2Q := MPolyCatFunctions2(Symbol,IES,IES,Q,R,PQ,PR)

             retract(pz:PZ) ==
               rif : Union($,"failed") := retractIfCan(pz)@Union($,"failed")
               (rif case "failed") => error"failed in retract: POLY Z -> $ from RPOLCAT"
               rif::$

             convert(pz:PZ) ==
               retract(pz)@$

             retract(pq:PQ) ==
               rif : Union($,"failed") := retractIfCan(pq)@Union($,"failed")
               (rif case "failed") => error"failed in retract: POLY Z -> $ from RPOLCAT"
               rif::$

             convert(pq:PQ) ==
               retract(pq)@$

             if not (R has QuotientFieldCategory(Integer))
               then
                 -- the only operation to implement is retractIfCan : PR -> Union($,"failed")
                 -- when V does not have Finite

                 if V has Finite
                   then
                     retractIfCan(pr:PR) ==
                       localRetractIfCan(pr)@Union($,"failed")

                     retractIfCan(pq:PQ) ==
                       localRetractIfCanPQ(pq)@Union($,"failed")
                   else
                     retractIfCan(pq:PQ) ==
                       pr : PR := map(#1::R,pq)$mpc2Q
                       retractIfCan(pr)@Union($,"failed")

                 retractIfCan(pz:PZ) ==
                   pr : PR := map(#1::R,pz)$mpc2Z
                   retractIfCan(pr)@Union($,"failed")

                 retract(pr:PR) ==
                   rif : Union($,"failed") := retractIfCan(pr)@Union($,"failed")
                   (rif case "failed") => error"failed in retract: POLY Z -> $ from RPOLCAT"
                   rif::$

                 convert(pr:PR) ==
                   retract(pr)@$

               else
                 -- the only operation to implement is retractIfCan : PQ -> Union($,"failed")
                 -- when V does not have Finite
                 mpc2ZQ := MPolyCatFunctions2(Symbol,IES,IES,Z,Q,PZ,PQ)
                 mpc2RQ := MPolyCatFunctions2(Symbol,IES,IES,R,Q,PR,PQ)
                 ZToQ(z:Z):Q == coerce(z)@Q
                 RToQ(r:R):Q == retract(r)@Q

                 PZToPQ (pz:PZ):PQ == map(ZToQ,pz)$mpc2ZQ
                 PRToPQ (pr:PR):PQ == map(RToQ,pr)$mpc2RQ

                 retractIfCan(pz:PZ) ==
                   pq : PQ := PZToPQ(pz)
                   retractIfCan(pq)@Union($,"failed")

                 if V has Finite
                   then
                     retractIfCan(pq:PQ) ==
                       localRetractIfCanPQ(pq)@Union($,"failed")

                     convert(pr:PR) ==
                       lrif : Union($,"failed") := localRetractIfCan(pr)@Union($,"failed")
                       (lrif case "failed") => error"failed in convert: PR->$ from RPOLCAT"
                       lrif::$
                   else
                     convert(pr:PR) ==
                       pq : PQ := PRToPQ(pr)
                       retract(pq)@$

         if (R has Algebra Integer) and not(R has Algebra Fraction Integer)
           then 

             mpc2Z := MPolyCatFunctions2(Symbol,IES,IES,Z,R,PZ,PR)

             retract(pz:PZ) ==
               rif : Union($,"failed") := retractIfCan(pz)@Union($,"failed")
               (rif case "failed") => error"failed in retract: POLY Z -> $ from RPOLCAT"
               rif::$

             convert(pz:PZ) ==
               retract(pz)@$

             if not (R has IntegerNumberSystem)
               then
                 -- the only operation to implement is retractIfCan : PR -> Union($,"failed")
                 -- when V does not have Finite

                 if V has Finite
                   then
                     retractIfCan(pr:PR) ==
                       localRetractIfCan(pr)@Union($,"failed")

                     retractIfCan(pz:PZ) ==
                       localRetractIfCanPZ(pz)@Union($,"failed")
                   else
                     retractIfCan(pz:PZ) ==
                       pr : PR := map(#1::R,pz)$mpc2Z
                       retractIfCan(pr)@Union($,"failed")

                 retract(pr:PR) ==
                   rif : Union($,"failed") := retractIfCan(pr)@Union($,"failed")
                   (rif case "failed") => error"failed in retract: POLY Z -> $ from RPOLCAT"
                   rif::$

                 convert(pr:PR) ==
                   retract(pr)@$

               else
                 -- the only operation to implement is retractIfCan : PZ -> Union($,"failed")
                 -- when V does not have Finite

                 mpc2RZ := MPolyCatFunctions2(Symbol,IES,IES,R,Z,PR,PZ)
                 RToZ(r:R):Z == retract(r)@Z
                 PRToPZ (pr:PR):PZ == map(RToZ,pr)$mpc2RZ

                 if V has Finite
                   then
                     convert(pr:PR) ==
                       lrif : Union($,"failed") := localRetractIfCan(pr)@Union($,"failed")
                       (lrif case "failed") => error"failed in convert: PR->$ from RPOLCAT"
                       lrif::$
                     retractIfCan(pz:PZ) ==
                       localRetractIfCanPZ(pz)@Union($,"failed")
                   else
                     convert(pr:PR) ==
                       pz : PZ := PRToPZ(pr)
                       retract(pz)@$


         if not(R has Algebra Integer) and not(R has Algebra Fraction Integer)
           then 
             -- the only operation to implement is retractIfCan : PR -> Union($,"failed")

             if V has Finite
               then
                 retractIfCan(pr:PR) ==
                   localRetractIfCan(pr)@Union($,"failed")

             retract(pr:PR) ==
               rif : Union($,"failed") := retractIfCan(pr)@Union($,"failed")
               (rif case "failed") => error"failed in retract: POLY Z -> $ from RPOLCAT"
               rif::$

             convert(pr:PR) ==
               retract(pr)@$

         if (R has RetractableTo(INT))
           then

             convert(pol:$):String ==
               ground?(pol) => string(retract(ground(pol))@INT)
               ipol : $ := init(pol)
               vpol : V := mvar(pol)
               dpol : NNI := mdeg(pol)
               tpol: $  := tail(pol)
               sipol,svpol,sdpol,stpol : String
               if one? ipol
                 then 
                   sipol := empty()$String
                 else
                   if one?(-ipol)
                     then 
                       sipol := "-"
                     else
                       sipol := convert(ipol)@String
                       if not monomial?(ipol)
                         then
                           sipol := concat(["(",sipol,")*"])$String
                         else 
                           sipol := concat(sipol,"*")$String
               svpol := string(convert(vpol)@Symbol)
               if one? dpol
                 then
                   sdpol :=  empty()$String
                 else
                   sdpol := concat("**",string(convert(dpol)@INT))$String 
               if zero? tpol
                 then
                   stpol :=  empty()$String
                 else
                   if ground?(tpol)
                     then
                       n := retract(ground(tpol))@INT
                       if positive? n
                         then
                           stpol :=  concat(" +", string n)$String
                         else
                           stpol := string n
                     else
                       stpol := convert(tpol)@String
                       if not member?((stpol.1)::String,["+","-"])$(List String)
                         then
                           stpol :=  concat(" + ",stpol)$String
               concat([sipol,svpol,sdpol,stpol])$String

@
Based on the {\bf PseudoRemainderSequence} package, the domain
constructor {\bf NewSparseMulitvariatePolynomial} extends
the constructor {\bf SparseMultivariatePolynomial}. It also provides
some additional operations related to polynomial system solving
by means of triangular sets.
\section{domain NSMP NewSparseMultivariatePolynomial}
<<domain NSMP NewSparseMultivariatePolynomial>>=
)abbrev domain NSMP NewSparseMultivariatePolynomial
++ Author: Marc Moreno Maza
++ Date Created: 22/04/94
++ Date Last Updated: 14/12/1998
++ Basic Operations: 
++ Related Domains:
++ Also See:
++ AMS Classifications:
++ Keywords: 
++ Examples:
++ References:
++ Description: A post-facto extension for \axiomType{SMP} in order
++ to speed up operations related to pseudo-division and gcd.
++ This domain is based on the \axiomType{NSUP} constructor which is 
++ itself a post-facto extension of the \axiomType{SUP} constructor.
++ Version: 2

NewSparseMultivariatePolynomial(R,VarSet) : Exports == Implementation where
  R:Ring
  VarSet:OrderedSet
  N ==> NonNegativeInteger
  Z ==> Integer
  SUP ==> NewSparseUnivariatePolynomial
  SMPR ==> SparseMultivariatePolynomial(R, VarSet)
  SUP2 ==> NewSparseUnivariatePolynomialFunctions2($,$)

  Exports == Join(RecursivePolynomialCategory(R,IndexedExponents VarSet, VarSet),
                  CoercibleTo(SMPR),RetractableTo(SMPR))

  Implementation ==  SparseMultivariatePolynomial(R, VarSet) add

     D := NewSparseUnivariatePolynomial($)
     VPoly:=  Record(v:VarSet,ts:D)
     Rep:= Union(R,VPoly)

    --local function
     PSimp: (D,VarSet) -> %

     PSimp(up,mv) ==
       if degree(up) = 0 then leadingCoefficient(up) else [mv,up]$VPoly

     coerce (p:$):SMPR == 
       p pretend SMPR

     coerce (p:SMPR):$ == 
       p pretend $

     retractIfCan (p:$) : Union(SMPR,"failed") == 
       (p pretend SMPR)::Union(SMPR,"failed")

     mvar p == 
       p case R => error" Error in mvar from NSMP : #1 has no variables."
       p.v

     mdeg p == 
       p case R => 0$N
       degree(p.ts)$D

     init p == 
       p case R => error" Error in init from NSMP : #1 has no variables."
       leadingCoefficient(p.ts)$D

     head p == 
       p case R => p
       ([p.v,leadingMonomial(p.ts)$D]$VPoly)::Rep

     tail p == 
       p case R => 0$$
       red := reductum(p.ts)$D
       ground?(red)$D => (ground(red)$D)::Rep
       ([p.v,red]$VPoly)::Rep

     iteratedInitials p == 
       p case R => [] 
       p := leadingCoefficient(p.ts)$D
       cons(p,iteratedInitials(p)) 

     localDeepestInitial (p : $) : $ == 
       p case R => p 
       localDeepestInitial leadingCoefficient(p.ts)$D 

     deepestInitial p == 
       p case R => error"Error in deepestInitial from NSMP : #1 has no variables."
       localDeepestInitial leadingCoefficient(p.ts)$D  

     mainMonomial p == 
       zero? p => error"Error in mainMonomial from NSMP : the argument is zero"
       p case R => 1$$ 
       monomial(1$$,p.v,degree(p.ts)$D)

     leastMonomial p == 
       zero? p => error"Error in leastMonomial from NSMP : the argument is zero"
       p case R => 1$$
       monomial(1$$,p.v,minimumDegree(p.ts)$D)

     mainCoefficients p == 
       zero? p => error"Error in mainCoefficients from NSMP : the argument is zero"
       p case R => [p]
       coefficients(p.ts)$D

     leadingCoefficient(p:$,x:VarSet):$ == 
       (p case R) => p
       p.v = x => leadingCoefficient(p.ts)$D
       zero? (d := degree(p,x)) => p
       coefficient(p,x,d)

     localMonicModulo(a:$,b:$):$ ==
       -- b is assumed to have initial 1
       a case R => a
       a.v < b.v => a 
       mM: $
       if a.v > b.v
         then 
           m : D := map(localMonicModulo(#1,b),a.ts)$SUP2
         else
           m : D := monicModulo(a.ts,b.ts)$D
       if ground?(m)$D 
          then 
            mM := (ground(m)$D)::Rep 
          else 
            mM := ([a.v,m]$VPoly)::Rep
       mM

     monicModulo (a,b) == 
       b case R => error"Error in monicModulo from NSMP : #2 is constant"
       ib : $ := init(b)@$
       not ground?(ib)$$ => 
         error"Error in monicModulo from NSMP : #2 is not monic"
       mM : $
       if not one?(ib)$$
         then
           r : R := ground(ib)$$
           rec : Union(R,"failed"):= recip(r)$R
           (rec case "failed") =>
             error"Error in monicModulo from NSMP : #2 is not monic"
           a case R => a
           a := (rec::R) * a
           b := (rec::R) * b
           mM := ib * localMonicModulo (a,b)
         else
           mM := localMonicModulo (a,b)
       mM

     prem(a:$, b:$): $ == 
       -- with pseudoRemainder$NSUP
       b case R =>
         error "in prem$NSMP: ground? #2"
       db: N := degree(b.ts)$D
       lcb: $ := leadingCoefficient(b.ts)$D
       test: Z := degree(a,b.v)::Z - db
       delta: Z := max(test + 1$Z, 0$Z)
       (a case R) or (a.v < b.v) => lcb ** (delta::N) * a
       a.v = b.v =>
         r: D := pseudoRemainder(a.ts,b.ts)$D
         ground?(r) => return (ground(r)$D)::Rep 
         ([a.v,r]$VPoly)::Rep
       while not zero?(a) and not negative?(test) repeat 
         term := monomial(leadingCoefficient(a,b.v),b.v,test::N)
         a := lcb * a - term * b
         delta := delta - 1$Z 
         test := degree(a,b.v)::Z - db
       lcb ** (delta::N) * a

     pquo (a:$, b:$)  : $ == 
       cPS := lazyPseudoDivide (a,b)
       c := (cPS.coef) ** (cPS.gap)
       c * cPS.quotient

     pseudoDivide(a:$, b:$): Record (quotient : $, remainder : $) ==
       -- from RPOLCAT
       cPS := lazyPseudoDivide(a,b)
       c := (cPS.coef) ** (cPS.gap)
       [c * cPS.quotient, c * cPS.remainder]

     lazyPrem(a:$, b:$): $ == 
       -- with lazyPseudoRemainder$NSUP
       -- Uses leadingCoefficient: ($, V) -> $
       b case R =>
         error "in lazyPrem$NSMP: ground? #2"
       (a case R) or (a.v < b.v) =>  a
       a.v = b.v => PSimp(lazyPseudoRemainder(a.ts,b.ts)$D,a.v)
       db: N := degree(b.ts)$D
       lcb: $ := leadingCoefficient(b.ts)$D
       test: Z := degree(a,b.v)::Z - db
       while not zero?(a) and not negative?(test) repeat 
         term := monomial(leadingCoefficient(a,b.v),b.v,test::N)
         a := lcb * a - term * b
         test := degree(a,b.v)::Z - db
       a

     lazyPquo (a:$, b:$) : $ ==
       -- with lazyPseudoQuotient$NSUP
       b case R =>
         error " in lazyPquo$NSMP: #2 is conctant"
       (a case R) or (a.v < b.v) => 0
       a.v = b.v => PSimp(lazyPseudoQuotient(a.ts,b.ts)$D,a.v)
       db: N := degree(b.ts)$D
       lcb: $ := leadingCoefficient(b.ts)$D
       test: Z := degree(a,b.v)::Z - db
       q := 0$$
       test: Z := degree(a,b.v)::Z - db
       while not zero?(a) and not negative?(test) repeat 
         term := monomial(leadingCoefficient(a,b.v),b.v,test::N)
         a := lcb * a - term * b
         q := lcb * q + term
         test := degree(a,b.v)::Z - db
       q

     lazyPseudoDivide(a:$, b:$): Record(coef:$, gap: N,quotient:$, remainder:$) == 
       -- with lazyPseudoDivide$NSUP
       b case R =>
         error " in lazyPseudoDivide$NSMP: #2 is conctant"
       (a case R) or (a.v < b.v) => [1$$,0$N,0$$,a]
       a.v = b.v =>
         cgqr := lazyPseudoDivide(a.ts,b.ts)
         [cgqr.coef, cgqr.gap, PSimp(cgqr.quotient,a.v), PSimp(cgqr.remainder,a.v)]
       db: N := degree(b.ts)$D
       lcb: $ := leadingCoefficient(b.ts)$D
       test: Z := degree(a,b.v)::Z - db
       q := 0$$
       delta: Z := max(test + 1$Z, 0$Z) 
       while not zero?(a) and not negative?(test) repeat 
         term := monomial(leadingCoefficient(a,b.v),b.v,test::N)
         a := lcb * a - term * b
         q := lcb * q + term
         delta := delta - 1$Z 
         test := degree(a,b.v)::Z - db
       [lcb, (delta::N), q, a]

     lazyResidueClass(a:$, b:$): Record(polnum:$, polden:$, power:N) == 
       -- with lazyResidueClass$NSUP
       b case R =>
         error " in lazyResidueClass$NSMP: #2 is conctant"
       lcb: $ := leadingCoefficient(b.ts)$D
       (a case R) or (a.v < b.v) => [a,lcb,0]
       a.v = b.v =>
         lrc := lazyResidueClass(a.ts,b.ts)$D
         [PSimp(lrc.polnum,a.v), lrc.polden, lrc.power]
       db: N := degree(b.ts)$D
       test: Z := degree(a,b.v)::Z - db
       pow: N := 0
       while not zero?(a) and not negative?(test) repeat 
         term := monomial(leadingCoefficient(a,b.v),b.v,test::N)
         a := lcb * a - term * b
         pow := pow + 1
         test := degree(a,b.v)::Z - db
       [a, lcb, pow]

     if R has IntegralDomain
     then

       packD := PseudoRemainderSequence($,D)

       exactQuo(x:$, y:$):$ == 
         ex: Union($,"failed") := x exquo$$ y
         (ex case $) => ex::$
         error "in exactQuotient$NSMP: bad args"

       LazardQuotient(x:$, y:$, n: N):$ == 
         zero?(n) => error("LazardQuotient$NSMP : n = 0")
         one?(n) => x
         a: N := 1
         while n >= (b := 2*a) repeat a := b
         c: $ := x
         n := (n - a)::N
         repeat       
           one?(a) => return c
           a := a quo 2
           c := exactQuo(c*c,y)
           if n >= a then ( c := exactQuo(c*x,y) ; n := (n - a)::N )

       LazardQuotient2(p:$, a:$, b:$, n: N) ==
         zero?(n) => error " in LazardQuotient2$NSMP: bad #4"
         one?(n) => p
         c: $  := LazardQuotient(a,b,(n-1)::N)
         exactQuo(c*p,b)

       next_subResultant2(p:$, q:$, z:$, s:$) ==
         PSimp(next_sousResultant2(p.ts,q.ts,z.ts,s)$packD,p.v)

       subResultantGcd(a:$, b:$): $ ==
         (a case R) or (b case R) => 
           error "subResultantGcd$NSMP: one arg is constant"
         a.v ~= b.v => 
           error "subResultantGcd$NSMP: mvar(#1) ~= mvar(#2)"
         PSimp(subResultantGcd(a.ts,b.ts),a.v)

       halfExtendedSubResultantGcd1(a:$,b:$): Record (gcd: $, coef1: $) ==
         (a case R) or (b case R) => 
           error "halfExtendedSubResultantGcd1$NSMP: one arg is constant"
         a.v ~= b.v => 
           error "halfExtendedSubResultantGcd1$NSMP: mvar(#1) ~= mvar(#2)"
         hesrg := halfExtendedSubResultantGcd1(a.ts,b.ts)$D
         [PSimp(hesrg.gcd,a.v), PSimp(hesrg.coef1,a.v)]

       halfExtendedSubResultantGcd2(a:$,b:$): Record (gcd: $, coef2: $) ==
         (a case R) or (b case R) => 
           error "halfExtendedSubResultantGcd2$NSMP: one arg is constant"
         a.v ~= b.v => 
           error "halfExtendedSubResultantGcd2$NSMP: mvar(#1) ~= mvar(#2)"
         hesrg := halfExtendedSubResultantGcd2(a.ts,b.ts)$D
         [PSimp(hesrg.gcd,a.v), PSimp(hesrg.coef2,a.v)]

       extendedSubResultantGcd(a:$,b:$): Record (gcd: $, coef1: $, coef2: $) ==
         (a case R) or (b case R) => 
           error "extendedSubResultantGcd$NSMP: one arg is constant"
         a.v ~= b.v => 
           error "extendedSubResultantGcd$NSMP: mvar(#1) ~= mvar(#2)"
         esrg := extendedSubResultantGcd(a.ts,b.ts)$D
         [PSimp(esrg.gcd,a.v),PSimp(esrg.coef1,a.v),PSimp(esrg.coef2,a.v)]  

       resultant(a:$, b:$): $ ==
         (a case R) or (b case R) => 
           error "resultant$NSMP: one arg is constant"
         a.v ~= b.v => 
           error "resultant$NSMP: mvar(#1) ~= mvar(#2)"
         resultant(a.ts,b.ts)$D

       subResultantChain(a:$, b:$): List $ ==
         (a case R) or (b case R) => 
           error "subResultantChain$NSMP: one arg is constant"
         a.v ~= b.v => 
           error "subResultantChain$NSMP: mvar(#1) ~= mvar(#2)"
         [PSimp(up,a.v) for up in subResultantsChain(a.ts,b.ts)]

       lastSubResultant(a:$, b:$): $ ==
         (a case R) or (b case R) => 
           error "lastSubResultant$NSMP: one arg is constant"
         a.v ~= b.v => 
           error "lastSubResultant$NSMP: mvar(#1) ~= mvar(#2)"
         PSimp(lastSubResultant(a.ts,b.ts),a.v)

       if R has EuclideanDomain
       then

         exactQuotient (a:$,b:R) ==
           one? b => a
           a case R => (a::R quo$R b)::$
           ([a.v, map(exactQuotient(#1,b),a.ts)$SUP2]$VPoly)::Rep

         exactQuotient! (a:$,b:R) ==
           one? b => a
           a case R => (a::R quo$R b)::$
           a.ts := map(exactQuotient!(#1,b),a.ts)$SUP2
           a

       else

         exactQuotient (a:$,b:R) ==
           one? b => a
           a case R => ((a::R exquo$R b)::R)::$
           ([a.v, map(exactQuotient(#1,b),a.ts)$SUP2]$VPoly)::Rep

         exactQuotient! (a:$,b:R) == 
           one? b => a
           a case R => ((a::R exquo$R b)::R)::$
           a.ts := map(exactQuotient!(#1,b),a.ts)$SUP2
           a

     if R has GcdDomain
     then

       localGcd(r:R,p:$):R ==
         p case R => gcd(r,p::R)$R
         gcd(r,content(p))$R         

       gcd(r:R,p:$):R ==
         one? r => r
         zero? p => r
         localGcd(r,p)

       content p ==
         p case R => p
         up : D := p.ts
         r := 0$R
         while (not zero? up) and (not one? r) repeat
           r := localGcd(r,leadingCoefficient(up))
           up := reductum up
         r

       primitivePart! p ==
         zero? p => p
         p case R => 1$$
         cp := content(p)
         p.ts := unitCanonical(map(exactQuotient!(#1,cp),p.ts)$SUP2)$D
         p

@
\section{License}
<<license>>=
--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.
@
<<*>>=
<<license>>

<<domain NSUP NewSparseUnivariatePolynomial>>
<<package NSUP2 NewSparseUnivariatePolynomialFunctions2>>
<<category RPOLCAT RecursivePolynomialCategory>>
<<domain NSMP NewSparseMultivariatePolynomial>>
@
\eject
\begin{thebibliography}{99}
\bibitem{1} nothing
\end{thebibliography}
\end{document}