# AGCD computes a GCD in Zp[a]/m(a)[b]/m(b,a)...[x,y,z,...]
# If m(a), the first field extension, splits into distinct linear factors
# then AGCD will try to compute the GCD by computing the GCD modulo
# the linear factors an use the CRT to reconstruct alpha.
# If so, AGCD calls PGCD k=deg[z](m(z)) times otherwise it calls PGCD once.
# The reason we do this is because computing mod m(a) when k is small
# is very innefficient because almost all of the time is spent doing
# storage management instead of a few arithmetic operations mod p.
# Author: Michael Monagan, December 2000, 2003.

AGCD := proc(A::POLYNOMIAL,B::POLYNOMIAL)
local R,M,X,E,N,Re,Ie,Me,Xe,Ee,ae,be,ext,EXT,x,degext,rt,ep,i,g,U,G,l,deggi,d,w;

  checkconsistency(A,B);
  R := getring(A);
  M,X,E := op(R);
  N := nops(X)-nops(E);
  if M=0 then ERROR("inputs cannot be over characteristic 0") fi;
  if nops(E)=0 then RETURN( PGCD(A,B) ) fi; # no extensions

  # if A or B is zero then return the primitive part of the other
  if iszerorpoly(A) then RETURN(pprpoly(B,X[1..N-1])) fi; 
  if iszerorpoly(B) then RETURN(pprpoly(A,X[1..N-1])) fi;  

  # this should get the last extension if more than one extension
  x,ext := X[-1],getalgext(A,-1); # for now this is okay
  # need some information about this extension so
  Re,Ie := op(ext); Me,Xe,Ee := op(Re);
  degext := degrpoly(ext,x); # degree of the extension in it's main var.
  # if degext>4 then no use to split just use PGCD
  if degext>4 then RETURN( PGCD(A,B) ); fi; # else we will try to split
  #xpol := monrpoly(Re,[1]);
  #if degrpoly(gcdrpoly(ext,diffrpoly(ext)))>0 # Gcd(m,m')<>1
  #   or degrpoly(gcdrpoly(ext,subrpoly(powmodrpoly(xpol,M,ext),xpol))))<degext
  EXT := mods( convert(ext,POLYNOMIAL,x), M );
  if Gcd(EXT,diff(EXT,x)) mod M <> 1 then RETURN(PGCD(A,B)) fi;
  w := Powmod(x,M,EXT,x) mod M;
  if degree( Gcd(EXT,w-x) mod M , x ) < degext then RETURN(PGCD(A,B)) fi;
  # else we could split into linear factors so lets get a list of the evaluation points 

  rt := Roots(EXT) mod M;
  if nops(rt) < degext then ERROR("should not happen") fi;
  ep := map(c -> op(1,c),rt);
  # make sure that the leading coeffs don't vanish
  for i to degext do
    l := evalrpoly(lcrpoly(A),Xe[-1]=op(i,ep));
    if type(l,POLYNOMIAL) then l := op(2,l) fi; # else just an integer
    if l=0 then 
       printf("   AGCD: lc bad evaluation: calling PGCD\n");
       RETURN( PGCD(A,B) ); fi; # Bad evaluation point 
  od;

printf("  AGCD: %a = %a mod %a\n",EXT,mul(modp(x-i,M),i=ep),M);

  # Evaluate A and B at the evaluation points
  for i to degext do
    printf("   AGCD: | %a = %a\n", Xe[-1], ep[i]);
    ae := evalrpoly(A,Xe[-1]=op(i,ep));
    be := evalrpoly(B,Xe[-1]=op(i,ep));
    g[i] := traperror( AGCD(ae,be) );
    # This recursive call to AGCD may hit a zero divisor.
    # In order for NGCD to be able to reconstruct a possible
    # factor of a ``minimal polynomial'' in R it will need the
    # zero divisor over the original ring R.  We do not attempt
    # to interpolate the zero divisor over R but instead we just
    # call PGCD which will generate the right error information.
    if g[i] = lasterror then 
       printf("   AGCD: hit zero divisor: calling PGCD\n");
       RETURN( PGCD(A,B) );
    fi;
    # degree check; the degree's of the gi's must be the same
    deggi := degrpoly(g[i],X[1..N]);
# IF g[i] has degree 0 then we can stop ONLY if coefficient ring is a field
# i.e. has no zero divisors in characteristic 0.  Otherwise we must test if
# g[i] = 1 for all evaluation points before we can stop.
    if i=1 then d := deggi;
    elif compdeg(d,deggi)=less or compdeg(d,deggi)=greater then
       # We have an unlucky evaluation, for example
       # let m = z^3-z+p = z*(z-1)*(z+1) mod p and
       # let A = x^2+x+z and B = x^2+x-z, then for
       # z =  0 we have  GCD(A,B) mod p = x^2+x and for
       # z =  1 we have GCD(A,B) mod p = 1 and for
       # z = -1 we have GCD(A,B) mod p = 1 .
       # There is no way here to recover so compute the GCD
       # over the finite ring Zp[alpha]/<m(alpha)> where p=M
       printf("   AGCD: unlucky evaluation: calling PGCD\n");
       RETURN(PGCD(A,B))
    fi;
    # drop the extensions in order to interpolate
    g[i] := subsop([1,3]=[],g[i]);
  od;

printf("  AGCD: images done for %a - interpolating ...\n",EXT);

  U := [seq(g[i],i=1..degext)];
  G := interprpoly(U,ep,Xe[-1]);
  # put the ring information back in
  G := subsop(1=R,G);
  # if divrpoly(A,G) and divrpoly(B,G) then RETURN(G) fi;
  # print(`trial divisions failed in AGCD`);
  # These trial divisions are not necessary!!
  # Given f = z-a and g = z-b in Zp[z] where a <> b then
  # G|A and G|B mod f and G|A and G|B mod g ==> G|A and G|B mod fg
  RETURN(G);

end:

# For testing: number field case
#
# AGCD depends on recdeneval and PGCD
#read recden;
#read PGCD;
## two extensions
#p := 9973:
#M1 := RootOf(x^2+3*x+2):
#M2 := RootOf(y^2+7*y+12):
##M2 := RootOf(y^2+M1+5):
#d := 4: X := [M1,M2,x,y]:
#g := randpoly(X,degree=d,dense):
#a := expand(randpoly(X,degree=d,dense)*g): a := a mod p: 
#b := expand(randpoly(X,degree=d,dense)*g): b := b mod p:
#X := [x,y,w2,w1]:
#A := convert(a,POLYNOMIAL,X,[w1=M1,w2=M2],p):
#B := convert(b,POLYNOMIAL,X,[w1=M1,w2=M2],p):
#st := time(): g := Gcd(a,b) mod p: mt := time()-st;
##g := subs( {M1=w1}, Primpart(g,X[1]) mod p ):
#read(profile):
#profile(AGCD,PGCD,ichremrpoly,divrpoly,irrrpoly,evalrpoly):
#st := time(): G := AGCD(A,B); at := time()-st; 
#GG :=  expand(convert(G,POLYNOMIAL)):
#showprofile();
##evalb(normal(g-GG)=0 or normal(g+GG)=0);
#print(MT=mt,AT=at);
#
## one extension
#M1 := RootOf(x^2+5):
##M2 := RootOf(y^2+5):
#p := 7: d := 4: X := [M1,x,y]:
#g := randpoly(X,degree=d,dense):
#a := expand(randpoly(X,degree=d,dense)*g): a := a mod p: 
#b := expand(randpoly(X,degree=d,dense)*g): b := b mod p:
#X := [x,y,w1]:
#A := convert(a,POLYNOMIAL,X,[w1=M1],p):
#B := convert(b,POLYNOMIAL,X,[w1=M1],p):
#st := time(): g := Gcd(a,b) mod p: time()-st;
#g := subs( {M1=w1}, Primpart(g,X[1]) mod p ):
##resetprofiler();
##profile(AGCD,PGCD,ichremrpoly,divrpoly,irrrpoly,evalrpoly):
#st := time(): G2 := AGCD(A,B): time()-st; GG :=  expand(convert(G2,POLYNOMIAL)):
##showprofile();
#evalb(normal(g-GG)=0 or normal(g+GG)=0);
#
## one extension
#p := 9973:
#M1 := RootOf( mul(x+i,i=1..8) ) mod p;
#d := 4: X := [M1,x,y]:
#g := randpoly(X,degree=d,dense):
#a := expand(randpoly(X,degree=d,dense)*g): a := a mod p: 
#b := expand(randpoly(X,degree=d,dense)*g): b := b mod p:
#X := [x,y,w1]:
#A := convert(a,POLYNOMIAL,X,[w1=M1 mod p],p):
#B := convert(b,POLYNOMIAL,X,[w1=M1 mod p],p):
#st := time(): G2 := AGCD(A,B): time()-st; GG :=  expand(convert(G2,POLYNOMIAL)):
#st := time(): g := Gcd(a,b) mod p: time()-st;
#g := subs( {M1=w1}, Primpart(g,[X[1],X[2]]) mod p ):
#evalb(normal(g-GG) mod p =0 or normal(g+GG) mod p =0);