# 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: