\documentclass{article} \usepackage{open-axiom} \begin{document} \title{\$SPAD/src/algebra modgcd.spad} \author{James Davenport, Patrizia Gianni} \maketitle \begin{abstract} \end{abstract} \eject \tableofcontents \eject \section{package INMODGCD InnerModularGcd} <>= )abbrev package INMODGCD InnerModularGcd ++ Author: J.H. Davenport and P. Gianni ++ Date Created: July 1990 ++ Date Last Updated: November 1991 ++ Description: ++ This file contains the functions for modular gcd algorithm ++ for univariate polynomials with coefficients in a ++ non-trivial euclidean domain (i.e. not a field). ++ The package parametrised by the coefficient domain, ++ the polynomial domain, a prime, ++ and a function for choosing the next prime Z ==> Integer NNI ==> NonNegativeInteger InnerModularGcd(R,BP,pMod,nextMod):C == T where R : EuclideanDomain BP : UnivariatePolynomialCategory(R) pMod : R nextMod : (R,NNI) -> R C == with modularGcdPrimitive : List BP -> BP ++ modularGcdPrimitive(f1,f2) computes the gcd of the two polynomials ++ f1 and f2 by modular methods. modularGcd : List BP -> BP ++ modularGcd(listf) computes the gcd of the list of polynomials ++ listf by modular methods. reduction : (BP,R) -> BP ++ reduction(f,p) reduces the coefficients of the polynomial f ++ modulo the prime p. T == add -- local functions -- height : BP -> NNI mbound : (BP,BP) -> NNI modGcdPrimitive : (BP,BP) -> BP test : (BP,BP,BP) -> Boolean merge : (R,R) -> Union(R,"failed") modInverse : (R,R) -> R exactquo : (BP,BP,R) -> Union(BP,"failed") constNotZero : BP -> Boolean constcase : (List NNI ,List BP ) -> BP lincase : (List NNI ,List BP ) -> BP if R has IntegerNumberSystem then reduction(u:BP,p:R):BP == p = 0 => u map(symmetricRemainder(#1,p),u) else reduction(u:BP,p:R):BP == p = 0 => u map(#1 rem p,u) FP:=EuclideanModularRing(R,BP,R,reduction,merge,exactquo) zeroChar : Boolean := R has CharacteristicZero -- exported functions -- -- modular Gcd for a list of primitive polynomials modularGcdPrimitive(listf : List BP) :BP == empty? listf => 0$BP g := first listf for f in rest listf | not zero? f while positive? degree g repeat g:=modGcdPrimitive(g,f) g -- gcd for univariate polynomials modularGcd(listf : List BP): BP == listf:=remove!(0$BP,listf) empty? listf => 0$BP # listf = 1 => first listf minpol:=1$BP -- extract a monomial gcd mdeg:= "min"/[minimumDegree f for f in listf] if positive? mdeg then minpol1:= monomial(1,mdeg) listf:= [(f exquo minpol1)::BP for f in listf] minpol:=minpol*minpol1 listdeg:=[degree f for f in listf ] -- make the polynomials primitive listCont := [content f for f in listf] contgcd:= gcd listCont -- make the polynomials primitive listf :=[(f exquo cf)::BP for f in listf for cf in listCont] minpol:=contgcd*minpol ans:BP := --one polynomial is constant member?(1,listf) => 1 --one polynomial is linear member?(1,listdeg) => lincase(listdeg,listf) modularGcdPrimitive listf minpol*ans -- local functions -- --one polynomial is linear, remark that they are primitive lincase(listdeg:List NNI ,listf:List BP ): BP == n:= position(1,listdeg) g:=listf.n for f in listf repeat if (f1:=f exquo g) case "failed" then return 1$BP g -- test if d is the gcd test(f:BP,g:BP,d:BP):Boolean == d0:=coefficient(d,0) coefficient(f,0) exquo d0 case "failed" => false coefficient(g,0) exquo d0 case "failed" => false f exquo d case "failed" => false g exquo d case "failed" => false true -- gcd and cofactors for PRIMITIVE univariate polynomials -- also assumes that constant terms are non zero modGcdPrimitive(f:BP,g:BP): BP == df:=degree f dg:=degree g dp:FP lcf:=leadingCoefficient f lcg:=leadingCoefficient g testdeg:NNI lcd:R:=gcd(lcf,lcg) prime:=pMod bound:=mbound(f,g) while zero? (lcd rem prime) repeat prime := nextMod(prime,bound) soFar:=gcd(reduce(f,prime),reduce(g,prime))::BP testdeg:=degree soFar zero? testdeg => return 1$BP ldp:FP:= ((lcdp:=leadingCoefficient(soFar::BP)) = 1) => reduce(lcd::BP,prime) reduce((modInverse(lcdp,prime)*lcd)::BP,prime) soFar:=reduce(ldp::BP *soFar,prime)::BP soFarModulus:=prime -- choose the prime while true repeat prime := nextMod(prime,bound) lcd rem prime =0 => "next prime" fp:=reduce(f,prime) gp:=reduce(g,prime) dp:=gcd(fp,gp) dgp :=euclideanSize dp if dgp =0 then return 1$BP if dgp=dg and not (f exquo g case "failed") then return g if dgp=df and not (g exquo f case "failed") then return f dgp > testdeg => "next prime" ldp:FP:= ((lcdp:=leadingCoefficient(dp::BP)) = 1) => reduce(lcd::BP,prime) reduce((modInverse(lcdp,prime)*lcd)::BP,prime) dp:=ldp *dp dgp=testdeg => correction:=reduce(dp::BP-soFar,prime)::BP zero? correction => ans:=reduce(lcd::BP*soFar,soFarModulus)::BP cont:=content ans ans:=(ans exquo cont)::BP test(f,g,ans) => return ans soFarModulus:=soFarModulus*prime correctionFactor:=modInverse(soFarModulus rem prime,prime) -- the initial rem is just for efficiency soFar:=soFar+soFarModulus*(correctionFactor*correction) soFarModulus:=soFarModulus*prime soFar:=reduce(soFar,soFarModulus)::BP dgp soFarModulus:=prime soFar:=dp::BP testdeg:=dgp if not zeroChar and euclideanSize(prime)>1 then result:=dp::BP test(f,g,result) => return result -- this is based on the assumption that the caller of this package, -- in non-zero characteristic, will use primes of the form -- x-alpha as long as possible, but, if these are exhausted, -- will switch to a prime of degree larger than the answer -- so the result can be used directly. merge(p:R,q:R):Union(R,"failed") == p = q => p p = 0 => q q = 0 => p "failed" modInverse(c:R,p:R):R == (extendedEuclidean(c,p,1)::Record(coef1:R,coef2:R)).coef1 exactquo(u:BP,v:BP,p:R):Union(BP,"failed") == invlcv:=modInverse(leadingCoefficient v,p) r:=monicDivide(u,reduction(invlcv*v,p)) not zero? reduction(r.remainder,p) => "failed" reduction(invlcv*r.quotient,p) -- compute the height of a polynomial -- height(f:BP):NNI == degf:=degree f "max"/[euclideanSize cc for cc in coefficients f] -- compute the bound mbound(f:BP,g:BP):NNI == hf:=height f hg:=height g 2*min(hf,hg) -- ForModularGcd(R,BP) : C == T -- where -- R : EuclideanDomain -- characteristic 0 -- BP : UnivariatePolynomialCategory(R) -- -- C == with -- nextMod : (R,NNI) -> R -- -- T == add -- nextMod(val:R,bound:NNI) : R == -- ival:Z:= val pretend Z -- (nextPrime(ival)$IntegerPrimesPackage(Z))::R -- -- ForTwoGcd(F) : C == T -- where -- F : Join(Finite,Field) -- SUP ==> SparseUnivariatePolynomial -- R ==> SUP F -- P ==> SUP R -- UPCF2 ==> UnivariatePolynomialCategoryFunctions2 -- -- C == with -- nextMod : (R,NNI) -> R -- -- T == add -- nextMod(val:R,bound:NNI) : R == -- ris:R:= nextItem(val) :: R -- euclideanSize ris < 2 => ris -- generateIrredPoly( -- (bound+1)::PositiveInteger)$IrredPolyOverFiniteField(F) -- -- -- ModularGcd(R,BP) == T -- where -- R : EuclideanDomain -- characteristic 0 -- BP : UnivariatePolynomialCategory(R) -- T ==> InnerModularGcd(R,BP,67108859::R,nextMod$ForModularGcd(R,BP)) -- -- TwoGcd(F) : C == T -- where -- F : Join(Finite,Field) -- SUP ==> SparseUnivariatePolynomial -- R ==> SUP F -- P ==> SUP R -- -- T ==> InnerModularGcd(R,P,nextMod(monomial(1,1)$R)$ForTwoGcd(F), -- nextMod$ForTwoGcd(F)) @ \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}