diff options
Diffstat (limited to 'src/algebra/poly.spad.pamphlet')
| -rw-r--r-- | src/algebra/poly.spad.pamphlet | 1249 |
1 files changed, 1249 insertions, 0 deletions
diff --git a/src/algebra/poly.spad.pamphlet b/src/algebra/poly.spad.pamphlet new file mode 100644 index 00000000..ad0206cd --- /dev/null +++ b/src/algebra/poly.spad.pamphlet @@ -0,0 +1,1249 @@ +\documentclass{article} +\usepackage{axiom} +\begin{document} +\title{\$SPAD/src/algebra poly.spad} +\author{Dave Barton, James Davenport, Barry Trager, Patrizia Gianni, Marc Moreno Maza} +\maketitle +\begin{abstract} +\end{abstract} +\eject +\tableofcontents +\eject +\section{domain FM FreeModule} +<<domain FM FreeModule>>= +)abbrev domain FM FreeModule +++ Author: Dave Barton, James Davenport, Barry Trager +++ Date Created: +++ Date Last Updated: +++ Basic Functions: BiModule(R,R) +++ Related Constructors: +++ Also See: +++ AMS Classifications: +++ Keywords: +++ References: +++ Description: +++ A bi-module is a free module +++ over a ring with generators indexed by an ordered set. +++ Each element can be expressed as a finite linear combination of +++ generators. Only non-zero terms are stored. + +FreeModule(R:Ring,S:OrderedSet): + Join(BiModule(R,R),IndexedDirectProductCategory(R,S)) with + if R has CommutativeRing then Module(R) + == IndexedDirectProductAbelianGroup(R,S) add + --representations + Term:= Record(k:S,c:R) + Rep:= List Term + --declarations + x,y: % + r: R + n: Integer + f: R -> R + s: S + --define + if R has EntireRing then + r * x == + zero? r => 0 +-- one? r => x + (r = 1) => x + --map(r*#1,x) + [[u.k,r*u.c] for u in x ] + else + r * x == + zero? r => 0 +-- one? r => x + (r = 1) => x + --map(r*#1,x) + [[u.k,a] for u in x | (a:=r*u.c) ^= 0$R] + if R has EntireRing then + x * r == + zero? r => 0 +-- one? r => x + (r = 1) => x + --map(r*#1,x) + [[u.k,u.c*r] for u in x ] + else + x * r == + zero? r => 0 +-- one? r => x + (r = 1) => x + --map(r*#1,x) + [[u.k,a] for u in x | (a:=u.c*r) ^= 0$R] + + coerce(x) : OutputForm == + null x => (0$R) :: OutputForm + le : List OutputForm := nil + for rec in reverse x repeat + rec.c = 1 => le := cons(rec.k :: OutputForm, le) + le := cons(rec.c :: OutputForm * rec.k :: OutputForm, le) + reduce("+",le) + +@ +\section{domain PR PolynomialRing} +<<domain PR PolynomialRing>>= +)abbrev domain PR PolynomialRing +++ Author: Dave Barton, James Davenport, Barry Trager +++ Date Created: +++ Date Last Updated: 14.08.2000. Improved exponentiation [MMM/TTT] +++ Basic Functions: Ring, degree, coefficient, monomial, reductum +++ Related Constructors: +++ Also See: +++ AMS Classifications: +++ Keywords: +++ References: +++ Description: +++ This domain represents generalized polynomials with coefficients +++ (from a not necessarily commutative ring), and terms +++ indexed by their exponents (from an arbitrary ordered abelian monoid). +++ This type is used, for example, +++ by the \spadtype{DistributedMultivariatePolynomial} domain where +++ the exponent domain is a direct product of non negative integers. + +PolynomialRing(R:Ring,E:OrderedAbelianMonoid): T == C + where + T == FiniteAbelianMonoidRing(R,E) with + --assertions + if R has IntegralDomain and E has CancellationAbelianMonoid then + fmecg: (%,E,R,%) -> % + ++ fmecg(p1,e,r,p2) finds X : p1 - r * X**e * p2 + if R has canonicalUnitNormal then canonicalUnitNormal + ++ canonicalUnitNormal guarantees that the function + ++ unitCanonical returns the same representative for all + ++ associates of any particular element. + + C == FreeModule(R,E) add + --representations + Term:= Record(k:E,c:R) + Rep:= List Term + + + --declarations + x,y,p,p1,p2: % + n: Integer + nn: NonNegativeInteger + np: PositiveInteger + e: E + r: R + --local operations + 1 == [[0$E,1$R]] + characteristic == characteristic$R + numberOfMonomials x == (# x)$Rep + degree p == if null p then 0 else p.first.k + minimumDegree p == if null p then 0 else (last p).k + leadingCoefficient p == if null p then 0$R else p.first.c + leadingMonomial p == if null p then 0 else [p.first] + reductum p == if null p then p else p.rest + retractIfCan(p:%):Union(R,"failed") == + null p => 0$R + not null p.rest => "failed" + zero?(p.first.k) => p.first.c + "failed" + coefficient(p,e) == + for tm in p repeat + tm.k=e => return tm.c + tm.k < e => return 0$R + 0$R + recip(p) == + null p => "failed" + p.first.k > 0$E => "failed" + (u:=recip(p.first.c)) case "failed" => "failed" + (u::R)::% + + coerce(r) == if zero? r then 0$% else [[0$E,r]] + coerce(n) == (n::R)::% + + ground?(p): Boolean == empty? p or (empty? rest p and zero? degree p) + + qsetrest!: (Rep, Rep) -> Rep + qsetrest!(l: Rep, e: Rep): Rep == RPLACD(l, e)$Lisp + + times!: (R, %) -> % + times: (R, E, %) -> % + + entireRing? := R has EntireRing + + times!(r: R, x: %): % == + res, endcell, newend, xx: Rep + if entireRing? then + for tx in x repeat tx.c := r*tx.c + else + xx := x + res := empty() + while not empty? xx repeat + tx := first xx + tx.c := r * tx.c + if zero? tx.c then + xx := rest xx + else + newend := xx + xx := rest xx + if empty? res then + res := newend + endcell := res + else + qsetrest!(endcell, newend) + endcell := newend + res; + + --- term * polynomial + termTimes: (R, E, Term) -> Term + termTimes(r: R, e: E, tx:Term): Term == [e+tx.k, r*tx.c] + times(tco: R, tex: E, rx: %): % == + if entireRing? then + map(termTimes(tco, tex, #1), rx::Rep) + else + [[tex + tx.k, r] for tx in rx::Rep | not zero? (r := tco * tx.c)] + + + + -- local addm! + addm!: (Rep, R, E, Rep) -> Rep + -- p1 + coef*x^E * p2 + -- `spare' (commented out) is for storage efficiency (not so good for + -- performance though. + addm!(p1:Rep, coef:R, exp: E, p2:Rep): Rep == + --local res, newend, last: Rep + res, newcell, endcell: Rep + spare: List Rep + res := empty() + endcell := empty() + --spare := empty() + while not empty? p1 and not empty? p2 repeat + tx := first p1 + ty := first p2 + exy := exp + ty.k + newcell := empty(); + if tx.k = exy then + newcoef := tx.c + coef * ty.c + if not zero? newcoef then + tx.c := newcoef + newcell := p1 + --else + -- spare := cons(p1, spare) + p1 := rest p1 + p2 := rest p2 + else if tx.k > exy then + newcell := p1 + p1 := rest p1 + else + newcoef := coef * ty.c + if not entireRing? and zero? newcoef then + newcell := empty() + --else if empty? spare then + -- ttt := [exy, newcoef] + -- newcell := cons(ttt, empty()) + --else + -- newcell := first spare + -- spare := rest spare + -- ttt := first newcell + -- ttt.k := exy + -- ttt.c := newcoef + else + ttt := [exy, newcoef] + newcell := cons(ttt, empty()) + p2 := rest p2 + if not empty? newcell then + if empty? res then + res := newcell + endcell := res + else + qsetrest!(endcell, newcell) + endcell := newcell + if not empty? p1 then -- then end is const * p1 + newcell := p1 + else -- then end is (coef, exp) * p2 + newcell := times(coef, exp, p2) + empty? res => newcell + qsetrest!(endcell, newcell) + res + pomopo! (p1, r, e, p2) == addm!(p1, r, e, p2) + p1 * p2 == + xx := p1::Rep + empty? xx => p1 + yy := p2::Rep + empty? yy => p2 + zero? first(xx).k => first(xx).c * p2 + zero? first(yy).k => p1 * first(yy).c + --if #xx > #yy then + -- (xx, yy) := (yy, xx) + -- (p1, p2) := (p2, p1) + xx := reverse xx + res : Rep := empty() + for tx in xx repeat res:=addm!(res,tx.c,tx.k,yy) + res + +-- if R has EntireRing then +-- p1 * p2 == +-- null p1 => 0 +-- null p2 => 0 +-- zero?(p1.first.k) => p1.first.c * p2 +-- one? p2 => p1 +-- +/[[[t1.k+t2.k,t1.c*t2.c]$Term for t2 in p2] +-- for t1 in reverse(p1)] +-- -- This 'reverse' is an efficiency improvement: +-- -- reduces both time and space [Abbott/Bradford/Davenport] +-- else +-- p1 * p2 == +-- null p1 => 0 +-- null p2 => 0 +-- zero?(p1.first.k) => p1.first.c * p2 +-- one? p2 => p1 +-- +/[[[t1.k+t2.k,r]$Term for t2 in p2 | (r:=t1.c*t2.c) ^= 0] +-- for t1 in reverse(p1)] +-- -- This 'reverse' is an efficiency improvement: +-- -- reduces both time and space [Abbott/Bradford/Davenport] + if R has CommutativeRing then + p ** np == p ** (np pretend NonNegativeInteger) + p ^ np == p ** (np pretend NonNegativeInteger) + p ^ nn == p ** nn + + + p ** nn == + null p => 0 + zero? nn => 1 +-- one? nn => p + (nn = 1) => p + empty? p.rest => + zero?(cc:=p.first.c ** nn) => 0 + [[nn * p.first.k, cc]] + binomThmExpt([p.first], p.rest, nn) + + if R has Field then + unitNormal(p) == + null p or (lcf:R:=p.first.c) = 1 => [1,p,1] + a := inv lcf + [lcf::%, [[p.first.k,1],:(a * p.rest)], a::%] + unitCanonical(p) == + null p or (lcf:R:=p.first.c) = 1 => p + a := inv lcf + [[p.first.k,1],:(a * p.rest)] + else if R has IntegralDomain then + unitNormal(p) == + null p or p.first.c = 1 => [1,p,1] + (u,cf,a):=unitNormal(p.first.c) + [u::%, [[p.first.k,cf],:(a * p.rest)], a::%] + unitCanonical(p) == + null p or p.first.c = 1 => p + (u,cf,a):=unitNormal(p.first.c) + [[p.first.k,cf],:(a * p.rest)] + if R has IntegralDomain then + associates?(p1,p2) == + null p1 => null p2 + null p2 => false + p1.first.k = p2.first.k and + associates?(p1.first.c,p2.first.c) and + ((p2.first.c exquo p1.first.c)::R * p1.rest = p2.rest) + p exquo r == + [(if (a:= tm.c exquo r) case "failed" + then return "failed" else [tm.k,a]) + for tm in p] :: Union(%,"failed") + if E has CancellationAbelianMonoid then + fmecg(p1:%,e:E,r:R,p2:%):% == -- p1 - r * X**e * p2 + rout:%:= [] + r:= - r + for tm in p2 repeat + e2:= e + tm.k + c2:= r * tm.c + c2 = 0 => "next term" + while not null p1 and p1.first.k > e2 repeat + (rout:=[p1.first,:rout]; p1:=p1.rest) --use PUSH and POP? + null p1 or p1.first.k < e2 => rout:=[[e2,c2],:rout] + if (u:=p1.first.c + c2) ^= 0 then rout:=[[e2, u],:rout] + p1:=p1.rest + NRECONC(rout,p1)$Lisp + if R has approximate then + p1 exquo p2 == + null p2 => error "Division by 0" + p2 = 1 => p1 + p1=p2 => 1 + --(p1.lastElt.c exquo p2.lastElt.c) case "failed" => "failed" + rout:= []@List(Term) + while not null p1 repeat + (a:= p1.first.c exquo p2.first.c) + a case "failed" => return "failed" + ee:= subtractIfCan(p1.first.k, p2.first.k) + ee case "failed" => return "failed" + p1:= fmecg(p1.rest, ee, a, p2.rest) + rout:= [[ee,a], :rout] + null p1 => reverse(rout)::% -- nreverse? + "failed" + else -- R not approximate + p1 exquo p2 == + null p2 => error "Division by 0" + p2 = 1 => p1 + --(p1.lastElt.c exquo p2.lastElt.c) case "failed" => "failed" + rout:= []@List(Term) + while not null p1 repeat + (a:= p1.first.c exquo p2.first.c) + a case "failed" => return "failed" + ee:= subtractIfCan(p1.first.k, p2.first.k) + ee case "failed" => return "failed" + p1:= fmecg(p1.rest, ee, a, p2.rest) + rout:= [[ee,a], :rout] + null p1 => reverse(rout)::% -- nreverse? + "failed" + if R has Field then + x/r == inv(r)*x + +@ +\section{domain SUP SparseUnivariatePolynomial} +<<domain SUP SparseUnivariatePolynomial>>= +)abbrev domain SUP SparseUnivariatePolynomial +++ Author: Dave Barton, Barry Trager +++ Date Created: +++ Date Last Updated: +++ Basic Functions: Ring, monomial, coefficient, reductum, differentiate, +++ elt, map, resultant, discriminant +++ Related Constructors: UnivariatePolynomial, Polynomial +++ Also See: +++ AMS Classifications: +++ Keywords: +++ References: +++ Description: +++ This domain represents univariate polynomials over arbitrary +++ (not necessarily commutative) coefficient rings. The variable is +++ unspecified so that the variable displays as \spad{?} on output. +++ If it is necessary to specify the variable name, use type \spadtype{UnivariatePolynomial}. +++ The representation is sparse +++ in the sense that only non-zero terms are represented. +++ Note: if the coefficient ring is a field, this domain forms a euclidean domain. + +SparseUnivariatePolynomial(R:Ring): UnivariatePolynomialCategory(R) with + outputForm : (%,OutputForm) -> OutputForm + ++ outputForm(p,var) converts the SparseUnivariatePolynomial p to + ++ an output form (see \spadtype{OutputForm}) printed as a polynomial in the + ++ output form variable. + fmecg: (%,NonNegativeInteger,R,%) -> % + ++ fmecg(p1,e,r,p2) finds X : p1 - r * X**e * p2 + == PolynomialRing(R,NonNegativeInteger) + add + --representations + Term := Record(k:NonNegativeInteger,c:R) + Rep := List Term + p:% + n:NonNegativeInteger + np: PositiveInteger + FP ==> SparseUnivariatePolynomial % + pp,qq: FP + lpp:List FP + + -- for karatsuba + kBound: NonNegativeInteger := 63 + upmp := UnivariatePolynomialMultiplicationPackage(R,%) + + + if R has FieldOfPrimeCharacteristic then + p ** np == p ** (np pretend NonNegativeInteger) + p ^ np == p ** (np pretend NonNegativeInteger) + p ^ n == p ** n + p ** n == + null p => 0 + zero? n => 1 +-- one? n => p + (n = 1) => p + empty? p.rest => + zero?(cc:=p.first.c ** n) => 0 + [[n * p.first.k, cc]] + -- not worth doing special trick if characteristic is too small + if characteristic()$R < 3 then return expt(p,n pretend PositiveInteger)$RepeatedSquaring(%) + y:%:=1 + -- break up exponent in qn * characteristic + rn + -- exponentiating by the characteristic is fast + rec := divide(n, characteristic()$R) + qn:= rec.quotient + rn:= rec.remainder + repeat + if rn = 1 then y := y * p + if rn > 1 then y:= y * binomThmExpt([p.first], p.rest, rn) + zero? qn => return y + -- raise to the characteristic power + p:= [[t.k * characteristic()$R , primeFrobenius(t.c)$R ]$Term for t in p] + rec := divide(qn, characteristic()$R) + qn:= rec.quotient + rn:= rec.remainder + y + + + + zero?(p): Boolean == empty?(p) +-- one?(p):Boolean == not empty? p and (empty? rest p and zero? first(p).k and one? first(p).c) + one?(p):Boolean == not empty? p and (empty? rest p and zero? first(p).k and (first(p).c = 1)) + ground?(p): Boolean == empty? p or (empty? rest p and zero? first(p).k) + multiplyExponents(p,n) == [ [u.k*n,u.c] for u in p] + divideExponents(p,n) == + null p => p + m:= (p.first.k :: Integer exquo n::Integer) + m case "failed" => "failed" + u:= divideExponents(p.rest,n) + u case "failed" => "failed" + [[m::Integer::NonNegativeInteger,p.first.c],:u] + karatsubaDivide(p, n) == + zero? n => [p, 0] + lowp: Rep := p + highp: Rep := [] + repeat + if empty? lowp then break + t := first lowp + if t.k < n then break + lowp := rest lowp + highp := cons([subtractIfCan(t.k,n)::NonNegativeInteger,t.c]$Term,highp) + [ reverse highp, lowp] + shiftRight(p, n) == + [[subtractIfCan(t.k,n)::NonNegativeInteger,t.c]$Term for t in p] + shiftLeft(p, n) == + [[t.k + n,t.c]$Term for t in p] + pomopo!(p1,r,e,p2) == + rout:%:= [] + for tm in p2 repeat + e2:= e + tm.k + c2:= r * tm.c + c2 = 0 => "next term" + while not null p1 and p1.first.k > e2 repeat + (rout:=[p1.first,:rout]; p1:=p1.rest) --use PUSH and POP? + null p1 or p1.first.k < e2 => rout:=[[e2,c2],:rout] + if (u:=p1.first.c + c2) ^= 0 then rout:=[[e2, u],:rout] + p1:=p1.rest + NRECONC(rout,p1)$Lisp + +-- implementation using karatsuba algorithm conditionally +-- +-- p1 * p2 == +-- xx := p1::Rep +-- empty? xx => p1 +-- yy := p2::Rep +-- empty? yy => p2 +-- zero? first(xx).k => first(xx).c * p2 +-- zero? first(yy).k => p1 * first(yy).c +-- (first(xx).k > kBound) and (first(yy).k > kBound) and (#xx > kBound) and (#yy > kBound) => +-- karatsubaOnce(p1,p2)$upmp +-- xx := reverse xx +-- res : Rep := empty() +-- for tx in xx repeat res:= rep pomopo!( res,tx.c,tx.k,p2) +-- res + + + univariate(p:%) == p pretend SparseUnivariatePolynomial(R) + multivariate(sup:SparseUnivariatePolynomial(R),v:SingletonAsOrderedSet) == + sup pretend % + univariate(p:%,v:SingletonAsOrderedSet) == + zero? p => 0 + monomial(leadingCoefficient(p)::%,degree p) + + univariate(reductum p,v) + multivariate(supp:SparseUnivariatePolynomial(%),v:SingletonAsOrderedSet) == + zero? supp => 0 + lc:=leadingCoefficient supp + degree lc > 0 => error "bad form polynomial" + monomial(leadingCoefficient lc,degree supp) + + multivariate(reductum supp,v) + if R has FiniteFieldCategory and R has PolynomialFactorizationExplicit then + RXY ==> SparseUnivariatePolynomial SparseUnivariatePolynomial R + squareFreePolynomial pp == + squareFree(pp)$UnivariatePolynomialSquareFree(%,FP) + factorPolynomial pp == + (generalTwoFactor(pp pretend RXY)$TwoFactorize(R)) + pretend Factored SparseUnivariatePolynomial % + factorSquareFreePolynomial pp == + (generalTwoFactor(pp pretend RXY)$TwoFactorize(R)) + pretend Factored SparseUnivariatePolynomial % + gcdPolynomial(pp,qq) == gcd(pp,qq)$FP + factor p == factor(p)$DistinctDegreeFactorize(R,%) + solveLinearPolynomialEquation(lpp,pp) == + solveLinearPolynomialEquation(lpp, pp)$FiniteFieldSolveLinearPolynomialEquation(R,%,FP) + else if R has PolynomialFactorizationExplicit then + import PolynomialFactorizationByRecursionUnivariate(R,%) + solveLinearPolynomialEquation(lpp,pp)== + solveLinearPolynomialEquationByRecursion(lpp,pp) + factorPolynomial(pp) == + factorByRecursion(pp) + factorSquareFreePolynomial(pp) == + factorSquareFreeByRecursion(pp) + + if R has IntegralDomain then + if R has approximate then + p1 exquo p2 == + null p2 => error "Division by 0" + p2 = 1 => p1 + p1=p2 => 1 + --(p1.lastElt.c exquo p2.lastElt.c) case "failed" => "failed" + rout:= []@List(Term) + while not null p1 repeat + (a:= p1.first.c exquo p2.first.c) + a case "failed" => return "failed" + ee:= subtractIfCan(p1.first.k, p2.first.k) + ee case "failed" => return "failed" + p1:= fmecg(p1.rest, ee, a, p2.rest) + rout:= [[ee,a], :rout] + null p1 => reverse(rout)::% -- nreverse? + "failed" + else -- R not approximate + p1 exquo p2 == + null p2 => error "Division by 0" + p2 = 1 => p1 + --(p1.lastElt.c exquo p2.lastElt.c) case "failed" => "failed" + rout:= []@List(Term) + while not null p1 repeat + (a:= p1.first.c exquo p2.first.c) + a case "failed" => return "failed" + ee:= subtractIfCan(p1.first.k, p2.first.k) + ee case "failed" => return "failed" + p1:= fmecg(p1.rest, ee, a, p2.rest) + rout:= [[ee,a], :rout] + null p1 => reverse(rout)::% -- nreverse? + "failed" + fmecg(p1,e,r,p2) == -- p1 - r * X**e * p2 + rout:%:= [] + r:= - r + for tm in p2 repeat + e2:= e + tm.k + c2:= r * tm.c + c2 = 0 => "next term" + while not null p1 and p1.first.k > e2 repeat + (rout:=[p1.first,:rout]; p1:=p1.rest) --use PUSH and POP? + null p1 or p1.first.k < e2 => rout:=[[e2,c2],:rout] + if (u:=p1.first.c + c2) ^= 0 then rout:=[[e2, u],:rout] + p1:=p1.rest + NRECONC(rout,p1)$Lisp + pseudoRemainder(p1,p2) == + null p2 => error "PseudoDivision by Zero" + null p1 => 0 + co:=p2.first.c; + e:=p2.first.k; + p2:=p2.rest; + e1:=max(p1.first.k:Integer-e+1,0):NonNegativeInteger + while not null p1 repeat + if (u:=subtractIfCan(p1.first.k,e)) case "failed" then leave + p1:=fmecg(co * p1.rest, u, p1.first.c, p2) + e1:= (e1 - 1):NonNegativeInteger + e1 = 0 => p1 + co ** e1 * p1 + toutput(t1:Term,v:OutputForm):OutputForm == + t1.k = 0 => t1.c :: OutputForm + if t1.k = 1 + then mon:= v + else mon := v ** t1.k::OutputForm + t1.c = 1 => mon + t1.c = -1 and + ((t1.c :: OutputForm) = (-1$Integer)::OutputForm)@Boolean => - mon + t1.c::OutputForm * mon + outputForm(p:%,v:OutputForm) == + l: List(OutputForm) + l:=[toutput(t,v) for t in p] + null l => (0$Integer)::OutputForm -- else FreeModule 0 problems + reduce("+",l) + + coerce(p:%):OutputForm == outputForm(p, "?"::OutputForm) + elt(p:%,val:R) == + null p => 0$R + co:=p.first.c + n:=p.first.k + for tm in p.rest repeat + co:= co * val ** (n - (n:=tm.k)):NonNegativeInteger + tm.c + n = 0 => co + co * val ** n + elt(p:%,val:%) == + null p => 0$% + coef:% := p.first.c :: % + n:=p.first.k + for tm in p.rest repeat + coef:= coef * val ** (n-(n:=tm.k)):NonNegativeInteger+(tm.c::%) + n = 0 => coef + coef * val ** n + + monicDivide(p1:%,p2:%) == + null p2 => error "monicDivide: division by 0" + leadingCoefficient p2 ^= 1 => error "Divisor Not Monic" + p2 = 1 => [p1,0] + null p1 => [0,0] + degree p1 < (n:=degree p2) => [0,p1] + rout:Rep := [] + p2 := p2.rest + while not null p1 repeat + (u:=subtractIfCan(p1.first.k, n)) case "failed" => leave + rout:=[[u, p1.first.c], :rout] + p1:=fmecg(p1.rest, rout.first.k, rout.first.c, p2) + [reverse_!(rout),p1] + + if R has IntegralDomain then + discriminant(p) == discriminant(p)$PseudoRemainderSequence(R,%) +-- discriminant(p) == +-- null p or zero?(p.first.k) => error "cannot take discriminant of constants" +-- dp:=differentiate p +-- corr:= p.first.c ** ((degree p - 1 - degree dp)::NonNegativeInteger) +-- (-1)**((p.first.k*(p.first.k-1)) quo 2):NonNegativeInteger +-- * (corr * resultant(p,dp) exquo p.first.c)::R + + subResultantGcd(p1,p2) == subResultantGcd(p1,p2)$PseudoRemainderSequence(R,%) +-- subResultantGcd(p1,p2) == --args # 0, non-coef, prim, ans not prim +-- --see algorithm 1 (p. 4) of Brown's latest (unpublished) paper +-- if p1.first.k < p2.first.k then (p1,p2):=(p2,p1) +-- p:=pseudoRemainder(p1,p2) +-- co:=1$R; +-- e:= (p1.first.k - p2.first.k):NonNegativeInteger +-- while not null p and p.first.k ^= 0 repeat +-- p1:=p2; p2:=p; p:=pseudoRemainder(p1,p2) +-- null p or p.first.k = 0 => "enuf" +-- co:=(p1.first.c ** e exquo co ** max(0, (e-1))::NonNegativeInteger)::R +-- e:= (p1.first.k - p2.first.k):NonNegativeInteger; c1:=co**e +-- p:=[[tm.k,((tm.c exquo p1.first.c)::R exquo c1)::R] for tm in p] +-- if null p then p2 else 1$% + + resultant(p1,p2) == resultant(p1,p2)$PseudoRemainderSequence(R,%) +-- resultant(p1,p2) == --SubResultant PRS Algorithm +-- null p1 or null p2 => 0$R +-- 0 = degree(p1) => ((first p1).c)**degree(p2) +-- 0 = degree(p2) => ((first p2).c)**degree(p1) +-- if p1.first.k < p2.first.k then +-- (if odd?(p1.first.k) then p1:=-p1; (p1,p2):=(p2,p1)) +-- p:=pseudoRemainder(p1,p2) +-- co:=1$R; e:=(p1.first.k-p2.first.k):NonNegativeInteger +-- while not null p repeat +-- if not odd?(e) then p:=-p +-- p1:=p2; p2:=p; p:=pseudoRemainder(p1,p2) +-- co:=(p1.first.c ** e exquo co ** max(e:Integer-1,0):NonNegativeInteger)::R +-- e:= (p1.first.k - p2.first.k):NonNegativeInteger; c1:=co**e +-- p:=(p exquo ((leadingCoefficient p1) * c1))::% +-- degree p2 > 0 => 0$R +-- (p2.first.c**e exquo co**((e-1)::NonNegativeInteger))::R + if R has GcdDomain then + content(p) == if null p then 0$R else "gcd"/[tm.c for tm in p] + --make CONTENT more efficient? + + primitivePart(p) == + null p => p + ct :=content(p) + unitCanonical((p exquo ct)::%) + -- exquo present since % is now an IntegralDomain + + gcd(p1,p2) == + gcdPolynomial(p1 pretend SparseUnivariatePolynomial R, + p2 pretend SparseUnivariatePolynomial R) pretend % + + if R has Field then + divide( p1, p2) == + zero? p2 => error "Division by 0" +-- one? p2 => [p1,0] + (p2 = 1) => [p1,0] + ct:=inv(p2.first.c) + n:=p2.first.k + p2:=p2.rest + rout:=empty()$List(Term) + while p1 ^= 0 repeat + (u:=subtractIfCan(p1.first.k, n)) case "failed" => leave + rout:=[[u, ct * p1.first.c], :rout] + p1:=fmecg(p1.rest, rout.first.k, rout.first.c, p2) + [reverse_!(rout),p1] + + p / co == inv(co) * p + +@ +\section{package SUP2 SparseUnivariatePolynomialFunctions2} +<<package SUP2 SparseUnivariatePolynomialFunctions2>>= +)abbrev package SUP2 SparseUnivariatePolynomialFunctions2 +++ Author: +++ Date Created: +++ Date Last Updated: +++ Basic Functions: +++ Related Constructors: +++ Also See: +++ AMS Classifications: +++ Keywords: +++ References: +++ Description: +++ This package lifts a mapping from coefficient rings R to S to +++ a mapping from sparse univariate polynomial over R to +++ a sparse univariate polynomial over S. +++ Note that the mapping is assumed +++ to send zero to zero, since it will only be applied to the non-zero +++ coefficients of the polynomial. + +SparseUnivariatePolynomialFunctions2(R:Ring, S:Ring): with + map:(R->S,SparseUnivariatePolynomial R) -> SparseUnivariatePolynomial S + ++ map(func, poly) creates a new polynomial by applying func to + ++ every non-zero coefficient of the polynomial poly. + == add + map(f, p) == map(f, p)$UnivariatePolynomialCategoryFunctions2(R, + SparseUnivariatePolynomial R, S, SparseUnivariatePolynomial S) + +@ +\section{domain UP UnivariatePolynomial} +<<domain UP UnivariatePolynomial>>= +)abbrev domain UP UnivariatePolynomial +++ Author: +++ Date Created: +++ Date Last Updated: +++ Basic Functions: Ring, monomial, coefficient, reductum, differentiate, +++ elt, map, resultant, discriminant +++ Related Constructors: SparseUnivariatePolynomial, MultivariatePolynomial +++ Also See: +++ AMS Classifications: +++ Keywords: +++ References: +++ Description: +++ This domain represents univariate polynomials in some symbol +++ over arbitrary (not necessarily commutative) coefficient rings. +++ The representation is sparse +++ in the sense that only non-zero terms are represented. +++ Note: if the coefficient ring is a field, then this domain forms a euclidean domain. + +UnivariatePolynomial(x:Symbol, R:Ring): + UnivariatePolynomialCategory(R) with + coerce: Variable(x) -> % + ++ coerce(x) converts the variable x to a univariate polynomial. + fmecg: (%,NonNegativeInteger,R,%) -> % + ++ fmecg(p1,e,r,p2) finds X : p1 - r * X**e * p2 + == SparseUnivariatePolynomial(R) add + Rep:=SparseUnivariatePolynomial(R) + coerce(p:%):OutputForm == outputForm(p, outputForm x) + coerce(v:Variable(x)):% == monomial(1, 1) + +@ +\section{package UP2 UnivariatePolynomialFunctions2} +<<package UP2 UnivariatePolynomialFunctions2>>= +)abbrev package UP2 UnivariatePolynomialFunctions2 +++ Author: +++ Date Created: +++ Date Last Updated: +++ Basic Functions: +++ Related Constructors: +++ Also See: +++ AMS Classifications: +++ Keywords: +++ References: +++ Description: +++ This package lifts a mapping from coefficient rings R to S to +++ a mapping from \spadtype{UnivariatePolynomial}(x,R) to +++ \spadtype{UnivariatePolynomial}(y,S). Note that the mapping is assumed +++ to send zero to zero, since it will only be applied to the non-zero +++ coefficients of the polynomial. + +UnivariatePolynomialFunctions2(x:Symbol, R:Ring, y:Symbol, S:Ring): with + map: (R -> S, UnivariatePolynomial(x,R)) -> UnivariatePolynomial(y,S) + ++ map(func, poly) creates a new polynomial by applying func to + ++ every non-zero coefficient of the polynomial poly. + == add + map(f, p) == map(f, p)$UnivariatePolynomialCategoryFunctions2(R, + UnivariatePolynomial(x, R), S, UnivariatePolynomial(y, S)) + +@ +\section{package POLY2UP PolynomialToUnivariatePolynomial} +<<package POLY2UP PolynomialToUnivariatePolynomial>>= +)abbrev package POLY2UP PolynomialToUnivariatePolynomial +++ Author: +++ Date Created: +++ Date Last Updated: +++ Basic Functions: +++ Related Constructors: +++ Also See: +++ AMS Classifications: +++ Keywords: +++ References: +++ Description: +++ This package is primarily to help the interpreter do coercions. +++ It allows you to view a polynomial as a +++ univariate polynomial in one of its variables with +++ coefficients which are again a polynomial in all the +++ other variables. + +PolynomialToUnivariatePolynomial(x:Symbol, R:Ring): with + univariate: (Polynomial R, Variable x) -> + UnivariatePolynomial(x, Polynomial R) + ++ univariate(p, x) converts the polynomial p to a one of type + ++ \spad{UnivariatePolynomial(x,Polynomial(R))}, ie. as a member of \spad{R[...][x]}. + == add + univariate(p, y) == + q:SparseUnivariatePolynomial(Polynomial R) := univariate(p, x) + map(#1, q)$UnivariatePolynomialCategoryFunctions2(Polynomial R, + SparseUnivariatePolynomial Polynomial R, Polynomial R, + UnivariatePolynomial(x, Polynomial R)) + +@ +\section{package UPSQFREE UnivariatePolynomialSquareFree} +<<package UPSQFREE UnivariatePolynomialSquareFree>>= +)abbrev package UPSQFREE UnivariatePolynomialSquareFree +++ Author: Dave Barton, Barry Trager +++ Date Created: +++ Date Last Updated: +++ Basic Functions: squareFree, squareFreePart +++ Related Constructors: +++ Also See: +++ AMS Classifications: +++ Keywords: +++ References: +++ Description: +++ This package provides for square-free decomposition of +++ univariate polynomials over arbitrary rings, i.e. +++ a partial factorization such that each factor is a product +++ of irreducibles with multiplicity one and the factors are +++ pairwise relatively prime. If the ring +++ has characteristic zero, the result is guaranteed to satisfy +++ this condition. If the ring is an infinite ring of +++ finite characteristic, then it may not be possible to decide when +++ polynomials contain factors which are pth powers. In this +++ case, the flag associated with that polynomial is set to "nil" +++ (meaning that that polynomials are not guaranteed to be square-free). + +UnivariatePolynomialSquareFree(RC:IntegralDomain,P):C == T + where + fUnion ==> Union("nil", "sqfr", "irred", "prime") + FF ==> Record(flg:fUnion, fctr:P, xpnt:Integer) + P:Join(UnivariatePolynomialCategory(RC),IntegralDomain) with + gcd: (%,%) -> % + ++ gcd(p,q) computes the greatest-common-divisor of p and q. + + C == with + squareFree: P -> Factored(P) + ++ squareFree(p) computes the square-free factorization of the + ++ univariate polynomial p. Each factor has no repeated roots, and the + ++ factors are pairwise relatively prime. + squareFreePart: P -> P + ++ squareFreePart(p) returns a polynomial which has the same + ++ irreducible factors as the univariate polynomial p, but each + ++ factor has multiplicity one. + BumInSepFFE: FF -> FF + ++ BumInSepFFE(f) is a local function, exported only because + ++ it has multiple conditional definitions. + + T == add + + if RC has CharacteristicZero then + squareFreePart(p:P) == (p exquo gcd(p, differentiate p))::P + else + squareFreePart(p:P) == + unit(s := squareFree(p)$%) * */[f.factor for f in factors s] + + if RC has FiniteFieldCategory then + BumInSepFFE(ffe:FF) == + ["sqfr", map(charthRoot,ffe.fctr), characteristic$P*ffe.xpnt] + else if RC has CharacteristicNonZero then + BumInSepFFE(ffe:FF) == + np := multiplyExponents(ffe.fctr,characteristic$P:NonNegativeInteger) + (nthrp := charthRoot(np)) case "failed" => + ["nil", np, ffe.xpnt] + ["sqfr", nthrp, characteristic$P*ffe.xpnt] + + else + BumInSepFFE(ffe:FF) == + ["nil", + multiplyExponents(ffe.fctr,characteristic$P:NonNegativeInteger), + ffe.xpnt] + + + if RC has CharacteristicZero then + squareFree(p:P) == --Yun's algorithm - see SYMSAC '76, p.27 + --Note ci primitive is, so GCD's don't need to %do contents. + --Change gcd to return cofctrs also? + ci:=p; di:=differentiate(p); pi:=gcd(ci,di) + degree(pi)=0 => + (u,c,a):=unitNormal(p) + makeFR(u,[["sqfr",c,1]]) + i:NonNegativeInteger:=0; lffe:List FF:=[] + lcp := leadingCoefficient p + while degree(ci)^=0 repeat + ci:=(ci exquo pi)::P + di:=(di exquo pi)::P - differentiate(ci) + pi:=gcd(ci,di) + i:=i+1 + degree(pi) > 0 => + lcp:=(lcp exquo (leadingCoefficient(pi)**i))::RC + lffe:=[["sqfr",pi,i],:lffe] + makeFR(lcp::P,lffe) + + else + squareFree(p:P) == --Musser's algorithm - see SYMSAC '76, p.27 + --p MUST BE PRIMITIVE, Any characteristic. + --Note ci primitive, so GCD's don't need to %do contents. + --Change gcd to return cofctrs also? + ci := gcd(p,differentiate(p)) + degree(ci)=0 => + (u,c,a):=unitNormal(p) + makeFR(u,[["sqfr",c,1]]) + di := (p exquo ci)::P + i:NonNegativeInteger:=0; lffe:List FF:=[] + dunit : P := 1 + while degree(di)^=0 repeat + diprev := di + di := gcd(ci,di) + ci:=(ci exquo di)::P + i:=i+1 + degree(diprev) = degree(di) => + lc := (leadingCoefficient(diprev) exquo leadingCoefficient(di))::RC + dunit := lc**i * dunit + pi:=(diprev exquo di)::P + lffe:=[["sqfr",pi,i],:lffe] + dunit := dunit * di ** (i+1) + degree(ci)=0 => makeFR(dunit*ci,lffe) + redSqfr:=squareFree(divideExponents(ci,characteristic$P)::P) + lsnil:= [BumInSepFFE(ffe) for ffe in factorList redSqfr] + lffe:=append(lsnil,lffe) + makeFR(dunit*(unit redSqfr),lffe) + +@ +\section{package PSQFR PolynomialSquareFree} +<<package PSQFR PolynomialSquareFree>>= +)abbrev package PSQFR PolynomialSquareFree +++ Author: +++ Date Created: +++ Date Last Updated: November 1993, (P.Gianni) +++ Basic Functions: +++ Related Constructors: +++ Also See: +++ AMS Classifications: +++ Keywords: +++ References: +++ Description: +++ This package computes square-free decomposition of multivariate +++ polynomials over a coefficient ring which is an arbitrary gcd domain. +++ The requirement on the coefficient domain guarantees that the \spadfun{content} can be +++ removed so that factors will be primitive as well as square-free. +++ Over an infinite ring of finite characteristic,it may not be possible to +++ guarantee that the factors are square-free. + +PolynomialSquareFree(VarSet:OrderedSet,E,RC:GcdDomain,P):C == T where + E:OrderedAbelianMonoidSup + P:PolynomialCategory(RC,E,VarSet) + + C == with + squareFree : P -> Factored P + ++ squareFree(p) returns the square-free factorization of the + ++ polynomial p. Each factor has no repeated roots, and the + ++ factors are pairwise relatively prime. + + T == add + SUP ==> SparseUnivariatePolynomial(P) + NNI ==> NonNegativeInteger + fUnion ==> Union("nil", "sqfr", "irred", "prime") + FF ==> Record(flg:fUnion, fctr:P, xpnt:Integer) + + finSqFr : (P,List VarSet) -> Factored P + pthPower : P -> Factored P + pPolRoot : P -> P + putPth : P -> P + + chrc:=characteristic$RC + + if RC has CharacteristicNonZero then + -- find the p-th root of a polynomial + pPolRoot(f:P) : P == + lvar:=variables f + empty? lvar => f + mv:=first lvar + uf:=univariate(f,mv) + uf:=divideExponents(uf,chrc)::SUP + uf:=map(pPolRoot,uf) + multivariate(uf,mv) + + -- substitute variables with their p-th power + putPth(f:P) : P == + lvar:=variables f + empty? lvar => f + mv:=first lvar + uf:=univariate(f,mv) + uf:=multiplyExponents(uf,chrc)::SUP + uf:=map(putPth,uf) + multivariate(uf,mv) + + -- the polynomial is a perfect power + pthPower(f:P) : Factored P == + proot : P := 0 + isSq : Boolean := false + if (g:=charthRoot f) case "failed" then proot:=pPolRoot(f) + else + proot := g :: P + isSq := true + psqfr:=finSqFr(proot,variables f) + isSq => + makeFR((unit psqfr)**chrc,[[u.flg,u.fctr, + (u.xpnt)*chrc] for u in factorList psqfr]) + makeFR((unit psqfr),[["nil",putPth u.fctr,u.xpnt] + for u in factorList psqfr]) + + -- compute the square free decomposition, finite characteristic case + finSqFr(f:P,lvar:List VarSet) : Factored P == + empty? lvar => pthPower(f) + mv:=first lvar + lvar:=lvar.rest + differentiate(f,mv)=0 => finSqFr(f,lvar) + uf:=univariate(f,mv) + cont := content uf + cont1:P:=1 + uf := (uf exquo cont)::SUP + squf := squareFree(uf)$UnivariatePolynomialSquareFree(P,SUP) + pfaclist:List FF :=[] + for u in factorList squf repeat + uexp:NNI:=(u.xpnt):NNI + u.flg = "sqfr" => -- the square free factor is OK + pfaclist:= cons([u.flg,multivariate(u.fctr,mv),uexp], + pfaclist) + --listfin1:= finSqFr(multivariate(u.fctr,mv),lvar) + listfin1:= squareFree multivariate(u.fctr,mv) + flistfin1:=[[uu.flg,uu.fctr,uu.xpnt*uexp] + for uu in factorList listfin1] + cont1:=cont1*((unit listfin1)**uexp) + pfaclist:=append(flistfin1,pfaclist) + cont:=cont*cont1 + cont ^= 1 => + sqp := squareFree cont + pfaclist:= append (factorList sqp,pfaclist) + makeFR(unit(sqp)*coefficient(unit squf,0),pfaclist) + makeFR(coefficient(unit squf,0),pfaclist) + + squareFree(p:P) == + mv:=mainVariable p + mv case "failed" => makeFR(p,[])$Factored(P) + characteristic$RC ^=0 => finSqFr(p,variables p) + up:=univariate(p,mv) + cont := content up + up := (up exquo cont)::SUP + squp := squareFree(up)$UnivariatePolynomialSquareFree(P,SUP) + pfaclist:List FF := + [[u.flg,multivariate(u.fctr,mv),u.xpnt] + for u in factorList squp] + cont ^= 1 => + sqp := squareFree cont + makeFR(unit(sqp)*coefficient(unit squp,0), + append(factorList sqp, pfaclist)) + makeFR(coefficient(unit squp,0),pfaclist) + +@ +\section{package UPMP UnivariatePolynomialMultiplicationPackage} +<<package UPMP UnivariatePolynomialMultiplicationPackage>>= +)abbrev package UPMP UnivariatePolynomialMultiplicationPackage +++ Author: Marc Moreno Maza +++ Date Created: 14.08.2000 +++ Description: +++ This package implements Karatsuba's trick for multiplying +++ (large) univariate polynomials. It could be improved with +++ a version doing the work on place and also with a special +++ case for squares. We've done this in Basicmath, but we +++ believe that this out of the scope of AXIOM. + +UnivariatePolynomialMultiplicationPackage(R: Ring, U: UnivariatePolynomialCategory(R)): C == T + where + HL ==> Record(quotient:U,remainder:U) + C == with + noKaratsuba: (U, U) -> U + ++ \spad{noKaratsuba(a,b)} returns \spad{a*b} without + ++ using Karatsuba's trick at all. + karatsubaOnce: (U, U) -> U + ++ \spad{karatsuba(a,b)} returns \spad{a*b} by applying + ++ Karatsuba's trick once. The other multiplications + ++ are performed by calling \spad{*} from \spad{U}. + karatsuba: (U, U, NonNegativeInteger, NonNegativeInteger) -> U; + ++ \spad{karatsuba(a,b,l,k)} returns \spad{a*b} by applying + ++ Karatsuba's trick provided that both \spad{a} and \spad{b} + ++ have at least \spad{l} terms and \spad{k > 0} holds + ++ and by calling \spad{noKaratsuba} otherwise. The other + ++ multiplications are performed by recursive calls with + ++ the same third argument and \spad{k-1} as fourth argument. + + T == add + noKaratsuba(a,b) == + zero? a => a + zero? b => b + zero?(degree(a)) => leadingCoefficient(a) * b + zero?(degree(b)) => a * leadingCoefficient(b) + lu: List(U) := reverse monomials(a) + res: U := 0; + for u in lu repeat + res := pomopo!(res, leadingCoefficient(u), degree(u), b) + res + karatsubaOnce(a:U,b:U): U == + da := minimumDegree(a) + db := minimumDegree(b) + if not zero? da then a := shiftRight(a,da) + if not zero? db then b := shiftRight(b,db) + d := da + db + n: NonNegativeInteger := min(degree(a),degree(b)) quo 2 + rec: HL := karatsubaDivide(a, n) + ha := rec.quotient + la := rec.remainder + rec := karatsubaDivide(b, n) + hb := rec.quotient + lb := rec.remainder + w: U := (ha - la) * (lb - hb) + u: U := la * lb + v: U := ha * hb + w := w + (u + v) + w := shiftLeft(w,n) + u + zero? d => shiftLeft(v,2*n) + w + shiftLeft(v,2*n + d) + shiftLeft(w,d) + karatsuba(a:U,b:U,l:NonNegativeInteger,k:NonNegativeInteger): U == + zero? k => noKaratsuba(a,b) + degree(a) < l => noKaratsuba(a,b) + degree(b) < l => noKaratsuba(a,b) + numberOfMonomials(a) < l => noKaratsuba(a,b) + numberOfMonomials(b) < l => noKaratsuba(a,b) + da := minimumDegree(a) + db := minimumDegree(b) + if not zero? da then a := shiftRight(a,da) + if not zero? db then b := shiftRight(b,db) + d := da + db + n: NonNegativeInteger := min(degree(a),degree(b)) quo 2 + k := subtractIfCan(k,1)::NonNegativeInteger + rec: HL := karatsubaDivide(a, n) + ha := rec.quotient + la := rec.remainder + rec := karatsubaDivide(b, n) + hb := rec.quotient + lb := rec.remainder + w: U := karatsuba(ha - la, lb - hb, l, k) + u: U := karatsuba(la, lb, l, k) + v: U := karatsuba(ha, hb, l, k) + w := w + (u + v) + w := shiftLeft(w,n) + u + zero? d => shiftLeft(v,2*n) + w + shiftLeft(v,2*n + d) + shiftLeft(w,d) + +@ +\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 FM FreeModule>> +<<domain PR PolynomialRing>> +<<package UPSQFREE UnivariatePolynomialSquareFree>> +<<package PSQFR PolynomialSquareFree>> +<<package UPMP UnivariatePolynomialMultiplicationPackage>> +<<domain SUP SparseUnivariatePolynomial>> +<<package SUP2 SparseUnivariatePolynomialFunctions2>> +<<domain UP UnivariatePolynomial>> +<<package UP2 UnivariatePolynomialFunctions2>> +<<package POLY2UP PolynomialToUnivariatePolynomial>> +@ +\eject +\begin{thebibliography}{99} +\bibitem{1} nothing +\end{thebibliography} +\end{document} |