Initial population.

1 files changed, 1888 insertions, 0 deletions

diff --git a/src/algebra/newpoly.spad.pamphlet b/src/algebra/newpoly.spad.pamphlet
new file mode 100644
index 00000000..db2974e2
--- /dev/null
+++ b/src/algebra/newpoly.spad.pamphlet
@@ -0,0 +1,1888 @@
+\documentclass{article}
+\usepackage{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) => 
+                not ((yy.first.c) = 1) => 
+		   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 "failed" 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 "failed" 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)
+                (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 "failed" 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 "failed" 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 "failed" 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 "failed" 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 "failed" 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)))
+       if (((ip := init(p))) = 1)
+         then
+           if zero?((tp := tail(p)))
+             then
+--               if one?((dp := mdeg(p)))
+               if (((dp := mdeg(p))) = 1)
+                 then
+                   return((mvar(p))::O)
+                 else
+                   return(((mvar(p))::O **$O (dp::O)))
+             else
+--               if one?((dp := mdeg(p)))
+               if (((dp := mdeg(p))) = 1)
+                 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)))
+               if (((dp := mdeg(p))) = 1)
+                 then
+                   return((ip::O) *$O  (mvar(p))::O)
+                 else
+                   return((ip::O) *$O ((mvar(p))::O **$O (dp::O)))
+             else
+--               if one?(mdeg(p))
+               if ((mdeg(p)) = 1)
+                 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
+
+       if R has EuclideanDomain
+         then
+           exactQuo(r:R,s:R):R ==
+             r quo$R s
+         else
+           exactQuo(r:R,s:R):R ==
+             (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))
+           if ((cp := content(p)) = 1)
+             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
+         (r = 1) => 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)
+             ZToR (z:Z):R == coerce(z)@R
+             QToR (q:Q):R == coerce(q)@R
+             PZToPR (pz:PZ):PR == map(ZToR,pz)$mpc2Z
+             PQToPR (pq:PQ):PR == map(QToR,pq)$mpc2Q
+
+             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 := PQToPR(pq)
+                       retractIfCan(pr)@Union($,"failed")
+
+                 retractIfCan(pz:PZ) ==
+                   pr : PR := PZToPR(pz)
+                   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)
+             ZToR (z:Z):R == coerce(z)@R
+             PZToPR (pz:PZ):PR == map(ZToR,pz)$mpc2Z
+
+             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 := PZToPR(pz)
+                       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) => convert(retract(ground(pol))@INT)@String
+               ipol : $ := init(pol)
+               vpol : V := mvar(pol)
+               dpol : NNI := mdeg(pol)
+               tpol: $  := tail(pol)
+               sipol,svpol,sdpol,stpol : String
+--               if one? ipol
+               if (ipol = 1)
+                 then 
+                   sipol := empty()$String
+                 else
+--                   if one?(-ipol)
+                   if ((-ipol) = 1)
+                     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
+               if (dpol = 1)
+                 then
+                   sdpol :=  empty()$String
+                 else
+                   sdpol := concat("**",convert(convert(dpol)@INT)@String )$String 
+               if zero? tpol
+                 then
+                   stpol :=  empty()$String
+                 else
+                   if ground?(tpol)
+                     then
+                       n := retract(ground(tpol))@INT
+                       if n > 0
+                         then
+                           stpol :=  concat(" +",convert(n)@String)$String
+                         else
+                           stpol := convert(n)@String
+                     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)$$
+       if not ((ib) = 1)$$
+         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
+         (n = 1) => x
+         a: N := 1
+         while n >= (b := 2*a) repeat a := b
+         c: $ := x
+         n := (n - a)::N
+         repeat       
+--           one?(a) => return c
+           (a = 1) => 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
+         (n = 1) => 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
+           (b = 1) => 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
+           (b = 1) => 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
+           (b = 1) => 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
+           (b = 1) => 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
+         (r = 1) => 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
+         while (not zero? up) and (not (r = 1)) 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}