\documentclass{article} \usepackage{open-axiom} \begin{document} \title{\$SPAD/src/algebra prs.spad} \author{Lionel Ducos} \maketitle \begin{abstract} \end{abstract} \eject \tableofcontents \eject \section{package PRS PseudoRemainderSequence} The package constructor {\bf PseudoRemainderSequence} provides efficient algorithms by Lionel Ducos (University of Poitiers, France) for computing sub-resultants. This leads to a speed up in many places in Axiom where sub-resultants are computed (polynomial system solving, algebraic factorization, integration). <>= )abbrev package PRS PseudoRemainderSequence ++ Author: Ducos Lionel ++ Date Created: january 1995 ++ Date Last Updated: 5 february 1999 ++ Basic Functions: ++ Related Constructors: ++ Also See: ++ AMS Classifications: ++ Keywords: ++ Description: This package contains some functions: ++ \axiomOpFrom{discriminant}{PseudoRemainderSequence}, ++ \axiomOpFrom{resultant}{PseudoRemainderSequence}, ++ \axiomOpFrom{subResultantGcd}{PseudoRemainderSequence}, ++ \axiomOpFrom{chainSubResultants}{PseudoRemainderSequence}, ++ \axiomOpFrom{degreeSubResultant}{PseudoRemainderSequence}, ++ \axiomOpFrom{lastSubResultant}{PseudoRemainderSequence}, ++ \axiomOpFrom{resultantEuclidean}{PseudoRemainderSequence}, ++ \axiomOpFrom{subResultantGcdEuclidean}{PseudoRemainderSequence}, ++ \axiomOpFrom{semiSubResultantGcdEuclidean1}{PseudoRemainderSequence}, ++ \axiomOpFrom{semiSubResultantGcdEuclidean2}{PseudoRemainderSequence}, etc. ++ This procedures are coming from improvements ++ of the subresultants algorithm. ++ Version : 7 ++ References : Lionel Ducos "Optimizations of the subresultant algorithm" ++ to appear in the Journal of Pure and Applied Algebra. ++ Author : Ducos Lionel \axiom{Lionel.Ducos@mathlabo.univ-poitiers.fr} PseudoRemainderSequence(R, polR) : Specification == Implementation where R : IntegralDomain polR : UnivariatePolynomialCategory(R) NNI ==> NonNegativeInteger LC ==> leadingCoefficient Specification == with resultant : (polR, polR) -> R ++ \axiom{resultant(P, Q)} returns the resultant ++ of \axiom{P} and \axiom{Q} resultantEuclidean : (polR, polR) -> Record(coef1 : polR, coef2 : polR, resultant : R) ++ \axiom{resultantEuclidean(P,Q)} carries out the equality ++ \axiom{coef1*P + coef2*Q = resultant(P,Q)} semiResultantEuclidean2 : (polR, polR) -> Record(coef2 : polR, resultant : R) ++ \axiom{semiResultantEuclidean2(P,Q)} carries out the equality ++ \axiom{...P + coef2*Q = resultant(P,Q)}. ++ Warning: \axiom{degree(P) >= degree(Q)}. semiResultantEuclidean1 : (polR, polR) -> Record(coef1 : polR, resultant : R) ++ \axiom{semiResultantEuclidean1(P,Q)} carries out the equality ++ \axiom{coef1.P + ? Q = resultant(P,Q)}. indiceSubResultant : (polR, polR, NNI) -> polR ++ \axiom{indiceSubResultant(P, Q, i)} returns ++ the subresultant of indice \axiom{i} indiceSubResultantEuclidean : (polR, polR, NNI) -> Record(coef1 : polR, coef2 : polR, subResultant : polR) ++ \axiom{indiceSubResultant(P, Q, i)} returns ++ the subresultant \axiom{S_i(P,Q)} and carries out the equality ++ \axiom{coef1*P + coef2*Q = S_i(P,Q)} semiIndiceSubResultantEuclidean : (polR, polR, NNI) -> Record(coef2 : polR, subResultant : polR) ++ \axiom{semiIndiceSubResultantEuclidean(P, Q, i)} returns ++ the subresultant \axiom{S_i(P,Q)} and carries out the equality ++ \axiom{...P + coef2*Q = S_i(P,Q)} ++ Warning: \axiom{degree(P) >= degree(Q)}. degreeSubResultant : (polR, polR, NNI) -> polR ++ \axiom{degreeSubResultant(P, Q, d)} computes ++ a subresultant of degree \axiom{d}. degreeSubResultantEuclidean : (polR, polR, NNI) -> Record(coef1 : polR, coef2 : polR, subResultant : polR) ++ \axiom{indiceSubResultant(P, Q, i)} returns ++ a subresultant \axiom{S} of degree \axiom{d} ++ and carries out the equality \axiom{coef1*P + coef2*Q = S_i}. semiDegreeSubResultantEuclidean : (polR, polR, NNI) -> Record(coef2 : polR, subResultant : polR) ++ \axiom{indiceSubResultant(P, Q, i)} returns ++ a subresultant \axiom{S} of degree \axiom{d} ++ and carries out the equality \axiom{...P + coef2*Q = S_i}. ++ Warning: \axiom{degree(P) >= degree(Q)}. lastSubResultant : (polR, polR) -> polR ++ \axiom{lastSubResultant(P, Q)} computes ++ the last non zero subresultant of \axiom{P} and \axiom{Q} lastSubResultantEuclidean : (polR, polR) -> Record(coef1 : polR, coef2 : polR, subResultant : polR) ++ \axiom{lastSubResultantEuclidean(P, Q)} computes ++ the last non zero subresultant \axiom{S} ++ and carries out the equality \axiom{coef1*P + coef2*Q = S}. semiLastSubResultantEuclidean : (polR, polR) -> Record(coef2 : polR, subResultant : polR) ++ \axiom{semiLastSubResultantEuclidean(P, Q)} computes ++ the last non zero subresultant \axiom{S} ++ and carries out the equality \axiom{...P + coef2*Q = S}. ++ Warning: \axiom{degree(P) >= degree(Q)}. subResultantGcd : (polR, polR) -> polR ++ \axiom{subResultantGcd(P, Q)} returns the gcd ++ of two primitive polynomials \axiom{P} and \axiom{Q}. subResultantGcdEuclidean : (polR, polR) -> Record(coef1 : polR, coef2 : polR, gcd : polR) ++ \axiom{subResultantGcdEuclidean(P,Q)} carries out the equality ++ \axiom{coef1*P + coef2*Q = +/- S_i(P,Q)} ++ where the degree (not the indice) ++ of the subresultant \axiom{S_i(P,Q)} is the smaller as possible. semiSubResultantGcdEuclidean2 : (polR, polR) -> Record(coef2 : polR, gcd : polR) ++ \axiom{semiSubResultantGcdEuclidean2(P,Q)} carries out the equality ++ \axiom{...P + coef2*Q = +/- S_i(P,Q)} ++ where the degree (not the indice) ++ of the subresultant \axiom{S_i(P,Q)} is the smaller as possible. ++ Warning: \axiom{degree(P) >= degree(Q)}. semiSubResultantGcdEuclidean1: (polR, polR)->Record(coef1: polR, gcd: polR) ++ \axiom{semiSubResultantGcdEuclidean1(P,Q)} carries out the equality ++ \axiom{coef1*P + ? Q = +/- S_i(P,Q)} ++ where the degree (not the indice) ++ of the subresultant \axiom{S_i(P,Q)} is the smaller as possible. discriminant : polR -> R ++ \axiom{discriminant(P, Q)} returns the discriminant ++ of \axiom{P} and \axiom{Q}. discriminantEuclidean : polR -> Record(coef1 : polR, coef2 : polR, discriminant : R) ++ \axiom{discriminantEuclidean(P)} carries out the equality ++ \axiom{coef1 * P + coef2 * D(P) = discriminant(P)}. semiDiscriminantEuclidean : polR -> Record(coef2 : polR, discriminant : R) ++ \axiom{discriminantEuclidean(P)} carries out the equality ++ \axiom{...P + coef2 * D(P) = discriminant(P)}. ++ Warning: \axiom{degree(P) >= degree(Q)}. chainSubResultants : (polR, polR) -> List(polR) ++ \axiom{chainSubResultants(P, Q)} computes the list ++ of non zero subresultants of \axiom{P} and \axiom{Q}. schema : (polR, polR) -> List(NNI) ++ \axiom{schema(P,Q)} returns the list of degrees of ++ non zero subresultants of \axiom{P} and \axiom{Q}. if R has GcdDomain then resultantReduit : (polR, polR) -> R ++ \axiom{resultantReduit(P,Q)} returns the "reduce resultant" ++ of \axiom{P} and \axiom{Q}. resultantReduitEuclidean : (polR, polR) -> Record(coef1 : polR, coef2 : polR, resultantReduit : R) ++ \axiom{resultantReduitEuclidean(P,Q)} returns ++ the "reduce resultant" and carries out the equality ++ \axiom{coef1*P + coef2*Q = resultantReduit(P,Q)}. semiResultantReduitEuclidean : (polR, polR) -> Record(coef2 : polR, resultantReduit : R) ++ \axiom{semiResultantReduitEuclidean(P,Q)} returns ++ the "reduce resultant" and carries out the equality ++ \axiom{...P + coef2*Q = resultantReduit(P,Q)}. gcd : (polR, polR) -> polR ++ \axiom{gcd(P, Q)} returns the gcd of \axiom{P} and \axiom{Q}. -- sub-routines exported for convenience ---------------------------- * : (R, Vector(polR)) -> Vector(polR) ++ \axiom{r * v} computes the product of \axiom{r} and \axiom{v} exquo : (Vector(polR), R) -> Vector(polR) ++ \axiom{v exquo r} computes ++ the exact quotient of \axiom{v} by \axiom{r} pseudoDivide : (polR, polR) -> Record(coef:R, quotient:polR, remainder:polR) ++ \axiom{pseudoDivide(P,Q)} computes the pseudoDivide ++ of \axiom{P} by \axiom{Q}. divide : (polR, polR) -> Record(quotient : polR, remainder : polR) ++ \axiom{divide(F,G)} computes quotient and rest ++ of the exact euclidean division of \axiom{F} by \axiom{G}. Lazard : (R, R, NNI) -> R ++ \axiom{Lazard(x, y, n)} computes \axiom{x**n/y**(n-1)} Lazard2 : (polR, R, R, NNI) -> polR ++ \axiom{Lazard2(F, x, y, n)} computes \axiom{(x/y)**(n-1) * F} next_sousResultant2 : (polR, polR, polR, R) -> polR ++ \axiom{nextsousResultant2(P, Q, Z, s)} returns ++ the subresultant \axiom{S_{e-1}} where ++ \axiom{P ~ S_d, Q = S_{d-1}, Z = S_e, s = lc(S_d)} resultant_naif : (polR, polR) -> R ++ \axiom{resultantEuclidean_naif(P,Q)} returns ++ the resultant of \axiom{P} and \axiom{Q} computed ++ by means of the naive algorithm. resultantEuclidean_naif : (polR, polR) -> Record(coef1 : polR, coef2 : polR, resultant : R) ++ \axiom{resultantEuclidean_naif(P,Q)} returns ++ the extended resultant of \axiom{P} and \axiom{Q} computed ++ by means of the naive algorithm. semiResultantEuclidean_naif : (polR, polR) -> Record(coef2 : polR, resultant : R) ++ \axiom{resultantEuclidean_naif(P,Q)} returns ++ the semi-extended resultant of \axiom{P} and \axiom{Q} computed ++ by means of the naive algorithm. Implementation == add X : polR := monomial(1$R,1) r : R * v : Vector(polR) == r::polR * v -- the instruction map(r * #1, v) is slower !? v : Vector(polR) exquo r : R == map((#1 exquo r)::polR, v) pseudoDivide(P : polR, Q : polR) : Record(coef:R,quotient:polR,remainder:polR) == -- computes the pseudoDivide of P by Q zero?(Q) => error("PseudoDivide$PRS : division by 0") zero?(P) => construct(1, 0, P) lcQ : R := LC(Q) (degP, degQ) := (degree(P), degree(Q)) degP < degQ => construct(1, 0, P) Q := reductum(Q) i : NNI := (degP - degQ + 1)::NNI co : R := lcQ**i quot : polR := 0$polR while (delta : Integer := degree(P) - degQ) >= 0 repeat i := (i - 1)::NNI mon := monomial(LC(P), delta::NNI)$polR quot := quot + lcQ**i * mon P := lcQ * reductum(P) - mon * Q P := lcQ**i * P return construct(co, quot, P) divide(F : polR, G : polR) : Record(quotient : polR, remainder : polR)== -- computes quotient and rest of the exact euclidean division of F by G lcG : R := LC(G) degG : NNI := degree(G) zero?(degG) => ( F := (F exquo lcG)::polR; return construct(F, 0)) G : polR := reductum(G) quot : polR := 0 while (delta := degree(F) - degG) >= 0 repeat mon : polR := monomial((LC(F) exquo lcG)::R, delta::NNI) quot := quot + mon F := reductum(F) - mon * G return construct(quot, F) resultant_naif(P : polR, Q : polR) : R == -- valid over a field a : R := 1 repeat zero?(Q) => return 0 (degP, degQ) := (degree(P), degree(Q)) if odd?(degP) and odd?(degQ) then a := - a zero?(degQ) => return (a * LC(Q)**degP) U : polR := divide(P, Q).remainder a := a * LC(Q)**(degP - degree(U))::NNI (P, Q) := (Q, U) resultantEuclidean_naif(P : polR, Q : polR) : Record(coef1 : polR, coef2 : polR, resultant : R) == -- valid over a field. a : R := 1 old_cf1 : polR := 1 ; cf1 : polR := 0 old_cf2 : polR := 0 ; cf2 : polR := 1 repeat zero?(Q) => return construct(0::polR, 0::polR, 0::R) (degP, degQ) := (degree(P), degree(Q)) if odd?(degP) and odd?(degQ) then a := -a if zero?(degQ) then a := a * LC(Q)**(degP-1)::NNI return construct(a*cf1, a*cf2, a*LC(Q)) divid := divide(P,Q) a := a * LC(Q)**(degP - degree(divid.remainder))::NNI (P, Q) := (Q, divid.remainder) (old_cf1, old_cf2, cf1, cf2) := (cf1, cf2, old_cf1 - divid.quotient * cf1, old_cf2 - divid.quotient * cf2) semiResultantEuclidean_naif(P : polR, Q : polR) : Record(coef2 : polR, resultant : R) == -- valid over a field a : R := 1 old_cf2 : polR := 0 ; cf2 : polR := 1 repeat zero?(Q) => return construct(0::polR, 0::R) (degP, degQ) := (degree(P), degree(Q)) if odd?(degP) and odd?(degQ) then a := -a if zero?(degQ) then a := a * LC(Q)**(degP-1)::NNI return construct(a*cf2, a*LC(Q)) divid := divide(P,Q) a := a * LC(Q)**(degP - degree(divid.remainder))::NNI (P, Q) := (Q, divid.remainder) (old_cf2, cf2) := (cf2, old_cf2 - divid.quotient * cf2) Lazard(x : R, y : R, n : NNI) : R == zero?(n) => error("Lazard$PRS : n = 0") one?(n) => x a : NNI := 1 while n >= (b := 2*a) repeat a := b c : R := x n := (n - a)::NNI repeat -- c = x**i / y**(i-1), i=n_0 quo a, a=2**? one?(a) => return c a := a quo 2 c := ((c * c) exquo y)::R if n >= a then ( c := ((c * x) exquo y)::R ; n := (n - a)::NNI ) Lazard2(F : polR, x : R, y : R, n : NNI) : polR == zero?(n) => error("Lazard2$PRS : n = 0") one?(n) => F x := Lazard(x, y, (n-1)::NNI) return ((x * F) exquo y)::polR Lazard3(V : Vector(polR), x : R, y : R, n : NNI) : Vector(polR) == -- computes x**(n-1) * V / y**(n-1) zero?(n) => error("Lazard2$prs : n = 0") one?(n) => V x := Lazard(x, y, (n-1)::NNI) return ((x * V) exquo y) next_sousResultant2(P : polR, Q : polR, Z : polR, s : R) : polR == (lcP, c, se) := (LC(P), LC(Q), LC(Z)) (d, e) := (degree(P), degree(Q)) (P, Q, H) := (reductum(P), reductum(Q), - reductum(Z)) A : polR := coefficient(P, e) * H for i in e+1..d-1 repeat H := if degree(H) = e-1 then X * reductum(H) - ((LC(H) * Q) exquo c)::polR else X * H -- H = s_e * X^i mod S_d-1 A := coefficient(P, i) * H + A while degree(P) >= e repeat P := reductum(P) A := A + se * P -- A = s_e * reductum(P_0) mod S_d-1 A := (A exquo lcP)::polR -- A = s_e * reductum(S_d) / s_d mod S_d-1 A := if degree(H) = e-1 then c * (X * reductum(H) + A) - LC(H) * Q else c * (X * H + A) A := (A exquo s)::polR -- A = +/- S_e-1 return (if odd?(d-e) then A else - A) next_sousResultant3(VP : Vector(polR), VQ : Vector(polR), s : R, ss : R) : Vector(polR) == -- P ~ S_d, Q = S_d-1, s = lc(S_d), ss = lc(S_e) (P, Q) := (VP.1, VQ.1) (lcP, c) := (LC(P), LC(Q)) e : NNI := degree(Q) if one?(delta := degree(P) - e) then -- algo_new VP := c * VP - coefficient(P, e) * VQ VP := VP exquo lcP VP := c * (VP - X * VQ) + coefficient(Q, (e-1)::NNI) * VQ VP := VP exquo s else -- algorithm of Lickteig - Roy (r, rr) := (s * lcP, ss * c) divid := divide(rr * P, Q) VP.1 := (divid.remainder exquo r)::polR for i in 2..#VP repeat VP.i := rr * VP.i - VQ.i * divid.quotient VP.i := (VP.i exquo r)::polR return (if odd?(delta) then VP else - VP) algo_new(P : polR, Q : polR) : R == delta : NNI := (degree(P) - degree(Q))::NNI s : R := LC(Q)**delta (P, Q) := (Q, pseudoRemainder(P, -Q)) repeat -- P = S_c-1 (except the first turn : P ~ S_c-1), -- Q = S_d-1, s = lc(S_d) zero?(Q) => return 0 delta := (degree(P) - degree(Q))::NNI Z : polR := Lazard2(Q, LC(Q), s, delta) -- Z = S_e ~ S_d-1 zero?(degree(Z)) => return LC(Z) (P, Q) := (Q, next_sousResultant2(P, Q, Z, s)) s := LC(Z) resultant(P : polR, Q : polR) : R == zero?(Q) or zero?(P) => 0 if degree(P) < degree(Q) then (P, Q) := (Q, P) if odd?(degree(P)) and odd?(degree(Q)) then Q := - Q zero?(degree(Q)) => LC(Q)**degree(P) -- degree(P) >= degree(Q) > 0 R has Finite => resultant_naif(P, Q) return algo_new(P, Q) subResultantEuclidean(P : polR, Q : polR) : Record(coef1 : polR, coef2 : polR, resultant : R) == s : R := LC(Q)**(degree(P) - degree(Q))::NNI VP : Vector(polR) := [Q, 0::polR, 1::polR] pdiv := pseudoDivide(P, -Q) VQ : Vector(polR) := [pdiv.remainder, pdiv.coef::polR, pdiv.quotient] repeat -- VP.1 = S_{c-1}, VQ.1 = S_{d-1}, s=lc(S_d) -- S_{c-1} = VP.2 P_0 + VP.3 Q_0, S_{d-1} = VQ.2 P_0 + VQ.3 Q_0 (P, Q) := (VP.1, VQ.1) zero?(Q) => return construct(0::polR, 0::polR, 0::R) e : NNI := degree(Q) delta : NNI := (degree(P) - e)::NNI if zero?(e) then l : Vector(polR) := Lazard3(VQ, LC(Q), s, delta) return construct(l.2, l.3, LC(l.1)) ss : R := Lazard(LC(Q), s, delta) (VP, VQ) := (VQ, next_sousResultant3(VP, VQ, s, ss)) s := ss resultantEuclidean(P : polR, Q : polR) : Record(coef1 : polR, coef2 : polR, resultant : R) == zero?(P) or zero?(Q) => construct(0::polR, 0::polR, 0::R) if degree(P) < degree(Q) then e : Integer := if odd?(degree(P)) and odd?(degree(Q)) then -1 else 1 l := resultantEuclidean(Q, e * P) return construct(e * l.coef2, l.coef1, l.resultant) if zero?(degree(Q)) then degP : NNI := degree(P) zero?(degP) => error("resultantEuclidean$PRS : constant polynomials") s : R := LC(Q)**(degP-1)::NNI return construct(0::polR, s::polR, s * LC(Q)) R has Finite => resultantEuclidean_naif(P, Q) return subResultantEuclidean(P,Q) semiSubResultantEuclidean(P : polR, Q : polR) : Record(coef2 : polR, resultant : R) == s : R := LC(Q)**(degree(P) - degree(Q))::NNI VP : Vector(polR) := [Q, 1::polR] pdiv := pseudoDivide(P, -Q) VQ : Vector(polR) := [pdiv.remainder, pdiv.quotient] repeat -- VP.1 = S_{c-1}, VQ.1 = S_{d-1}, s=lc(S_d) -- S_{c-1} = ...P_0 + VP.3 Q_0, S_{d-1} = ...P_0 + VQ.3 Q_0 (P, Q) := (VP.1, VQ.1) zero?(Q) => return construct(0::polR, 0::R) e : NNI := degree(Q) delta : NNI := (degree(P) - e)::NNI if zero?(e) then l : Vector(polR) := Lazard3(VQ, LC(Q), s, delta) return construct(l.2, LC(l.1)) ss : R := Lazard(LC(Q), s, delta) (VP, VQ) := (VQ, next_sousResultant3(VP, VQ, s, ss)) s := ss semiResultantEuclidean2(P : polR, Q : polR) : Record(coef2 : polR, resultant : R) == zero?(P) or zero?(Q) => construct(0::polR, 0::R) degree(P) < degree(Q) => error("semiResultantEuclidean2 : bad degrees") if zero?(degree(Q)) then degP : NNI := degree(P) zero?(degP) => error("semiResultantEuclidean2 : constant polynomials") s : R := LC(Q)**(degP-1)::NNI return construct(s::polR, s * LC(Q)) R has Finite => semiResultantEuclidean_naif(P, Q) return semiSubResultantEuclidean(P,Q) semiResultantEuclidean1(P : polR, Q : polR) : Record(coef1 : polR, resultant : R) == result := resultantEuclidean(P,Q) [result.coef1, result.resultant] indiceSubResultant(P : polR, Q : polR, i : NNI) : polR == zero?(Q) or zero?(P) => 0 if degree(P) < degree(Q) then (P, Q) := (Q, P) if odd?(degree(P)-i) and odd?(degree(Q)-i) then Q := - Q if i = degree(Q) then delta : NNI := (degree(P)-degree(Q))::NNI zero?(delta) => error("indiceSubResultant$PRS : bad degrees") s : R := LC(Q)**(delta-1)::NNI return s*Q i > degree(Q) => 0 s : R := LC(Q)**(degree(P) - degree(Q))::NNI (P, Q) := (Q, pseudoRemainder(P, -Q)) repeat -- P = S_{c-1} ~ S_d , Q = S_{d-1}, s = lc(S_d), i < d (degP, degQ) := (degree(P), degree(Q)) i = degP-1 => return Q zero?(Q) or (i > degQ) => return 0 Z : polR := Lazard2(Q, LC(Q), s, (degP - degQ)::NNI) -- Z = S_e ~ S_d-1 i = degQ => return Z (P, Q) := (Q, next_sousResultant2(P, Q, Z, s)) s := LC(Z) indiceSubResultantEuclidean(P : polR, Q : polR, i : NNI) : Record(coef1 : polR, coef2 : polR, subResultant : polR) == zero?(Q) or zero?(P) => construct(0::polR, 0::polR, 0::polR) if degree(P) < degree(Q) then e := if odd?(degree(P)-i) and odd?(degree(Q)-i) then -1 else 1 l := indiceSubResultantEuclidean(Q, e * P, i) return construct(e * l.coef2, l.coef1, l.subResultant) if i = degree(Q) then delta : NNI := (degree(P)-degree(Q))::NNI zero?(delta) => error("indiceSubResultantEuclidean$PRS : bad degrees") s : R := LC(Q)**(delta-1)::NNI return construct(0::polR, s::polR, s * Q) i > degree(Q) => construct(0::polR, 0::polR, 0::polR) s : R := LC(Q)**(degree(P) - degree(Q))::NNI VP : Vector(polR) := [Q, 0::polR, 1::polR] pdiv := pseudoDivide(P, -Q) VQ : Vector(polR) := [pdiv.remainder, pdiv.coef::polR, pdiv.quotient] repeat -- VP.1 = S_{c-1}, VQ.1 = S_{d-1}, s=lc(S_d), i < d -- S_{c-1} = VP.2 P_0 + VP.3 Q_0, S_{d-1} = VQ.2 P_0 + VQ.3 Q_0 (P, Q) := (VP.1, VQ.1) zero?(Q) => return construct(0::polR, 0::polR, 0::polR) (degP, degQ) := (degree(P), degree(Q)) i = degP-1 => return construct(VQ.2, VQ.3, VQ.1) (i > degQ) => return construct(0::polR, 0::polR, 0::polR) VZ := Lazard3(VQ, LC(Q), s, (degP - degQ)::NNI) i = degQ => return construct(VZ.2, VZ.3, VZ.1) ss : R := LC(VZ.1) (VP, VQ) := (VQ, next_sousResultant3(VP, VQ, s, ss)) s := ss semiIndiceSubResultantEuclidean(P : polR, Q : polR, i : NNI) : Record(coef2 : polR, subResultant : polR) == zero?(Q) or zero?(P) => construct(0::polR, 0::polR) degree(P) < degree(Q) => error("semiIndiceSubResultantEuclidean$PRS : bad degrees") if i = degree(Q) then delta : NNI := (degree(P)-degree(Q))::NNI zero?(delta) => error("semiIndiceSubResultantEuclidean$PRS : bad degrees") s : R := LC(Q)**(delta-1)::NNI return construct(s::polR, s * Q) i > degree(Q) => construct(0::polR, 0::polR) s : R := LC(Q)**(degree(P) - degree(Q))::NNI VP : Vector(polR) := [Q, 1::polR] pdiv := pseudoDivide(P, -Q) VQ : Vector(polR) := [pdiv.remainder, pdiv.quotient] repeat -- VP.1 = S_{c-1}, VQ.1 = S_{d-1}, s = lc(S_d), i < d -- S_{c-1} = ...P_0 + VP.2 Q_0, S_{d-1} = ...P_0 + ...Q_0 (P, Q) := (VP.1, VQ.1) zero?(Q) => return construct(0::polR, 0::polR) (degP, degQ) := (degree(P), degree(Q)) i = degP-1 => return construct(VQ.2, VQ.1) (i > degQ) => return construct(0::polR, 0::polR) VZ := Lazard3(VQ, LC(Q), s, (degP - degQ)::NNI) i = degQ => return construct(VZ.2, VZ.1) ss : R := LC(VZ.1) (VP, VQ) := (VQ, next_sousResultant3(VP, VQ, s, ss)) s := ss degreeSubResultant(P : polR, Q : polR, i : NNI) : polR == zero?(Q) or zero?(P) => 0 if degree(P) < degree(Q) then (P, Q) := (Q, P) if i = degree(Q) then delta : NNI := (degree(P)-degree(Q))::NNI zero?(delta) => error("degreeSubResultant$PRS : bad degrees") s : R := LC(Q)**(delta-1)::NNI return s*Q i > degree(Q) => 0 s : R := LC(Q)**(degree(P) - degree(Q))::NNI (P, Q) := (Q, pseudoRemainder(P, -Q)) repeat -- P = S_{c-1}, Q = S_{d-1}, s = lc(S_d) zero?(Q) or (i > degree(Q)) => return 0 i = degree(Q) => return Q Z : polR := Lazard2(Q, LC(Q), s, (degree(P) - degree(Q))::NNI) -- Z = S_e ~ S_d-1 (P, Q) := (Q, next_sousResultant2(P, Q, Z, s)) s := LC(Z) degreeSubResultantEuclidean(P : polR, Q : polR, i : NNI) : Record(coef1 : polR, coef2 : polR, subResultant : polR) == zero?(Q) or zero?(P) => construct(0::polR, 0::polR, 0::polR) if degree(P) < degree(Q) then l := degreeSubResultantEuclidean(Q, P, i) return construct(l.coef2, l.coef1, l.subResultant) if i = degree(Q) then delta : NNI := (degree(P)-degree(Q))::NNI zero?(delta) => error("degreeSubResultantEuclidean$PRS : bad degrees") s : R := LC(Q)**(delta-1)::NNI return construct(0::polR, s::polR, s * Q) i > degree(Q) => construct(0::polR, 0::polR, 0::polR) s : R := LC(Q)**(degree(P) - degree(Q))::NNI VP : Vector(polR) := [Q, 0::polR, 1::polR] pdiv := pseudoDivide(P, -Q) VQ : Vector(polR) := [pdiv.remainder, pdiv.coef::polR, pdiv.quotient] repeat -- VP.1 = S_{c-1}, VQ.1 = S_{d-1}, s=lc(S_d) -- S_{c-1} = ...P_0 + VP.3 Q_0, S_{d-1} = ...P_0 + VQ.3 Q_0 (P, Q) := (VP.1, VQ.1) zero?(Q) or (i > degree(Q)) => return construct(0::polR, 0::polR, 0::polR) i = degree(Q) => return construct(VQ.2, VQ.3, VQ.1) ss : R := Lazard(LC(Q), s, (degree(P)-degree(Q))::NNI) (VP, VQ) := (VQ, next_sousResultant3(VP, VQ, s, ss)) s := ss semiDegreeSubResultantEuclidean(P : polR, Q : polR, i : NNI) : Record(coef2 : polR, subResultant : polR) == zero?(Q) or zero?(P) => construct(0::polR, 0::polR) degree(P) < degree(Q) => error("semiDegreeSubResultantEuclidean$PRS : bad degrees") if i = degree(Q) then delta : NNI := (degree(P)-degree(Q))::NNI zero?(delta) => error("semiDegreeSubResultantEuclidean$PRS : bad degrees") s : R := LC(Q)**(delta-1)::NNI return construct(s::polR, s * Q) i > degree(Q) => construct(0::polR, 0::polR) s : R := LC(Q)**(degree(P) - degree(Q))::NNI VP : Vector(polR) := [Q, 1::polR] pdiv := pseudoDivide(P, -Q) VQ : Vector(polR) := [pdiv.remainder, pdiv.quotient] repeat -- VP.1 = S_{c-1}, VQ.1 = S_{d-1}, s=lc(S_d) -- S_{c-1} = ...P_0 + VP.3 Q_0, S_{d-1} = ...P_0 + VQ.3 Q_0 (P, Q) := (VP.1, VQ.1) zero?(Q) or (i > degree(Q)) => return construct(0::polR, 0::polR) i = degree(Q) => return construct(VQ.2, VQ.1) ss : R := Lazard(LC(Q), s, (degree(P)-degree(Q))::NNI) (VP, VQ) := (VQ, next_sousResultant3(VP, VQ, s, ss)) s := ss lastSubResultant(P : polR, Q : polR) : polR == zero?(Q) or zero?(P) => 0 if degree(P) < degree(Q) then (P, Q) := (Q, P) zero?(degree(Q)) => (LC(Q)**degree(P))::polR s : R := LC(Q)**(degree(P) - degree(Q))::NNI (P, Q) := (Q, pseudoRemainder(P, -Q)) Z : polR := P repeat -- Z = S_d (except the first turn : Z = P) -- P = S_{c-1} ~ S_d, Q = S_{d-1}, s = lc(S_d) zero?(Q) => return Z Z := Lazard2(Q, LC(Q), s, (degree(P) - degree(Q))::NNI) -- Z = S_e ~ S_{d-1} zero?(degree(Z)) => return Z (P, Q) := (Q, next_sousResultant2(P, Q, Z, s)) s := LC(Z) lastSubResultantEuclidean(P : polR, Q : polR) : Record(coef1 : polR, coef2 : polR, subResultant : polR) == zero?(Q) or zero?(P) => construct(0::polR, 0::polR, 0::polR) if degree(P) < degree(Q) then l := lastSubResultantEuclidean(Q, P) return construct(l.coef2, l.coef1, l.subResultant) if zero?(degree(Q)) then degP : NNI := degree(P) zero?(degP) => error("lastSubResultantEuclidean$PRS : constant polynomials") s : R := LC(Q)**(degP-1)::NNI return construct(0::polR, s::polR, s * Q) s : R := LC(Q)**(degree(P) - degree(Q))::NNI VP : Vector(polR) := [Q, 0::polR, 1::polR] pdiv := pseudoDivide(P, -Q) VQ : Vector(polR) := [pdiv.remainder, pdiv.coef::polR, pdiv.quotient] VZ : Vector(polR) := copy(VP) repeat -- VZ.1 = S_d, VP.1 = S_{c-1}, VQ.1 = S_{d-1}, s = lc(S_d) -- S_{c-1} = VP.2 P_0 + VP.3 Q_0 -- S_{d-1} = VQ.2 P_0 + VQ.3 Q_0 -- S_d = VZ.2 P_0 + VZ.3 Q_0 (Q, Z) := (VQ.1, VZ.1) zero?(Q) => return construct(VZ.2, VZ.3, VZ.1) VZ := Lazard3(VQ, LC(Q), s, (degree(Z) - degree(Q))::NNI) zero?(degree(Q)) => return construct(VZ.2, VZ.3, VZ.1) ss : R := LC(VZ.1) (VP, VQ) := (VQ, next_sousResultant3(VP, VQ, s, ss)) s := ss semiLastSubResultantEuclidean(P : polR, Q : polR) : Record(coef2 : polR, subResultant : polR) == zero?(Q) or zero?(P) => construct(0::polR, 0::polR) degree(P) < degree(Q) => error("semiLastSubResultantEuclidean$PRS : bad degrees") if zero?(degree(Q)) then degP : NNI := degree(P) zero?(degP) => error("semiLastSubResultantEuclidean$PRS : constant polynomials") s : R := LC(Q)**(degP-1)::NNI return construct(s::polR, s * Q) s : R := LC(Q)**(degree(P) - degree(Q))::NNI VP : Vector(polR) := [Q, 1::polR] pdiv := pseudoDivide(P, -Q) VQ : Vector(polR) := [pdiv.remainder, pdiv.quotient] VZ : Vector(polR) := copy(VP) repeat -- VZ.1 = S_d, VP.1 = S_{c-1}, VQ.1 = S_{d-1}, s = lc(S_d) -- S_{c-1} = ... P_0 + VP.2 Q_0 -- S_{d-1} = ... P_0 + VQ.2 Q_0 -- S_d = ... P_0 + VZ.2 Q_0 (Q, Z) := (VQ.1, VZ.1) zero?(Q) => return construct(VZ.2, VZ.1) VZ := Lazard3(VQ, LC(Q), s, (degree(Z) - degree(Q))::NNI) zero?(degree(Q)) => return construct(VZ.2, VZ.1) ss : R := LC(VZ.1) (VP, VQ) := (VQ, next_sousResultant3(VP, VQ, s, ss)) s := ss chainSubResultants(P : polR, Q : polR) : List(polR) == zero?(Q) or zero?(P) => [] if degree(P) < degree(Q) then (P, Q) := (Q, P) if odd?(degree(P)) and odd?(degree(Q)) then Q := - Q L : List(polR) := [] zero?(degree(Q)) => L L := [Q] s : R := LC(Q)**(degree(P) - degree(Q))::NNI (P, Q) := (Q, pseudoRemainder(P, -Q)) repeat -- P = S_{c-1}, Q = S_{d-1}, s = lc(S_d) -- L = [S_d,....,S_{q-1}] zero?(Q) => return L L := concat(Q, L) -- L = [S_{d-1},....,S_{q-1}] delta : NNI := (degree(P) - degree(Q))::NNI Z : polR := Lazard2(Q, LC(Q), s, delta) -- Z = S_e ~ S_d-1 if delta > 1 then L := concat(Z, L) -- L = [S_e,....,S_{q-1}] zero?(degree(Z)) => return L (P, Q) := (Q, next_sousResultant2(P, Q, Z, s)) s := LC(Z) schema(P : polR, Q : polR) : List(NNI) == zero?(Q) or zero?(P) => [] if degree(P) < degree(Q) then (P, Q) := (Q, P) zero?(degree(Q)) => [0] L : List(NNI) := [] s : R := LC(Q)**(degree(P) - degree(Q))::NNI (P, Q) := (Q, pseudoRemainder(P, Q)) repeat -- P = S_{c-1} ~ S_d, Q = S_{d-1}, s = lc(S_d) zero?(Q) => return L e : NNI := degree(Q) L := concat(e, L) delta : NNI := (degree(P) - e)::NNI Z : polR := Lazard2(Q, LC(Q), s, delta) -- Z = S_e ~ S_d-1 if delta > 1 then L := concat(e, L) zero?(e) => return L (P, Q) := (Q, next_sousResultant2(P, Q, Z, s)) s := LC(Z) subResultantGcd(P : polR, Q : polR) : polR == zero?(P) and zero?(Q) => 0 zero?(P) => Q zero?(Q) => P if degree(P) < degree(Q) then (P, Q) := (Q, P) zero?(degree(Q)) => 1$polR s : R := LC(Q)**(degree(P) - degree(Q))::NNI (P, Q) := (Q, pseudoRemainder(P, -Q)) repeat -- P = S_{c-1}, Q = S_{d-1}, s = lc(S_d) zero?(Q) => return P zero?(degree(Q)) => return 1$polR Z : polR := Lazard2(Q, LC(Q), s, (degree(P) - degree(Q))::NNI) -- Z = S_e ~ S_d-1 (P, Q) := (Q, next_sousResultant2(P, Q, Z, s)) s := LC(Z) subResultantGcdEuclidean(P : polR, Q : polR) : Record(coef1 : polR, coef2 : polR, gcd : polR) == zero?(P) and zero?(Q) => construct(0::polR, 0::polR, 0::polR) zero?(P) => construct(0::polR, 1::polR, Q) zero?(Q) => construct(1::polR, 0::polR, P) if degree(P) < degree(Q) then l := subResultantGcdEuclidean(Q, P) return construct(l.coef2, l.coef1, l.gcd) zero?(degree(Q)) => construct(0::polR, 1::polR, Q) s : R := LC(Q)**(degree(P) - degree(Q))::NNI VP : Vector(polR) := [Q, 0::polR, 1::polR] pdiv := pseudoDivide(P, -Q) VQ : Vector(polR) := [pdiv.remainder, pdiv.coef::polR, pdiv.quotient] repeat -- VP.1 = S_{c-1}, VQ.1 = S_{d-1}, s=lc(S_d) -- S_{c-1} = VP.2 P_0 + VP.3 Q_0, S_{d-1} = VQ.2 P_0 + VQ.3 Q_0 (P, Q) := (VP.1, VQ.1) zero?(Q) => return construct(VP.2, VP.3, P) e : NNI := degree(Q) zero?(e) => return construct(VQ.2, VQ.3, Q) ss := Lazard(LC(Q), s, (degree(P) - e)::NNI) (VP,VQ) := (VQ, next_sousResultant3(VP, VQ, s, ss)) s := ss semiSubResultantGcdEuclidean2(P : polR, Q : polR) : Record(coef2 : polR, gcd : polR) == zero?(P) and zero?(Q) => construct(0::polR, 0::polR) zero?(P) => construct(1::polR, Q) zero?(Q) => construct(0::polR, P) degree(P) < degree(Q) => error("semiSubResultantGcdEuclidean2$PRS : bad degrees") zero?(degree(Q)) => construct(1::polR, Q) s : R := LC(Q)**(degree(P) - degree(Q))::NNI VP : Vector(polR) := [Q, 1::polR] pdiv := pseudoDivide(P, -Q) VQ : Vector(polR) := [pdiv.remainder, pdiv.quotient] repeat -- P=S_{c-1}, Q=S_{d-1}, s=lc(S_d) -- S_{c-1} = ? P_0 + old_cf2 Q_0, S_{d-1} = ? P_0 + cf2 Q_0 (P, Q) := (VP.1, VQ.1) zero?(Q) => return construct(VP.2, P) e : NNI := degree(Q) zero?(e) => return construct(VQ.2, Q) ss := Lazard(LC(Q), s, (degree(P) - e)::NNI) (VP,VQ) := (VQ, next_sousResultant3(VP, VQ, s, ss)) s := ss semiSubResultantGcdEuclidean1(P : polR, Q : polR) : Record(coef1 : polR, gcd : polR) == result := subResultantGcdEuclidean(P,Q) [result.coef1, result.gcd] discriminant(P : polR) : R == d : Integer := degree(P) zero?(d) => error "cannot take discriminant of constants" a : Integer := (d * (d-1)) quo 2 a := (-1)**a::NonNegativeInteger dP : polR := differentiate P r : R := resultant(P, dP) d := d - degree(dP) - 1 return (if zero?(d) then a * (r exquo LC(P))::R else a * r * LC(P)**(d-1)::NNI) discriminantEuclidean(P : polR) : Record(coef1 : polR, coef2 : polR, discriminant : R) == d : Integer := degree(P) zero?(d) => error "cannot take discriminant of constants" a : Integer := (d * (d-1)) quo 2 a := (-1)**a::NonNegativeInteger dP : polR := differentiate P rE := resultantEuclidean(P, dP) d := d - degree(dP) - 1 zero? d => [a * (rE.coef1 exquo LC(P))::polR, a * (rE.coef2 exquo LC(P))::polR, a * (rE.resultant exquo LC(P))::R] [a * rE.coef1 * LC(P)**(d-1)::NNI, a * rE.coef2 * LC(P)**(d-1)::NNI, a * rE.resultant * LC(P)**(d-1)::NNI] semiDiscriminantEuclidean(P : polR) : Record(coef2 : polR, discriminant : R) == d : Integer := degree(P) zero?(d) => error "cannot take discriminant of constants" a : Integer := (d * (d-1)) quo 2 a := (-1)**a::NonNegativeInteger dP : polR := differentiate P rE := semiResultantEuclidean2(P, dP) d := d - degree(dP) - 1 zero? d => [a * (rE.coef2 exquo LC(P))::polR, a * (rE.resultant exquo LC(P))::R] [a * rE.coef2 * LC(P)**(d-1)::NNI, a * rE.resultant * LC(P)**(d-1)::NNI] if R has GcdDomain then resultantReduit(P : polR, Q : polR) : R == UV := subResultantGcdEuclidean(P, Q) UVs : polR := UV.gcd positive? degree(UVs) => 0 l : List(R) := concat(coefficients(UV.coef1), coefficients(UV.coef2)) return (LC(UVs) exquo gcd(l))::R resultantReduitEuclidean(P : polR, Q : polR) : Record(coef1 : polR, coef2 : polR, resultantReduit : R) == UV := subResultantGcdEuclidean(P, Q) UVs : polR := UV.gcd positive? degree(UVs) => construct(0::polR, 0::polR, 0::R) l : List(R) := concat(coefficients(UV.coef1), coefficients(UV.coef2)) gl : R := gcd(l) c1 : polR := (UV.coef1 exquo gl)::polR c2 : polR := (UV.coef2 exquo gl)::polR rr : R := (LC(UVs) exquo gl)::R return construct(c1, c2, rr) semiResultantReduitEuclidean(P : polR, Q : polR) : Record(coef2 : polR, resultantReduit : R) == UV := subResultantGcdEuclidean(P, Q) UVs : polR := UV.gcd positive? degree(UVs) => construct(0::polR, 0::R) l : List(R) := concat(coefficients(UV.coef1), coefficients(UV.coef2)) gl : R := gcd(l) c2 : polR := (UV.coef2 exquo gl)::polR rr : R := (LC(UVs) exquo gl)::R return construct(c2, rr) gcd_naif(P : polR, Q : polR) : polR == -- valid over a field zero?(P) => (Q exquo LC(Q))::polR repeat zero?(Q) => return (P exquo LC(P))::polR zero?(degree(Q)) => return 1$polR (P, Q) := (Q, divide(P, Q).remainder) gcd(P : polR, Q : polR) : polR == R has Finite => gcd_naif(P,Q) zero?(P) => Q zero?(Q) => P cP : R := content(P) cQ : R := content(Q) P := (P exquo cP)::polR Q := (Q exquo cQ)::polR G : polR := subResultantGcd(P, Q) return gcd(cP,cQ) * primitivePart(G) @ \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}