\documentclass{article} \usepackage{open-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} <>= )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:OrderedType): Join(BiModule(R,R),IndexedDirectProductCategory(R,S)) with if R has CommutativeRing then Module(R) == IndexedDirectProductAbelianGroup(R,S) add if R has EntireRing then r:R * x:% == zero? r => 0$% one? r => x convert [term(index u, r * coefficient u) for u in terms x ] else r:R * x:% == zero? r => 0$% one? r => x convert [term(index u,a) for u in terms x | not zero?(a := r * coefficient u)] if R has EntireRing then x:% * r:R == zero? r => 0$% one? r => x convert [term(index u, coefficient(u) * r) for u in terms x ] else x:% * r:R == zero? r => 0$% one? r => x convert [term(index u,a) for u in terms x | not zero?(a := coefficient(u) * r)] if S has CoercibleTo OutputForm then coerce(x) : OutputForm == zero? x => (0$R) :: OutputForm le : List OutputForm := nil for rec in reverse terms x repeat one? coefficient rec => le := cons(index rec :: OutputForm, le) le := cons(coefficient rec :: OutputForm * index rec :: OutputForm, le) reduce("+",le) @ \section{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" positive? p.first.k => "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) import %tail: Rep -> Rep from Foreign Builtin qsetrest!: (Rep, Rep) -> Rep qsetrest!(l: Rep, e: Rep): Rep == %store(%tail l,e)$Foreign(Builtin) entireRing? := R has EntireRing --- 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 ** nn == null p => 0 zero? nn => 1 one? nn => 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 nothing => 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 nothing => 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} <>= import NonNegativeInteger import OutputForm )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 ** n == null p => 0 zero? n => 1 one? n => 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) 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 positive? degree lc => 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 nothing => 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 nothing => 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 nothing 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" not one? leadingCoefficient p2 => 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 nothing => 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] 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 nothing => 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} <>= )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} <>= )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): Join(UnivariatePolynomialCategory(R),CoercibleFrom Variable x) with 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} <>= )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} <>= )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} <>= )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 nothing => ["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 not zero? degree(ci) repeat ci:=(ci exquo pi)::P di:=(di exquo pi)::P - differentiate(ci) pi:=gcd(ci,di) i:=i+1 positive? degree(pi) => 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 not zero? degree(di) 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} <>= )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 nothing 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 not one? cont => 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) not zero?(characteristic$RC) => 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] not one? cont => sqp := squareFree cont makeFR(unit(sqp)*coefficient(unit squp,0), append(factorList sqp, pfaclist)) makeFR(coefficient(unit squp,0),pfaclist) @ \section{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} <>= --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. @ <<*>>= <> <> <> <> <> <> <> <> <> <> <> @ \eject \begin{thebibliography}{99} \bibitem{1} nothing \end{thebibliography} \end{document}