# MGCD depends on recden and PGCD # This version does rational reconstruction + special case for monic GCD with integer coefficients. # MBM Jul/03 NEXTPRIME := proc(p) option inline; modp1(Prime(p)) end: MGCD := proc(A::POLYNOMIAL,B::POLYNOMIAL) local R,I,M,X,E,N,a,b,conta,contb,gcdCont,gamma,d,m,gm,p,ap,bp,gp,h,zerodeg,deggp,i,k,g; checkconsistency(A,B); R,I := op(A); M,X,E := op(R); N := nops(X)-nops(E); ASSERT(M=0); # if M<>0 then PGCD is the procedure to use # if A or B is zero then return the other if op(2,A)=0 then RETURN(B) fi; if op(2,B)=0 then RETURN(A) fi; zerodeg := [seq(0,i=1..N)]; a := ipprpoly(A); b := ipprpoly(B); gamma := igcd(ilcrpoly(a),ilcrpoly(b)); printf("MGCD: Gcd in Q%a, log[10](gamma) = %d.\n", X, length(gamma)); m := 1; # the product of the moduli p := 1; for k do # Choose p such p does NOT divide gamma p := NEXTPRIME(p); while irem(gamma,p)=0 do p := NEXTPRIME(p) od; printf("MGCD: Prime %d = %d ... ", k, p); gp := PGCD(phirpoly(a,p),phirpoly(b,p)); printf(" Gcd has (total) degree %d\n",tdegrpoly(gp)); deggp := degrpoly(gp,X[1..N]); if deggp=zerodeg then # gp is a constant RETURN(rpoly(1,R)) elif m=1 then # first time, don't chrem, just update m,gm,d := p,gp,deggp; elif compdeg(deggp,d)=greater then # do nothing, unlucky image elif compdeg(deggp,d)=less then # discard previous images and update d m,gm,d := p,gp,deggp; else h := ichremrpoly([gm,gp]); # Chinese Remaindering if op(2,modsrpoly(h))=op(2,modsrpoly(gm)) then # we can begin our "termination test" g := liftrpoly(modsrpoly(h)); printf("MGCD: Trial divisions over Z starting after %d primes\n",k); if divrpoly(a,g) and divrpoly(b,g) then RETURN(g) fi; printf("MGCD: trial divisions failed\n"); fi; # else keep going gm := h; m := m*p; if gamma = 1 then next fi; h := irrrpoly(h,gamma); # Rational reconstruction if h <> FAIL then g := ipprpoly(liftrpoly(h)); # clear fractions printf("MGCD: Trial divisions over Q starting after %d primes\n",k); if divrpoly(a,g) and divrpoly(b,g) then RETURN(g) fi; printf("MGCD: trial divisions failed\n"); fi; fi; od; end: #--------------------------------------------------------------------------- # compdeg(p::list,q::list)::{less,greater,equal} # # Compares two vector degrees p & q and returns: # "less" if p < q, # "greater" if p > q and # "equal" if p = q compdeg := proc(p,q) if nops(p)<>nops(q) then ERROR(`Comparing different size vector degrees`,p,q) fi; if p=[] and q=[] then RETURN(equal) elif p[1] < q[1] then RETURN(less) elif p[1] > q[1] then RETURN(greater) else compdeg(p[2..-1],q[2..-1]) fi; end: