# PGCD depends on recden readlib(ilog): # delete for Maple 6 # #--> PGCD(A,B) compute the GCD of A and B mod p. # # Here A and B are multivariate polynomials over Fp or # over Fp with one or more algebraic extensions. # This is an implementation of the dense modular algorithm of # Brown with trial division. A description of the basic algorithm # may be found in ``Algorithms for Computer Algebra'' by # Geddes, Labahn and Czapor, Kluwer Academic Publishers. # When the coefficient ring is a small finite field, # i.e. the input is in GF(q)[x1,...,xn], q < deg[x1](G), G = GCD(A,B), # then we compute the GCD mod m[1](x1), m[2](x1), ..., where # m[i] is irreducible over GF(q), and apply the Chinese remainder theorem. # This is more efficient than introducing a field extension using # a new (dummy) variable. # # Author: Michael Monagan, 2000, 2001. # # References # [1] W. S. Brown. On Euclid's Algorithm and the # Computation of Polynomial Greatest Common Divisors. # Journal of the ACM, 18, pp. 476-504, 1971 # [2] E. Kaltofen, M.B. Monagan. # On the Genericity of the Modular Polynomial GCD Algorithm. # Proceedings of ISSAC '99, ACM Press, pp. 59-66, 1999. # [3] M. B. Monagan, A. D. Wittkopf, # On the Design and Implementation of Brown's Algorithm # over the Integers and Number Fields, # Proceedings of ISSAC '2000, ACM Press, pp. 225-233. # PGCD := proc(A::POLYNOMIAL, B::POLYNOMIAL) local R,M,X,E,N,a,b,F,zerodeg,zero,one,conta,contb,gcdCont,gamma,gammap,m, alpha,alphapol,C,ap,bp,gp,deggp,xn,mp,gm,d,h,g,i,extend,l,deg1,deg2, bound,k,q,terminating,minpolys,format; R := getring(A); M,X,E := op(R); ASSERT(M<>0); # else MGCD is the procedure to use N := nops(X)-nops(E); if N<2 then RETURN(gcdrpoly(A,B)) fi; # single variable gcd # if A or B is zero then return the other checkconsistency(A,B); if iszerorpoly(A) then RETURN(pprpoly(B)) fi; if iszerorpoly(B) then RETURN(pprpoly(A)) fi; if nops(E)>0 then minpolys := [seq(mods(convert(getalgext(R,-i),POLYNOMIAL),M),i=1..nops(E))]; format := cat(" PGCD: GCD in GF(%a)%a%a where ", "%a, "\$(nops(E)-1), "%a\n"); printf(format,M,X[N+1..-1],X[1..N],op(minpolys)); else printf(" PGCD: Gcd in GF(%a)%a\n", M,X); fi; a,b := A,B; extend := false; if E=[] then F := M else F := [M,X[-nops(E)..-1],E] fi; zerodeg := [0\$(N-1)]; # Calculate the content of a & b considered multivariate polynomials conta, contb := contrpoly(a,X[1..N-1],'a'), contrpoly(b,X[1..N-1],'b'); gcdCont := gcdrpoly(conta,contb); # univariate gcd in X[N] over F # Calculate the leading coefficient correction coefficient gamma. gamma := gcdrpoly(lcrpoly(a,X[1..N-1]),lcrpoly(b,X[1..N-1])); m := 1; # the product of polynomial moduli if E=[] then zero := 0; alpha := 0; alphapol := 0; C := M; else alphapol := POLYNOMIAL(subsop([3,1]=NULL,F),0); alpha := subsop(1=F,alphapol); zero := POLYNOMIAL(F,0); C := monrpoly(subsop([3,1]=NULL,F),[nops(E[1])-1,0\$nops(E)-1]); fi; terminating := false; do # Choose alpha such that gamma(alpha) mod p <> 0 while alphapol <> C and evalrpoly(gamma,X[N]=alpha)=zero do if E=[] then alpha := alpha+1; alphapol := alpha; else alphapol := nrpoly(alphapol); alpha := subsop(1=F,alphapol); fi; od; if alphapol=C then # EXTEND THE FIELD if extend=false then # first time to extend the field extend := true; # hmm how big should we extend the field # we need to be able to interpolate the current variable X[N] bound := min(degrpoly(a,X[N]), degrpoly(b,X[N])); if m<>1 then bound := bound - degrpoly(m) fi; # we would like to interpolate X[2],...,X[N-1] without further extensions bound := 1 + max(bound,seq(min(degrpoly(a,X[i]),degrpoly(b,X[i])),i=2..N-1)); q := M; for i to nops(E) do q := q^(nops(E[-i])-1) od; k := max(2,1+ilog[q](bound+9)); # k = ceil(log[q](bound+10)); #if q=M then k := max(k,8) fi; mp := monrpoly([M,X[N..-1],E],[k,0\$nops(E)]); # mp = x^k fi; # make sure phirpoly(gamma,mp)<>0 mp := nmirpoly(mp); while op(2,phirpoly(gamma,mp))=0 do mp := nmirpoly(mp) od; printf(sprintf(" ... mod <%a>, a field extension of GF(%d)\n", mods(convert(mp,POLYNOMIAL),M),q)); ap := phirpoly(a, mp); bp := phirpoly(b, mp); gp := PGCD(ap,bp); # normalize gp to avoid the leading coefficient problem gammap := phirpoly(gamma,mp); gp := scarpoly(gammap,gp); else # evaluate the polynomials at alpha if type(alpha,integer) then printf(" ... mod <%a-%a>\n",X[N],alpha); else printf(" ... mod <%a-%a>\n",X[N],mods(convert(alpha,POLYNOMIAL),M)); fi; ap := evalrpoly(a, X[N]=alpha); bp := evalrpoly(b, X[N]=alpha); gp := PGCD(ap,bp); # normalize gp to avoid the leading coefficient problem gammap := evalrpoly(gamma, X[N]=alpha); gp := scarpoly(gammap, gp); if E=[] then mp := POLYNOMIAL([M,[X[N]],E],[-alpha,1]); else xn := monrpoly([M,[X[N],op(N+1..nops(X),X)],E],[1,0\$nops(E)]); mp := subrpoly(xn,POLYNOMIAL([M,[X[N],op(N+1..nops(X),X)],E],[op(2,alpha)])); fi; # we lost a variable when we evaluated so add X[N] back into the polynomial # do this by using extrpoly gp := extrpoly(gp,mp); fi; deggp := degrpoly(gp,X[1..N-1]); if deggp=zerodeg then # gp is a constant one := monrpoly(R,[0\$N,0\$nops(E)]); RETURN(scarpoly(gcdCont,one)) fi; if m=1 then # first time around.. no need to chrem... just update m,gm := mp,gp; d := deggp; elif compdeg(deggp,d)=greater then # we have an unlucky image... ignore it elif compdeg(deggp,d)=less then # discard previous images and update d m,gm := mp,gp; d := deggp; #elif terminating and phirpoly(liftrpoly(h),mp) = gp then # g := pprpoly(liftrpoly(h),X[1..N-1]); # make g primitive # RETURN( scarpoly(gcdCont,g) ); else # Chinese Remaindering h := chremrpoly([gm,gp]); if op(2,h)=op(2,gm) then # we can begin our "termination test" #terminating := true; g := pprpoly(liftrpoly(h),X[1..N-1]); # make g primitive if divrpoly(a,g) and divrpoly(b,g) then RETURN( scarpoly(gcdCont,g) ) else print("Trial Divisions failed: continuing") fi; fi; # continue interpolating # update m and gm m := getalgext(h); gm := h; fi; if extend=false then # go to the next alpha value if E=[] then alpha := alpha+1; alphapol := alpha; else alphapol := nrpoly(alphapol); alpha := subsop(1=F,alphapol); fi; fi; od; end: