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