diff options
Diffstat (limited to 'src/algebra/newpoly.spad.pamphlet')
-rw-r--r-- | src/algebra/newpoly.spad.pamphlet | 1888 |
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} |