\documentclass{article} \usepackage{open-axiom} \begin{document} \title{\$SPAD/src/algebra multpoly.spad} \author{Dave Barton, Barry Trager, James Davenport} \maketitle \begin{abstract} \end{abstract} \eject \tableofcontents \eject \section{domain POLY Polynomial} <>= )abbrev domain POLY Polynomial ++ Author: Dave Barton, Barry Trager ++ Date Created: ++ Date Last Updated: ++ Basic Functions: Ring, degree, eval, coefficient, monomial, differentiate, ++ resultant, gcd ++ Related Constructors: SparseMultivariatePolynomial, MultivariatePolynomial ++ Also See: ++ AMS Classifications: ++ Keywords: polynomial, multivariate ++ References: ++ Description: ++ This type is the basic representation of sparse recursive multivariate ++ polynomials whose variables are arbitrary symbols. The ordering ++ is alphabetic determined by the Symbol type. ++ The coefficient ring may be non commutative, ++ but the variables are assumed to commute. Polynomial(R:Ring): PolynomialCategory(R, IndexedExponents Symbol, Symbol) with if R has Algebra Fraction Integer then integrate: (%, Symbol) -> % ++ integrate(p,x) computes the integral of \spad{p*dx}, i.e. ++ integrates the polynomial p with respect to the variable x. == SparseMultivariatePolynomial(R, Symbol) add import UserDefinedPartialOrdering(Symbol) coerce(p:%):OutputForm == (r:= retractIfCan(p)@Union(R,"failed")) case R => r::R::OutputForm a := userOrdered?() => largest variables p mainVariable(p)::Symbol outputForm(univariate(p, a), a::OutputForm) if R has Algebra Fraction Integer then integrate(p, x) == (integrate univariate(p, x)) (x::%) @ \section{package POLY2 PolynomialFunctions2} <>= )abbrev package POLY2 PolynomialFunctions2 ++ Author: ++ Date Created: ++ Date Last Updated: ++ Basic Functions: ++ Related Constructors: ++ Also See: ++ AMS Classifications: ++ Keywords: ++ References: ++ Description: ++ This package takes a mapping between coefficient rings, and lifts ++ it to a mapping between polynomials over those rings. PolynomialFunctions2(R:Ring, S:Ring): with map: (R -> S, Polynomial R) -> Polynomial S ++ map(f, p) produces a new polynomial as a result of applying ++ the function f to every coefficient of the polynomial p. == add map(f, p) == map(#1::Polynomial(S), f(#1)::Polynomial(S), p)$PolynomialCategoryLifting(IndexedExponents Symbol, Symbol, R, Polynomial R, Polynomial S) @ \section{domain MPOLY MultivariatePolynomial} <>= )abbrev domain MPOLY MultivariatePolynomial ++ Author: Dave Barton, Barry Trager ++ Date Created: ++ Date Last Updated: ++ Basic Functions: Ring, degree, eval, coefficient, monomial, differentiate, ++ resultant, gcd ++ Related Constructors: SparseMultivariatePolynomial, Polynomial ++ Also See: ++ AMS Classifications: ++ Keywords: polynomial, multivariate ++ References: ++ Description: ++ This type is the basic representation of sparse recursive multivariate ++ polynomials whose variables are from a user specified list of symbols. ++ The ordering is specified by the position of the variable in the list. ++ The coefficient ring may be non commutative, ++ but the variables are assumed to commute. MultivariatePolynomial(vl:List Symbol, R:Ring) == SparseMultivariatePolynomial(--SparseUnivariatePolynomial, R, OrderedVariableList vl) @ \section{domain SMP SparseMultivariatePolynomial} <>= )abbrev domain SMP SparseMultivariatePolynomial ++ Author: Dave Barton, Barry Trager ++ Date Created: ++ Date Last Updated: 30 November 1994 ++ Fix History: ++ 30 Nov 94: added gcdPolynomial for float-type coefficients ++ Basic Functions: Ring, degree, eval, coefficient, monomial, differentiate, ++ resultant, gcd ++ Related Constructors: Polynomial, MultivariatePolynomial ++ Also See: ++ AMS Classifications: ++ Keywords: polynomial, multivariate ++ References: ++ Description: ++ This type is the basic representation of sparse recursive multivariate ++ polynomials. It is parameterized by the coefficient ring and the ++ variable set which may be infinite. The variable ordering is determined ++ by the variable set parameter. The coefficient ring may be non-commutative, ++ but the variables are assumed to commute. SparseMultivariatePolynomial(R: Ring,VarSet: OrderedSet): C == T where pgcd ==> PolynomialGcdPackage(IndexedExponents VarSet,VarSet,R,%) C == PolynomialCategory(R,IndexedExponents(VarSet),VarSet) SUP ==> SparseUnivariatePolynomial T == add --constants --D := F(%) replaced by next line until compiler support completed --representations D := SparseUnivariatePolynomial(%) VPoly:= Record(v:VarSet,ts:D) Rep:= Union(R,VPoly) --local function --declarations fn: R -> R n: Integer k: NonNegativeInteger kp:PositiveInteger k1:NonNegativeInteger c: R mvar: VarSet val : R var:VarSet up: D p,p1,p2,pval: % Lval : List(R) Lpval : List(%) Lvar : List(VarSet) --define 0 == 0$R::% 1 == 1$R::% zero? p == p case R and zero?(p)$R one? p == p case R and one?(p)$R -- a local function red(p:%):% == p case R => 0 if ground?(reductum p.ts) then leadingCoefficient(reductum p.ts) else [p.v,reductum p.ts]$VPoly numberOfMonomials(p): NonNegativeInteger == p case R => zero?(p)$R => 0 1 +/[numberOfMonomials q for q in coefficients(p.ts)] coerce(mvar):% == [mvar,monomial(1,1)$D]$VPoly monomial? p == p case R => true sup : D := p.ts 1 ~= numberOfMonomials(sup) => false monomial? leadingCoefficient(sup)$D -- local moreThanOneVariable?: % -> Boolean moreThanOneVariable? p == p case R => false q:=p.ts any?(not ground? #1 ,coefficients q) => true false -- if we already know we use this (slighlty) faster function univariateKnown: % -> SparseUnivariatePolynomial R univariateKnown p == p case R => (leadingCoefficient p) :: SparseUnivariatePolynomial(R) monomial( leadingCoefficient p,degree p.ts)+ univariateKnown(red p) univariate p == p case R =>(leadingCoefficient p) :: SparseUnivariatePolynomial(R) moreThanOneVariable? p => error "not univariate" monomial( leadingCoefficient p,degree p.ts)+ univariate(red p) multivariate (u:SparseUnivariatePolynomial(R),var:VarSet) == ground? u => (leadingCoefficient u) ::% [var,monomial(leadingCoefficient u,degree u)$D]$VPoly + multivariate(reductum u,var) univariate(p:%,mvar:VarSet):SparseUnivariatePolynomial(%) == p case R or mvar>p.v => monomial(p,0)$D pt:=p.ts mvar=p.v => pt monomial(1,p.v,degree pt)*univariate(leadingCoefficient pt,mvar)+ univariate(red p,mvar) -- a local functions, used in next definition unlikeUnivReconstruct(u:SparseUnivariatePolynomial(%),mvar:VarSet):% == zero? (d:=degree u) => coefficient(u,0) monomial(leadingCoefficient u,mvar,d)+ unlikeUnivReconstruct(reductum u,mvar) multivariate(u:SparseUnivariatePolynomial(%),mvar:VarSet):% == ground? u => coefficient(u,0) uu:=u while not zero? uu repeat cc:=leadingCoefficient uu cc case R or mvar > cc.v => uu:=reductum uu return unlikeUnivReconstruct(u,mvar) [mvar,u]$VPoly ground?(p:%):Boolean == p case R => true false -- const p == -- p case R => p -- error "the polynomial is not a constant" monomial(p,mvar,k1) == zero? k1 or zero? p => p p case R or mvar>p.v => [mvar,monomial(p,k1)$D]$VPoly p*[mvar,monomial(1,k1)$D]$VPoly monomial(c:R,e:IndexedExponents(VarSet)):% == zero? e => (c::%) monomial(1,leadingSupport e, leadingCoefficient e) * monomial(c,reductum e) coefficient(p:%, e:IndexedExponents(VarSet)) : R == zero? e => p case R => p::R coefficient(coefficient(p.ts,0),e) p case R => 0 ve := leadingSupport e vp := p.v ve < vp => coefficient(coefficient(p.ts,0),e) ve > vp => 0 coefficient(coefficient(p.ts,leadingCoefficient e),reductum e) -- coerce(e:IndexedExponents(VarSet)) : % == -- e = 0 => 1 -- monomial(1,leadingSupport e, leadingCoefficient e) * -- (reductum e)::% -- retract(p:%):IndexedExponents(VarSet) == -- q:Union(IndexedExponents(VarSet),"failed"):=retractIfCan p -- q :: IndexedExponents(VarSet) -- retractIfCan(p:%):Union(IndexedExponents(VarSet),"failed") == -- p = 0 => degree p -- reductum(p)=0 and leadingCoefficient(p)=1 => degree p -- "failed" coerce(n) == n::R::% coerce(c) == c::% characteristic == characteristic$R recip(p) == p case R => (uu:=recip(p::R);uu case "failed" => "failed"; uu::%) "failed" - p == p case R => -$R p [p.v, - p.ts]$VPoly n * p == p case R => n * p::R mvar:=p.v up:=n*p.ts if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly c * p == c = 1 => p p case R => c * p::R mvar:=p.v up:=c*p.ts if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly p1 + p2 == p1 case R and p2 case R => p1 +$R p2 p1 case R => [p2.v, p1::D + p2.ts]$VPoly p2 case R => [p1.v, p1.ts + p2::D]$VPoly p1.v = p2.v => mvar:=p1.v up:=p1.ts+p2.ts if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly p1.v < p2.v => [p2.v, p1::D + p2.ts]$VPoly [p1.v, p1.ts + p2::D]$VPoly p1 - p2 == p1 case R and p2 case R => p1 -$R p2 p1 case R => [p2.v, p1::D - p2.ts]$VPoly p2 case R => [p1.v, p1.ts - p2::D]$VPoly p1.v = p2.v => mvar:=p1.v up:=p1.ts-p2.ts if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly p1.v < p2.v => [p2.v, p1::D - p2.ts]$VPoly [p1.v, p1.ts - p2::D]$VPoly p1 = p2 == p1 case R => p2 case R => p1 =$R p2 false p2 case R => false p1.v = p2.v => p1.ts = p2.ts false p1 * p2 == p1 case R => p1::R * p2 p2 case R => mvar:=p1.v up:=p1.ts*p2 if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly p1.v = p2.v => mvar:=p1.v up:=p1.ts*p2.ts if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly p1.v > p2.v => mvar:=p1.v up:=p1.ts*p2 if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly --- p1.v < p2.v mvar:=p2.v up:=p1*p2.ts if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly p ** kp == p ** (kp pretend NonNegativeInteger ) p ** k == p case R => p::R ** k -- univariate special case not moreThanOneVariable? p => multivariate( (univariateKnown p) ** k , p.v) mvar:=p.v up:=p.ts ** k if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly if R has IntegralDomain then UnitCorrAssoc ==> Record(unit:%,canonical:%,associate:%) unitNormal(p) == u,c,a:R p case R => (u,c,a):= unitNormal(p::R)$R [u::%,c::%,a::%]$UnitCorrAssoc (u,c,a):= unitNormal(leadingCoefficient(p))$R [u::%,(a*p)::%,a::%]$UnitCorrAssoc unitCanonical(p) == p case R => unitCanonical(p::R)$R (u,c,a):= unitNormal(leadingCoefficient(p))$R a*p unit? p == p case R => unit?(p::R)$R false associates?(p1,p2) == p1 case R => p2 case R and associates?(p1,p2)$R p2 case VPoly and p1.v = p2.v and associates?(p1.ts,p2.ts) if R has approximate then p1 exquo p2 == p1 case R and p2 case R => a:= (p1::R exquo p2::R) if a case "failed" then "failed" else a::% zero? p1 => p1 one? p2 => p1 p1 case R or p2 case VPoly and p1.v < p2.v => "failed" p2 case R or p1.v > p2.v => a:= (p1.ts exquo p2::D) a case "failed" => "failed" [p1.v,a]$VPoly::% -- The next test is useful in the case that R has inexact -- arithmetic (in particular when it is Interval(...)). -- In the case where the test succeeds, empirical evidence -- suggests that it can speed up the computation several times, -- but in other cases where there are a lot of variables -- and p1 and p2 differ only in the low order terms (e.g. p1=p2+1) -- it slows exquo down by about 15-20%. p1 = p2 => 1 a:= p1.ts exquo p2.ts a case "failed" => "failed" mvar:=p1.v up:SUP %:=a if ground? (up) then leadingCoefficient(up) else [mvar,up]$VPoly::% else p1 exquo p2 == p1 case R and p2 case R => a:= (p1::R exquo p2::R) if a case "failed" then "failed" else a::% zero? p1 => p1 one? p2 => p1 p1 case R or p2 case VPoly and p1.v < p2.v => "failed" p2 case R or p1.v > p2.v => a:= (p1.ts exquo p2::D) a case "failed" => "failed" [p1.v,a]$VPoly::% a:= p1.ts exquo p2.ts a case "failed" => "failed" mvar:=p1.v up:SUP %:=a if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly::% map(fn,p) == p case R => fn(p) mvar:=p.v up:=map(map(fn,#1),p.ts) if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly if R has Field then (p : %) / (r : R) == inv(r) * p if R has GcdDomain then content(p) == p case R => p c :R :=0 up:=p.ts while not(zero? up) and not(one? c) repeat c:=gcd(c,content leadingCoefficient(up)) up := reductum up c if R has EuclideanDomain and R has CharacteristicZero and not(R has FloatingPointSystem) then content(p,mvar) == p case R => p gcd(coefficients univariate(p,mvar))$pgcd gcd(p1,p2) == gcd(p1,p2)$pgcd gcd(lp:List %) == gcd(lp)$pgcd gcdPolynomial(a:SUP $,b:SUP $):SUP $ == gcd(a,b)$pgcd else if R has GcdDomain then content(p,mvar) == p case R => p content univariate(p,mvar) gcd(p1,p2) == p1 case R => p2 case R => gcd(p1,p2)$R::% zero? p1 => p2 gcd(p1, content(p2.ts)) p2 case R => zero? p2 => p1 gcd(p2, content(p1.ts)) p1.v < p2.v => gcd(p1, content(p2.ts)) p1.v > p2.v => gcd(content(p1.ts), p2) mvar:=p1.v up:=gcd(p1.ts, p2.ts) if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly if R has FloatingPointSystem then -- eventually need a better notion of gcd's over floats -- this essentially computes the gcds of the monomial contents gcdPolynomial(a:SUP $,b:SUP $):SUP $ == ground? (a) => zero? a => b gcd(leadingCoefficient a, content b)::SUP $ ground?(b) => zero? b => b gcd(leadingCoefficient b, content a)::SUP $ conta := content a mona:SUP $ := monomial(conta, minimumDegree a) if not one? mona then a := (a exquo mona)::SUP $ contb := content b monb:SUP $ := monomial(contb, minimumDegree b) if not one? monb then b := (b exquo monb)::SUP $ mong:SUP $ := monomial(gcd(conta, contb), min(degree mona, degree monb)) degree(a) >= degree b => not((a exquo b) case "failed") => mong * b mong not((b exquo a) case "failed") => mong * a mong coerce(p):OutputForm == p case R => (p::R)::OutputForm outputForm(p.ts,p.v::OutputForm) coefficients p == p case R => list(p :: R)$List(R) "append"/[coefficients(p1)$% for p1 in coefficients(p.ts)] retract(p:%):R == p case R => p :: R error "cannot retract nonconstant polynomial" retractIfCan(p:%):Union(R, "failed") == p case R => p::R "failed" -- leadingCoefficientRecursive(p:%):% == -- p case R => p -- leadingCoefficient p.ts mymerge:(List VarSet,List VarSet) ->List VarSet mymerge(l:List VarSet,m:List VarSet):List VarSet == empty? l => m empty? m => l first l = first m => empty? rest l => setrest!(l,rest m) l empty? rest m => l setrest!(l, mymerge(rest l, rest m)) l first l > first m => empty? rest l => setrest!(l,m) l setrest!(l, mymerge(rest l, m)) l empty? rest m => setrest!(m,l) m setrest!(m,mymerge(l,rest m)) m variables p == p case R => empty() lv:List VarSet:=empty() q := p.ts while not zero? q repeat lv:=mymerge(lv,variables leadingCoefficient q) q := reductum q cons(p.v,lv) mainVariable p == p case R => "failed" p.v eval(p,mvar,pval) == univariate(p,mvar)(pval) eval(p,mvar,val) == univariate(p,mvar)(val) evalSortedVarlist(p,Lvar,Lpval):% == p case R => p empty? Lvar or empty? Lpval => p mvar := Lvar.first mvar > p.v => evalSortedVarlist(p,Lvar.rest,Lpval.rest) pval := Lpval.first pts := map(evalSortedVarlist(#1,Lvar,Lpval),p.ts) mvar=p.v => pval case R => pts (pval::R) pts pval multivariate(pts,p.v) eval(p,Lvar,Lpval) == empty? rest Lvar => evalSortedVarlist(p,Lvar,Lpval) sorted?(#1 > #2, Lvar) => evalSortedVarlist(p,Lvar,Lpval) nlvar := sort(#1 > #2,Lvar) nlpval := Lvar = nlvar => Lpval [Lpval.position(mvar,Lvar) for mvar in nlvar] evalSortedVarlist(p,nlvar,nlpval) eval(p,Lvar,Lval) == eval(p,Lvar,[val::% for val in Lval]$(List %)) -- kill? degree(p,mvar) == p case R => 0 mvar= p.v => degree p.ts mvar > p.v => 0 -- might as well take advantage of the order max(degree(leadingCoefficient p.ts,mvar),degree(red p,mvar)) degree(p,Lvar) == [degree(p,mvar) for mvar in Lvar] degree p == p case R => 0 degree(leadingCoefficient(p.ts)) + monomial(degree(p.ts), p.v) minimumDegree p == p case R => 0 md := minimumDegree p.ts minimumDegree(coefficient(p.ts,md)) + monomial(md, p.v) minimumDegree(p,mvar) == p case R => 0 mvar = p.v => minimumDegree p.ts md:=minimumDegree(leadingCoefficient p.ts,mvar) zero? (p1:=red p) => md min(md,minimumDegree(p1,mvar)) minimumDegree(p,Lvar) == [minimumDegree(p,mvar) for mvar in Lvar] totalDegree(p, Lvar) == ground? p => 0 null setIntersection(Lvar, variables p) => 0 u := univariate(p, mv := mainVariable(p)::VarSet) weight:NonNegativeInteger := (member?(mv,Lvar) => 1; 0) tdeg:NonNegativeInteger := 0 while u ~= 0 repeat termdeg:NonNegativeInteger := weight*degree u + totalDegree(leadingCoefficient u, Lvar) tdeg := max(tdeg, termdeg) u := reductum u tdeg if R has CommutativeRing then differentiate(p,mvar) == p case R => 0 mvar=p.v => up:=differentiate p.ts if ground? up then leadingCoefficient(up) else [mvar,up]$VPoly up:=map(differentiate(#1,mvar),p.ts) if ground? up then leadingCoefficient(up) else [p.v,up]$VPoly leadingCoefficient(p) == p case R => p leadingCoefficient(leadingCoefficient(p.ts)) -- trailingCoef(p) == -- p case R => p -- coef(p.ts,0) case R => coef(p.ts,0) -- trailingCoef(red p) -- TrailingCoef(p) == trailingCoef(p) leadingMonomial p == p case R => p monomial(leadingMonomial leadingCoefficient(p.ts), p.v, degree(p.ts)) reductum(p) == p case R => 0 p - leadingMonomial p -- if R is Integer then -- pgcd := PolynomialGcdPackage(%,VarSet) -- gcd(p1,p2) == -- gcd(p1,p2)$pgcd -- -- else if R is RationalNumber then -- gcd(p1,p2) == -- mrat:= MRationalFactorize(VarSet,%) -- gcd(p1,p2)$mrat -- -- else gcd(p1,p2) == -- p1 case R => -- p2 case R => gcd(p1,p2)$R::% -- p1 = 0 => p2 -- gcd(p1, content(p2.ts)) -- p2 case R => -- p2 = 0 => p1 -- gcd(p2, content(p1.ts)) -- p1.v < p2.v => gcd(p1, content(p2.ts)) -- p1.v > p2.v => gcd(content(p1.ts), p2) -- PSimp(p1.v, gcd(p1.ts, p2.ts)) @ \section{domain INDE IndexedExponents} <>= )abbrev domain INDE IndexedExponents ++ Author: James Davenport ++ Date Created: ++ Date Last Updated: ++ Basic Functions: ++ Related Constructors: ++ Also See: ++ AMS Classifications: ++ Keywords: ++ References: ++ Description: ++ IndexedExponents of an ordered set of variables gives a representation ++ for the degree of polynomials in commuting variables. It gives an ordered ++ pairing of non negative integer exponents with variables IndexedExponents(Varset:OrderedSet): C == T where C == Join(OrderedAbelianMonoidSup, IndexedDirectProductCategory(NonNegativeInteger,Varset)) T == IndexedDirectProductOrderedAbelianMonoidSup(NonNegativeInteger,Varset) add Term:= Record(k:Varset,c:NonNegativeInteger) Rep:= List Term coerceOF(t: Term):OutputForm == -- converts term to OutputForm t.c = 1 => (t.k)::OutputForm (t.k)::OutputForm ** (t.c)::OutputForm coerce(x: %):OutputForm == -- converts entire exponents to OutputForm null x => 1::Integer::OutputForm null rest x => coerceOF(first x) reduce("*",[coerceOF t for t in x]) @ \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}