# 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)))) 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]/ 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);