\documentclass{article} \usepackage{open-axiom} \begin{document} \title{\$SPAD/src/algebra gpgcd.spad} \author{The Axiom Team} \maketitle \begin{abstract} \end{abstract} \eject \tableofcontents \eject \section{package GENPGCD GeneralPolynomialGcdPackage} <>= )abbrev package GENPGCD GeneralPolynomialGcdPackage ++ Description: ++ This package provides operations for GCD computations ++ on polynomials GeneralPolynomialGcdPackage(E,OV,R,P):C == T where R : PolynomialFactorizationExplicit P : PolynomialCategory(R,E,OV) OV : OrderedSet E : OrderedAbelianMonoidSup SUPP ==> SparseUnivariatePolynomial P --JHD ContPrim ==> Record(cont:P,prim:P) C == with gcdPolynomial : (SUPP,SUPP) -> SUPP ++ gcdPolynomial(p,q) returns the GCD of p and q randomR : () ->R ++ randomR() should be local but conditional --JHD gcd : (P,P) -> P --JHD gcd : List P -> P --JHD gcdprim : (P,P) -> P --JHD gcdprim : List P -> P --JHD gcdcofact : List P -> List P --JHD gcdcofactprim : List P -> List P --JHD primitate : (P,OV) -> P --JHD primitate : SUPP -> SUPP --JHD content : P -> P --JHD content : List P -> List P --JHD contprim : List P -> List ContPrim --JHD monomContent : (P,OV) -> P --JHD monomContent : SUPP -> SUPP T == add SUPR ==> SparseUnivariatePolynomial R --JHD SUPLGcd ==> Record(locgcd:SUPP,goodint:List R) --JHD LGcd ==> Record(locgcd:P,goodint:List R) --JHD UTerm ==> Record(lpol:List SUPR,lint:List R,mpol:P) --JHD--JHD pmod:R := (prevPrime(2**26)$IntegerPrimesPackage(Integer))::R --JHD import MultivariateLifting(E,OV,R,P,pmod) import UnivariatePolynomialCategoryFunctions2(R,SUPR,P,SUPP) import UnivariatePolynomialCategoryFunctions2(P,SUPP,R,SUPR) -------- Local Functions -------- --JHD abs : P -> P better : (P,P) -> Boolean --JHD failtest : (P,P,P) -> Boolean --JHD gcdMonom : (P,P,OV) -> P --JHD gcdTermList : (P,P) -> P --JHD gcdPrim : (P,P,OV) -> P --JHD gcdSameMainvar : (P,P,OV) -> P --JHD internal : (P,P,OV) -> P --JHD good : (P,List OV) -> Record(upol:SUPR,inval:List R) --JHD gcdPrs : (P,P,NNI,OV) -> Union(P,"failed") --JHD --JHD chooseVal : (P,P,List OV) -> UTerm --JHD localgcd : (P,P,List OV) -> LGcd --JHD notCoprime : (P,P, List NNI,List OV) -> P --JHD imposelc : (List SUPR,List OV,List R,List P) -> List SUPR --JHD lift? :(P,P,UTerm,List NNI,List OV) -> Union("failed",P) -- lift :(P,SUPR,SUPR,P,List OV,List NNI,List R) -> P lift : (SUPR,SUPP,SUPR,List OV,List R) -> Union(SUPP,"failed") -- lifts first and third arguments as factors of the second -- fourth is number of variables. --JHD monomContent : (P,OV) -> P monomContentSup : SUPP -> SUPP -- --JHD gcdcofact : List P -> List P gcdTrivial : (SUPP,SUPP) -> SUPP gcdSameVariables: (SUPP,SUPP,List OV) -> SUPP recursivelyGCDCoefficients: (SUPP,List OV,SUPP,List OV) -> SUPP flatten : (SUPP,List OV) -> SUPP -- evaluates out all variables in the second -- argument, leaving a polynomial of the same -- degree -- eval : (SUPP,List OV,List R) -> SUPP variables : SUPP -> List OV ---- JHD's exported functions --- gcdPolynomial(p1:SUPP,p2:SUPP) == zero? p1 => p2 zero? p2 => p1 0=degree p1 => gcdTrivial(p1,p2) 0=degree p2 => gcdTrivial(p2,p1) if degree p1 < degree p2 then (p1,p2):=(p2,p1) p1 exquo p2 case SUPP => (unitNormal p2).canonical c1:= monomContentSup(p1) c2:= monomContentSup(p2) p1:= (p1 exquo c1)::SUPP p2:= (p2 exquo c2)::SUPP (p1 exquo p2) case SUPP => (unitNormal p2).canonical * gcd(c1,c2) vp1:=variables p1 vp2:=variables p2 v1:=setDifference(vp1,vp2) v2:=setDifference(vp2,vp1) #v1 = 0 and #v2 = 0 => gcdSameVariables(p1,p2,vp1)*gcd(c1,c2) -- all variables are in common v:=setDifference(vp1,v1) pp1:=flatten(p1,v1) pp2:=flatten(p2,v2) g:=gcdSameVariables(pp1,pp2,v) one? g => gcd(c1,c2)::SUPP (#v1 = 0 or not (p1 exquo g) case "failed") and -- if #vi = 0 then pp1 = p1, so we know g divides (#v2 = 0 or not (p2 exquo g) case "failed") => g*gcd(c1,c2) -- divdes them both, so is the gcd -- OK, so it's not the gcd: try again v:=variables g -- there can be at most these variables in answer v1:=setDifference(vp1,v) v2:=setDifference(vp2,v) if (#v1 = 0) then g:= gcdSameVariables(g,flatten(p2,v2),v) else if (#v2=0) then g:=gcdSameVariables(g,flatten(p1,v1),v) else g:=gcdSameVariables(g,flatten(p1,v1)-flatten(p2,v2),v) one? g => gcd(c1,c2)::SUPP (#v1 = 0 or not (p1 exquo g) case "failed") and (#v2 = 0 or not (p2 exquo g) case "failed") => g*gcd(c1,c2)::SUPP -- divdes them both, so is the gcd v:=variables g -- there can be at most these variables in answer v1:=setDifference(vp1,v) if #v1 ~= 0 then g:=recursivelyGCDCoefficients(g,v,p1,v1) one? g => return gcd(c1,c2)::SUPP v:=variables g -- there can be at most these variables in answer v2:=setDifference(vp2,v) recursivelyGCDCoefficients(g,v,p2,v2)*gcd(c1,c2) if R has StepThrough then randomCount:R := init() randomR() == (v:=nextItem(randomCount)) case R => randomCount:=v v SAY("Taking next stepthrough range in GeneralPolynomialGcdPackage")$Lisp randomCount:=init() randomCount else randomR() == (random(100)$Integer)::R ---- JHD's local functions --- gcdSameVariables(p1:SUPP,p2:SUPP,lv:List OV) == -- two non-trivial primitive (or, at least, we don't care -- about content) -- polynomials with precisely the same degree #lv = 0 => map(#1::P,gcdPolynomial(map(ground,p1), map(ground,p2))) degree p2 = 1 => p1 exquo p2 case SUPP => p2 1 gcdLC:=gcd(leadingCoefficient p1,leadingCoefficient p2) lr:=[randomR() for vv in lv] count:NonNegativeInteger:=0 while count<10 repeat while zero? eval(gcdLC,lv,lr) and count<10 repeat lr:=[randomR() for vv in lv] count:=count+1 count = 10 => error "too many evaluations in GCD code" up1:SUPR:=map(ground eval(#1,lv,lr),p1) up2:SUPR:=map(ground eval(#1,lv,lr),p2) u:=gcdPolynomial(up1,up2) degree u = 0 => return 1 -- let's pick a second one, just to check lrr:=[randomR() for vv in lv] while zero? eval(gcdLC,lv,lrr) and count<10 repeat lrr:=[randomR() for vv in lv] count:=count+1 count = 10 => error "too many evaluations in GCD code" vp1:SUPR:=map(ground eval(#1,lv,lrr),p1) vp2:SUPR:=map(ground eval(#1,lv,lrr),p2) v:=gcdPolynomial(vp1,vp2) degree v = 0 => return 1 if degree v < degree u then u:=v up1:=vp1 up2:=vp2 lr:=lrr up1:=(up1 exquo u)::SUPR degree gcd(u,up1) = 0 => ans:=lift(u,p1,up1,lv,lr) ans case SUPP => return ans "next" up2:=(up2 exquo u)::SUPR degree gcd(u,up2) = 0 => ans:=lift(u,p2,up2,lv,lr) ans case SUPP => return ans "next" -- so neither cofactor is relatively prime count:=0 while count < 10 repeat r:=randomR() uu:=up1+r*up2 degree gcd(u,uu)=0 => ans:= lift(u,p1+r::P *p2,uu,lv,lr) ans case SUPP => return ans "next" error "too many evaluations in GCD code" count >= 10 => error "too many evaluations in GCD code" lift(gR:SUPR,p:SUPP,cfR:SUPR,lv:List OV,lr:List R) == -- lift the coprime factorisation gR*cfR = (univariate of p) -- where the variables lv have been evaluated at lr lcp:=leadingCoefficient p g:=monomial(lcp,degree gR)+map(#1::P,reductum gR) cf:=monomial(lcp,degree cfR)+map(#1::P,reductum cfR) p:=lcp*p -- impose leaidng coefficient of p on each factor while lv ~= [] repeat v:=first lv r:=first lr lv:=rest lv lr:=rest lr thisp:=map(eval(#1,lv,lr),p) d:="max"/[degree(c,v) for c in coefficients p] prime:=v::P - r::P pn:=prime origFactors:=[g,cf]::List SUPP for n in 1..d repeat Ecart:=(thisp- g*cf) exquo pn Ecart case "failed" => error "failed lifting in hensel in Complex Polynomial GCD" zero? Ecart => leave step:=solveLinearPolynomialEquation(origFactors, map(eval(#1,v,r),Ecart::SUPP)) step case "failed" => return "failed" g:=g+pn*first step cf:=cf+pn*second step pn:=pn*prime thisp ~= g*cf => return "failed" g recursivelyGCDCoefficients(g:SUPP,v:List OV,p:SUPP,pv:List OV) == mv:=first pv -- take each coefficient w.r.t. mv pv:=rest pv -- and recurse on pv as necessary d:="max"/[degree(u,mv) for u in coefficients p] for i in 0..d repeat p1:=map(coefficient(#1,mv,i),p) oldg:=g if pv = [] then g:=gcdSameVariables(g,p1,v) else g:=recursivelyGCDCoefficients(p,v,p1,pv) one? g => return 1 g~=oldg => oldv:=v v:=variables g pv:=setUnion(pv,setDifference(v,oldv)) g flatten(p1:SUPP,lv:List OV) == #lv = 0 => p1 lr:=[ randomR() for vv in lv] dg:=degree p1 ans : SUPP while dg ~= degree (ans:= map(eval(#1,lv,lr),p1)) repeat lr:=[ randomR() for vv in lv] ans -- eval(p1:SUPP,lv:List OV,lr:List R) == map(eval(#1,lv,lr),p1) variables(p1:SUPP) == removeDuplicates ("concat"/[variables u for u in coefficients p1]) gcdTrivial(p1:SUPP,p2:SUPP) == -- p1 is non-zero, but has degree zero -- p2 is non-zero cp1:=leadingCoefficient p1 one? cp1 => 1 degree p2 = 0 => gcd(cp1,leadingCoefficient p2)::SUPP un?:=unit? cp1 while not zero? p2 and not un? repeat cp1:=gcd(leadingCoefficient p2,cp1) un?:=unit? cp1 p2:=reductum p2 un? => 1 cp1::SUPP ---- Local functions ---- --JHD -- test if something wrong happened in the gcd --JHD failtest(f:P,p1:P,p2:P) : Boolean == --JHD (p1 exquo f) case "failed" or (p2 exquo f) case "failed" --JHD --JHD -- Choose the integers --JHD chooseVal(p1:P,p2:P,lvar:List OV):UTerm == --JHD x:OV:=lvar.first --JHD lvr:=lvar.rest --JHD d1:=degree(p1,x) --JHD d2:=degree(p2,x) --JHD dd:NNI:=0$NNI --JHD nvr:NNI:=#lvr --JHD lval:List R :=[] --JHD range:I:=8 --JHD for i in 1.. repeat --JHD range:=2*range --JHD lval:=[(random(2*range)$I - range)::R for i in 1..nvr] --JHD uf1:SUPR:=univariate eval(p1,lvr,lval) --JHD degree uf1 ~= d1 => "new point" --JHD uf2:SUPR:=univariate eval(p2,lvr,lval) --JHD degree uf2 ~= d2 => "new point" --JHD u:=gcd(uf1,uf2) --JHD du:=degree u --JHD --the univariate gcd is 1 --JHD if du=0 then return [[1$SUPR],lval,0$P]$UTerm --JHD --JHD ugcd:List SUPR:=[u,(uf1 exquo u)::SUPR,(uf2 exquo u)::SUPR] --JHD uterm:=[ugcd,lval,0$P]$UTerm --JHD dd=0 => dd:=du --JHD --JHD --the degree is not changed --JHD du=dd => --JHD --JHD --test if one of the polynomials is the gcd --JHD dd=d1 => --JHD if ^((f:=p2 exquo p1) case "failed") then --JHD return [[u],lval,p1]$UTerm --JHD if dd~=d2 then dd:=(dd-1)::NNI --JHD --JHD dd=d2 => --JHD if ^((f:=p1 exquo p2) case "failed") then --JHD return [[u],lval,p2]$UTerm --JHD dd:=(dd-1)::NNI --JHD return uterm --JHD --JHD --the new gcd has degree less --JHD du
dd:=du --JHD --JHD good(f:P,lvr:List OV):Record(upol:SUPR,inval:List R) == --JHD nvr:NNI:=#lvr --JHD range:I:=1 --JHD ltry:List List R:=[] --JHD while true repeat --JHD range:=2*range --JHD lval:=[(random(2*range)$I -range)::R for i in 1..nvr] --JHD member?(lval,ltry) => "new point" --JHD ltry:=cons(lval,ltry) --JHD uf:=univariate eval(f,lvr,lval) --JHD if degree gcd(uf,differentiate uf)=0 then return [uf,lval] --JHD --JHD -- impose the right lc --JHD imposelc(lipol:List SUPR, --JHD lvar:List OV,lval:List R,leadc:List P):List SUPR == --JHD result:List SUPR :=[] --JHD lvar:=lvar.rest --JHD for pol in lipol for leadpol in leadc repeat --JHD p1:= univariate eval(leadpol,lvar,lval) * pol --JHD result:= cons((p1 exquo leadingCoefficient pol)::SUPR,result) --JHD reverse result --JHD --JHD --Compute the gcd between not coprime polynomials --JHD notCoprime(g:P,p2:P,ldeg:List NNI,lvar:List OV) : P == --JHD x:OV:=lvar.first --JHD lvar1:List OV:=lvar.rest --JHD lg1:=gcdcofact([g,differentiate(g,x)]) --JHD g1:=lg1.1 --JHD lg:LGcd:=localgcd(g1,p2,lvar) --JHD (l,lval):=(lg.locgcd,lg.goodint) --JHD p2:=(p2 exquo l)::P --JHD (gd1,gd2):=(l,l) --JHD ul:=univariate(eval(l,lvar1,lval)) --JHD dl:=degree ul --JHD if degree gcd(ul,differentiate ul) ~=0 then --JHD newchoice:=good(l,lvar.rest) --JHD ul:=newchoice.upol --JHD lval:=newchoice.inval --JHD ug1:=univariate(eval(g1,lvar1,lval)) --JHD ulist:=[ug1,univariate eval(p2,lvar1,lval)] --JHD lcpol:=[leadingCoefficient univariate(g1,x), --JHD leadingCoefficient univariate(p2,x)] --JHD while true repeat --JHD d:SUPR:=gcd(cons(ul,ulist)) --JHD if degree d =0 then return gd1 --JHD lquo:=(ul exquo d)::SUPR --JHD if degree lquo ~=0 then --JHD lgcd:=gcd(cons(leadingCoefficient univariate(l,x),lcpol)) --JHD gd2:=lift(l,d,lquo,lgcd,lvar,ldeg,lval) --JHD l:=gd2 --JHD ul:=univariate(eval(l,lvar1,lval)) --JHD dl:=degree ul --JHD gd1:=gd1*gd2 --JHD ulist:=[(uf exquo d)::SUPR for uf in ulist] --JHD --JHD -- we suppose that the poly have the same mainvar, deg p1#2,setUnion(variables p1,variables p2)) --JHD d1:=degree(p1,x) --JHD d2:=degree(p2,x) --JHD result: P:=localgcd(p1,p2,lvar).locgcd --JHD -- special cases --JHD result=1 => 1$P --JHD (dr:=degree(result,x))=d1 or dr=d2 => result --JHD while failtest(result,p1,p2) repeat --JHD SAY$Lisp "retrying gcd" --JHD result:=localgcd(p1,p2,lvar).locgcd --JHD result --JHD --JHD --local function for the gcd : it returns the evaluation point too --JHD localgcd(p1:P,p2:P,lvar:List(OV)) : LGcd == --JHD x:OV:=lvar.first --JHD uterm:=chooseVal(p1,p2,lvar) --JHD listpol:= uterm.lpol --JHD ud:=listpol.first --JHD dd:= degree ud --JHD --JHD --the univariate gcd is 1 --JHD dd=0 => [1$P,uterm.lint]$LGcd --JHD --JHD --one of the polynomials is the gcd --JHD dd=degree(p1,x) or dd=degree(p2,x) => --JHD [uterm.mpol,uterm.lint]$LGcd --JHD ldeg:List NNI:=map(min,degree(p1,lvar),degree(p2,lvar)) --JHD --JHD -- if there is a polynomial g s.t. g/gcd and gcd are coprime ... --JHD -- I can lift --JHD (h:=lift?(p1,p2,uterm,ldeg,lvar)) case "failed" => --JHD [notCoprime(p1,p2,ldeg,lvar),uterm.lint]$LGcd --JHD [h::P,uterm.lint]$LGcd --JHD --JHD --JHD -- content, internal functions return the poly if it is a monomial --JHD monomContent(p:P,var:OV):P == --JHD ground? p => 1$P --JHD md:= minimumDegree(p,var) --JHD ((var::P)**md)*(gcd sort(better,coefficients univariate(p,var))) monomContentSup(u:SUPP):SUPP == degree(u) = 0$NonNegativeInteger => 1$SUPP md:= minimumDegree u gcd(sort(better,coefficients u)) * monomial(1$P,md)$SUPP --JHD -- change the polynomials to have positive lc --JHD abs(p:P): P == unitNormal(p).canonical -- Ordering for gcd purposes better(p1:P,p2:P):Boolean == ground? p1 => true ground? p2 => false degree(p1,mainVariable(p1)::OV) < degree(p2,mainVariable(p2)::OV) -- PRS algorithm -- gcdPrs(p1:P,p2:P,d:NNI,var:OV):Union(P,"failed") == -- u1:= univariate(p1,var) -- u2:= univariate(p2,var) -- finished:Boolean:= false -- until finished repeat -- dd:NNI:=(degree u1 - degree u2)::NNI -- lc1:SUPP:=leadingCoefficient u2 * reductum u1 -- lc2:SUPP:=leadingCoefficient u1 * reductum u2 -- u3:SUPP:= primitate((lc1-lc2)*monomial(1$P,dd))$% -- (d3:=degree(u3)) <= d => finished:= true -- u1:= u2 -- u2:= u3 -- if d3 > degree(u1) then (u1,u2):= (u2,u1) -- g:= (u2 exquo u3) -- g case SUPP => abs multivariate(u3,var) -- "failed" -- Gcd between polynomial p1 and p2 with -- mainVariable p1 < x=mainVariable p2 --JHD gcdTermList(p1:P,p2:P) : P == --JHD termList:=sort(better, --JHD cons(p1,coefficients univariate(p2,(mainVariable p2)::OV))) --JHD q:P:=termList.first --JHD for term in termList.rest until q = 1$P repeat q:= gcd(q,term) --JHD q --JHD --JHD -- Gcd between polynomials with the same mainVariable --JHD gcdSameMainvar(p1:P,p2:P,mvar:OV): P == --JHD if degree(p1,mvar) < degree(p2,mvar) then (p1,p2):= (p2,p1) --JHD (p1 exquo p2) case P => abs p2 --JHD c1:= monomContent(p1,mvar)$% --JHD c1 = p1 => gcdMonom(p1,p2,mvar) --JHD c2:= monomContent(p2,mvar)$% --JHD c2 = p2 => gcdMonom(p2,p1,mvar) --JHD p1:= (p1 exquo c1)::P --JHD p2:= (p2 exquo c2)::P --JHD if degree(p1,mvar) < degree(p2,mvar) then (p1,p2):= (p2,p1) --JHD (p1 exquo p2) case P => abs(p2) * gcd(c1,c2) --JHD abs(gcdPrim(p1,p2,mvar)) * gcd(c1,c2) --JHD --JHD -- make the polynomial primitive with respect to var --JHD primitate(p:P,var:OV):P == (p exquo monomContent(p,var))::P --JHD --JHD primitate(u:SUPP):SUPP == (u exquo monomContentSup u)::SUPP --JHD --JHD -- gcd between primitive polynomials with the same mainVariable --JHD gcdPrim(p1:P,p2:P,mvar:OV):P == --JHD vars:= removeDuplicates append(variables p1,variables p2) --JHD #vars=1 => multivariate(gcd(univariate p1,univariate p2),mvar) --JHD vars:=delete(vars,position(mvar,vars)) --JHD --d:= degModGcd(p1,p2,mvar,vars) --JHD --d case "failed" => internal(p2,p1,mvar) --JHD --deg:= d:NNI --JHD --deg = 0$NNI => 1$P --JHD --deg = degree(p1,mvar) => --JHD -- (p2 exquo p1) case P => abs(p1) -- already know that --JHD -- ^(p1 exquo p2) --JHD -- internal(p2,p1,mvar) --JHD --cheapPrs?(p1,p2,deg,mvar) => --JHD -- g:= gcdPrs(p1,p2,deg,mvar) --JHD -- g case P => g::P --JHD -- internal(p2,p1,mvar) --JHD internal(p2,p1,mvar) --JHD --JHD -- gcd between a monomial and a polynomial --JHD gcdMonom(m:P,p:P,var:OV):P == --JHD ((var::P) ** min(minimumDegree(m,var),minimumDegree(p,var))) * --JHD gcdTermList(leadingCoefficient(univariate(m,var)),p) --JHD --JHD --If there is a pol s.t. pol/gcd and gcd are coprime I can lift --JHD lift?(p1:P,p2:P,uterm:UTerm,ldeg:List NNI, --JHD lvar:List OV) : Union("failed",P) == --JHD x:OV:=lvar.first --JHD leadpol:Boolean:=false --JHD (listpol,lval):=(uterm.lpol,uterm.lint) --JHD d:=listpol.first --JHD listpol:=listpol.rest --JHD nolift:Boolean:=true --JHD for uf in listpol repeat --JHD --note uf and d not necessarily primitive --JHD degree gcd(uf,d) =0 => nolift:=false --JHD nolift => "failed" --JHD f:P:=([p1,p2]$List(P)).(position(uf,listpol)) --JHD lgcd:=gcd(leadingCoefficient univariate(p1,x), --JHD leadingCoefficient univariate(p2,x)) --JHD lift(f,d,uf,lgcd,lvar,ldeg,lval) --JHD --JHD -- interface with the general "lifting" function --JHD lift(f:P,d:SUPR,uf:SUPR,lgcd:P,lvar:List OV, --JHD ldeg:List NNI,lval:List R):P == --JHD x:OV:=lvar.first --JHD leadpol:Boolean:=false --JHD lcf:P --JHD lcf:=leadingCoefficient univariate(f,x) --JHD df:=degree(f,x) --JHD leadlist:List(P):=[] --JHD --JHD if lgcd~=1$P then --JHD leadpol:=true --JHD f:=lgcd*f --JHD ldeg:=[n0+n1 for n0 in ldeg for n1 in degree(lgcd,lvar)] --JHD lcd:R:=leadingCoefficient d --JHD if ground? lgcd then d:=((retract lgcd) *d exquo lcd)::SUPR --JHD else d:=(retract(eval(lgcd,lvar.rest,lval)) * d exquo lcd)::SUPR --JHD uf:=lcd*uf --JHD leadlist:=[lgcd,lcf] --JHD lg:=imposelc([d,uf],lvar,lval,leadlist) --JHD plist:=lifting(univariate(f,x),lvar,lg,lval,leadlist,ldeg)::List P --JHD (p0:P,p1:P):=(plist.first,plist.2) --JHD if univariate eval(p0,rest lvar,lval) ~= lg.first then --JHD (p0,p1):=(p1,p0) --JHD ^leadpol => p0 --JHD cprim:=contprim([p0]) --JHD cprim.first.prim --JHD --JHD -- Gcd for two multivariate polynomials --JHD gcd(p1:P,p2:P) : P == --JHD (p1:= abs(p1)) = (p2:= abs(p2)) => p1 --JHD ground? p1 => --JHD p1 = 1$P => p1 --JHD p1 = 0$P => p2 --JHD ground? p2 => gcd((retract p1)@R,(retract p2)@R)::P --JHD gcdTermList(p1,p2) --JHD ground? p2 => --JHD p2 = 1$P => p2 --JHD p2 = 0$P => p1 --JHD gcdTermList(p2,p1) --JHD mv1:= mainVariable(p1)::OV --JHD mv2:= mainVariable(p2)::OV --JHD mv1 = mv2 => gcdSameMainvar(p1,p2,mv1) --JHD mv1 < mv2 => gcdTermList(p1,p2) --JHD gcdTermList(p2,p1) --JHD --JHD -- Gcd for a list of multivariate polynomials --JHD gcd(listp:List P) : P == --JHD lf:=sort(better,listp) --JHD f:=lf.first --JHD for g in lf.rest repeat --JHD f:=gcd(f,g) --JHD if f=1$P then return f --JHD f --JHD -- Gcd and cofactors for a list of polynomials --JHD gcdcofact(listp : List P) : List P == --JHD h:=gcd listp --JHD cons(h,[(f exquo h) :: P for f in listp]) --JHD --JHD -- Gcd for primitive polynomials --JHD gcdprim(p1:P,p2:P):P == --JHD (p1:= abs(p1)) = (p2:= abs(p2)) => p1 --JHD ground? p1 => --JHD ground? p2 => gcd((retract p1)@R,(retract p2)@R)::P --JHD p1 = 0$P => p2 --JHD 1$P --JHD ground? p2 => --JHD p2 = 0$P => p1 --JHD 1$P --JHD mv1:= mainVariable(p1)::OV --JHD mv2:= mainVariable(p2)::OV --JHD mv1 = mv2 => --JHD md:=min(minimumDegree(p1,mv1),minimumDegree(p2,mv1)) --JHD mp:=1$P --JHD if md>1 then --JHD mp:=(mv1::P)**md --JHD p1:=(p1 exquo mp)::P --JHD p2:=(p2 exquo mp)::P --JHD mp*gcdPrim(p1,p2,mv1) --JHD 1$P --JHD --JHD -- Gcd for a list of primitive multivariate polynomials --JHD gcdprim(listp:List P) : P == --JHD lf:=sort(better,listp) --JHD f:=lf.first --JHD for g in lf.rest repeat --JHD f:=gcdprim(f,g) --JHD if f=1$P then return f --JHD f --JHD -- Gcd and cofactors for a list of primitive polynomials --JHD gcdcofactprim(listp : List P) : List P == --JHD h:=gcdprim listp --JHD cons(h,[(f exquo h) :: P for f in listp]) --JHD --JHD -- content of a polynomial (with respect to its main var) --JHD content(f:P):P == --JHD ground? f => f --JHD x:OV:=(mainVariable f)::OV --JHD gcd sort(better,coefficients univariate(f,x)) --JHD --JHD -- contents of a list of polynomials --JHD content(listf:List P) : List P == [content f for f in listf] --JHD --JHD -- contents and primitive parts of a list of polynomials --JHD contprim(listf:List P) : List ContPrim == --JHD prelim :List P := content listf --JHD [[q,(f exquo q)::P]$ContPrim for q in prelim for f in listf] --JHD @ \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}