diff --git a/books/bookvol10.4.pamphlet b/books/bookvol10.4.pamphlet index 603c8ec..5a0db45 100644 --- a/books/bookvol10.4.pamphlet +++ b/books/bookvol10.4.pamphlet @@ -139486,8 +139486,8 @@ constraint the remainders will have huge numerators and denominators. reverse result --Compute the gcd between not coprime polynomials - notCoprime(g:SUPP,p2:SUPP,ldeg:List NNI,_ - lvar1:List OV,ltry:List List R) : SUPP == + notCoprime(g:SUPP, p2:SUPP, ldeg:List NNI,_ + lvar1:List OV, ltry:List List R) : SUPP == g1:=gcd(g,differentiate g) l1 := (g exquo g1)::SUPP lg:LGcd:=localgcd(l1,p2,lvar1,ltry) @@ -139579,10 +139579,21 @@ constraint the remainders will have huge numerators and denominators. ground? p2 => false degree(p1,mainVariable(p1)::OV) < degree(p2,mainVariable(p2)::OV) + best_to_front(l : List P) : List P == + ress := [] + best := first(l) + for p in rest l repeat + if better(p, best) then + ress := cons(best, ress) + best := p + else + ress := cons(p, ress) + cons(best, ress) + -- Gcd between polynomial p1 and p2 with -- mainVariable p1 < x=mainVariable p2 gcdTermList(p1:P,p2:P) : P == - termList:=sort(better, + termList := best_to_front( cons(p1,coefficients univariate(p2,(mainVariable p2)::OV))) q:P:=termList.first for term in termList.rest until q = 1$P repeat q:= gcd(q,term) @@ -139617,52 +139628,58 @@ by iteratively ``lifting'' the solution modulo successive powers of $p$. See Volume 10.1 for more details. \begin{chunk}{package PGCD PolynomialGcdPackage} + --If there is a pol s.t. pol/gcd and gcd are coprime I can lift lift?(p1:SUPP,p2:SUPP,uterm:UTerm,ldeg:List NNI, _ lvar:List OV) : _ Union(s:SUPP,failed:"failed",notCoprime:"notCoprime") == - leadpol:Boolean:=false - (listpol,lval):=(uterm.lpol,uterm.lint.first) - d:=listpol.first - listpol:=listpol.rest - nolift:Boolean:=true - for uf in listpol repeat - --note uf and d not necessarily primitive - degree gcd(uf,d) =0 => nolift:=false - nolift => ["notCoprime"] - f:SUPP:=([p1,p2]$List(SUPP)).(position(uf,listpol)) - lgcd:=gcd(leadingCoefficient p1, leadingCoefficient p2) - (l:=lift(f,d,uf,lgcd,lvar,ldeg,lval)) case "failed" => ["failed"] + (listpol, lval) := (uterm.lpol, first(uterm.lint)) + d := first(listpol) + listpol := rest(listpol) + uf := listpol(1) + f := p1 + --note uf and d not necessarily primitive + if degree gcd(uf, d) ~= 0 then + uf := listpol(2) + f := p2 + if degree gcd(uf, d) ~= 0 then return ["notCoprime"] + lgcd := gcd(leadingCoefficient p1, leadingCoefficient p2) + l := lift(f, d, uf, lgcd, lvar, ldeg, lval) + l case "failed" => ["failed"] [l :: SUPP] -- interface with the general "lifting" function lift(f:SUPP,d:SUP,uf:SUP,lgcd:P,lvar:List OV, ldeg:List NNI,lval:List R):Union(SUPP,"failed") == - leadpol:Boolean:=false - lcf:P - lcf:=leadingCoefficient f - df:=degree f - leadlist:List(P):=[] - - if lgcd^=1 then - leadpol:=true - f:=lgcd*f - ldeg:=[n0+n1 for n0 in ldeg for n1 in degree(lgcd,lvar)] - lcd:R:=leadingCoefficient d - if degree(lgcd)=0 then d:=((retract lgcd) *d exquo lcd)::SUP - else d:=(retract(eval(lgcd,lvar,lval)) * d exquo lcd)::SUP - uf:=lcd*uf - leadlist:=[lgcd,lcf] + leadpol : Boolean := false + lcf : P + lcf := leadingCoefficient f + df := degree f + leadlist : List(P) := [] + + if lgcd ^= 1 then + leadpol := true + f := lgcd*f + ldeg := [n0+n1 for n0 in ldeg for n1 in degree(lgcd, lvar)] + lcd : R := leadingCoefficient d + lgcd1 := + degree(lgcd) = 0 => retract lgcd + retract(eval(lgcd, lvar, lval)) + du := (lgcd1*d) exquo lcd + du case "failed" => "failed" + d := du::SUP + uf := lcd*uf + leadlist := [lgcd, lcf] lgu := imposelc([d, uf], lvar, lval, leadlist) lgu case "failed" => "failed" lg := lgu::List(SUP) - (pl:=lifting(f,lvar,lg,lval,leadlist,ldeg,pmod)) case "failed" => + (pl := lifting(f,lvar,lg,lval,leadlist,ldeg,pmod)) case "failed" => "failed" plist := pl :: List SUPP - (p0:SUPP,p1:SUPP):=(plist.first,plist.2) - if completeEval(p0,lvar,lval) ^= lg.first then - (p0,p1):=(p1,p0) - ^leadpol => p0 + (p0 : SUPP, p1 : SUPP) := (plist.first, plist.2) + if completeEval(p0, lvar, lval) ^= lg.first then + (p0, p1) := (p1, p0) + not leadpol => p0 p0 exquo content(p0) -- Gcd for two multivariate polynomials @@ -139688,7 +139705,7 @@ See Volume 10.1 for more details. -- Gcd for a list of multivariate polynomials gcd(listp:List P) : P == - lf:=sort(better,listp) + lf := best_to_front(listp) f:=lf.first for g in lf.rest repeat f:=gcd(f,g) diff --git a/changelog b/changelog index e6dcdc3..1f06820 100644 --- a/changelog +++ b/changelog @@ -1,3 +1,7 @@ +20140603 tpd src/axiom-website/patches.html 20140603.03.tpd.patch +20140603 tpd books/bookvol10.4 fix "bad reduction" in multivariate poly gcd +20140603 tpd src/input/pgcd.input add pgcd test case +20140603 tpd src/input/Makefile add pgcd.input 20140603 tpd src/axiom-website/patches.html 20140603.02.tpd.patch 20140603 tpd books/bookvolbib add published references from/to Axiom 20140603 tpd src/axiom-website/patches.html 20140603.01.tpd.patch diff --git a/src/axiom-website/patches.html b/src/axiom-website/patches.html index be17339..c97d73f 100644 --- a/src/axiom-website/patches.html +++ b/src/axiom-website/patches.html @@ -4366,6 +4366,8 @@ books/bookvol10.4 apply Waldek's changes to imposelc in PGCD buglist, add todo 335: add packages to )d op gcd 20140603.02.tpd.patch books/bookvolbib add published references from/to Axiom +20140603.03.tpd.patch +books/bookvol10.4 fix "bad reduction" in multivariate poly gcd diff --git a/src/input/Makefile.pamphlet b/src/input/Makefile.pamphlet index 98205ac..573340d 100644 --- a/src/input/Makefile.pamphlet +++ b/src/input/Makefile.pamphlet @@ -381,7 +381,8 @@ REGRESSTESTS= ackermann.regress \ parabola.regress pascal1.regress pascal.regress \ patch51.regress \ patmatch.regress pat.regress perman.regress perm.regress \ - pfaffian.regress pfr1.regress pfr.regress pmint.regress \ + pfaffian.regress pfr1.regress pfr.regress pgcd.regress \ + pmint.regress \ poly1.regress polycoer.regress poly.regress psgenfcn.regress \ quat1.regress quat.regress r20abugs.regress r20bugs.regress \ r21bugsbig.regress r21bugs.regress radff.regress radix.regress \ @@ -817,7 +818,7 @@ FILES= ${OUT}/ackermann.input \ ${OUT}/patch51.input \ ${OUT}/patmatch.input ${OUT}/perman.input \ ${OUT}/perm.input ${OUT}/pfaffian.input \ - ${OUT}/pfr.input ${OUT}/pfr1.input \ + ${OUT}/pfr.input ${OUT}/pfr1.input ${OUT}/pgcd.input \ ${OUT}/pinch.input ${OUT}/plotfile.input ${OUT}/pollevel.input \ ${OUT}/pmint.input ${OUT}/polycoer.input \ ${OUT}/poly1.input ${OUT}/psgenfcn.input \ @@ -1257,7 +1258,7 @@ DOCFILES= \ ${DOC}/patmatch.input.dvi \ ${DOC}/pdecomp0.as.dvi ${DOC}/perman.input.dvi \ ${DOC}/perm.input.dvi ${DOC}/pfaffian.input.dvi \ - ${DOC}/pfr1.input.dvi \ + ${DOC}/pfr1.input.dvi ${DOC}/pgcd.input.dvi \ ${DOC}/pfr.input.dvi ${DOC}/pinch.input.dvi \ ${DOC}/plotfile.input.dvi ${DOC}/plotlist.input.dvi \ ${DOC}/pollevel.input.dvi ${DOC}/poly1.input.dvi \ diff --git a/src/input/pgcd.input.pamphlet b/src/input/pgcd.input.pamphlet new file mode 100644 index 0000000..dfde7d6 --- /dev/null +++ b/src/input/pgcd.input.pamphlet @@ -0,0 +1,149 @@ +\documentclass{article} +\usepackage{axiom} +\setlength{\textwidth}{400pt} +\begin{document} +\title{\$SPAD/src/input pgcd.input} +\author{Timothy Daly} +\maketitle +\begin{abstract} +This is the test case for polynomial GCDs +\end{abstract} +\eject +\tableofcontents +\eject +\begin{chunk}{*} +)set break resume +)sys rm -f pgcd.output +)spool pgcd.output +)set message test on +)set message auto off +)clear all + +--S 1 of 5 +p1:=x^6 * (50*y^8 - 68*y^6 - 104*y^4 + 36*y^2 +22) + _ + x^4 * (12*y^8 - 10*y^6 - 64*y^4 - 6*y^2 + 4) + _ + x^2 * ( 2*y^6 - 14*y^4 - 8*y^2) + _ + ( -y^4 - y^2) +--R +--R +--R (1) +--R 6 4 8 6 4 2 6 6 4 2 4 +--R (50x + 12x )y + (- 68x - 10x + 2x )y + (- 104x - 64x - 14x - 1)y +--R + +--R 6 4 2 2 6 4 +--R (36x - 6x - 8x - 1)y + 22x + 4x +--R Type: Polynomial(Integer) +--E 1 + +--S 2 of 5 +p2:=x^8 * ( 16*y^10 + 16*y^8 - 96*y^6 - 224*y^4 - 176*y^2 - 48) + _ + x^6 * (-72*y^10 + 188*y^8 - 188*y^6 - 572*y^4 - 108*y^2 + 16) + _ + x^4 * (-64*y^10 + 188*y^8 - 372*y^6 - 60*y^4 - 124*y^2) + _ + x^2 * ( 68*y^8 + 172*y^6 + 97*y^4 - 15*y^2) + _ + ( 8*y^8 + 22*y^6 + 14*y^4) +--R +--R +--R (2) +--R 8 6 4 10 8 6 4 2 8 +--R (16x - 72x - 64x )y + (16x + 188x + 188x + 68x + 8)y +--R + +--R 8 6 4 2 6 +--R (- 96x - 188x - 372x + 172x + 22)y +--R + +--R 8 6 4 2 4 8 6 4 2 2 +--R (- 224x - 572x - 60x + 97x + 14)y + (- 176x - 108x - 124x - 15x )y +--R + +--R 8 6 +--R - 48x + 16x +--R Type: Polynomial(Integer) +--E 2 + +--S 3 of 5 +p1u:=univariate(p1,x) +--R +--R +--R (3) +--R 8 6 4 2 6 8 6 4 2 4 +--R (50y - 68y - 104y + 36y + 22)? + (12y - 10y - 64y - 6y + 4)? +--R + +--R 6 4 2 2 4 2 +--R (2y - 14y - 8y )? - y - y +--R Type: SparseUnivariatePolynomial(Polynomial(Integer)) +--E 3 + +--S 4 of 5 +p2u:=univariate(p2,x) +--R +--R +--R (4) +--R 10 8 6 4 2 8 +--R (16y + 16y - 96y - 224y - 176y - 48)? +--R + +--R 10 8 6 4 2 6 +--R (- 72y + 188y - 188y - 572y - 108y + 16)? +--R + +--R 10 8 6 4 2 4 8 6 4 2 2 +--R (- 64y + 188y - 372y - 60y - 124y )? + (68y + 172y + 97y - 15y )? +--R + +--R 8 6 4 +--R 8y + 22y + 14y +--R Type: SparseUnivariatePolynomial(Polynomial(Integer)) +--E 4 + +--S 5 of 5 +lg:=[gcd(p1u,p2u) for i in 1..1000] +--R +--R +--R (5) +--R [1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, +--R 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1] +--R Type: List(SparseUnivariatePolynomial(Polynomial(Integer))) +--E 5 + +)spool +)lisp (bye) + +\end{chunk} +\eject +\begin{thebibliography}{99} +\bibitem{1} nothing +\end{thebibliography} +\end{document}