# #--> NGCD(A,B) GCD of A and B in Q(alpha)[X]. # # The inputs are multivariate polynomials over an algebraic number # field K = Q(alpha1,...,alphan) with zero or more field extensions. # (Algorithm MGCD can be used if n = 0, i.e. there are no field extensions). # Let G = GCD(A,B) be the monic GCD over A and B. # The output of NGCD(A,B) is the primitive associate G~ of G over Z. # E.g. if G = x + y/3 - sqrt(2)/2 which is monic (in x), # then G~ = 6 x + 2 y - 3 sqrt(2). # # This is an implementation of the dense modular GCD algorithm # of Brown [1] with rational reconstruction and field splitting. # We compute the GCD of A and B modulo primes p[1], p[2], ..., and apply # Chinese remaindering to the images and rational reconstruction # to recover G, the monic GCD of A and B, and use trial division # to estabilish the correctness of the GCD. # If the first field extension is of low degree, e.g. a quadratic # or cubic, the primes are chosen so that the minimal polynomial # splits into distinct linear factors. This allows us to eliminate # this field extension, which otherwise makes the modular algorithm # innefficient. # # NGCD calls AGCD which calls PGCD which calls the Euclidean algorithm. # NGCD chooses primes, computes images and uses rational reconstruction. # AGCD splits extensions of low degree and interpolates the images. # PGCD computes the GCD mod the primes via evaluation/interpolation # # Authors: M. Monagan (2000, 2001, 2002, 2003, 2004), C. Pastro (2000). # # Extension: MBM May/03 to handle reducible defining polynomials. # Extension: MBM Nov/04 to handle non sqrfree defining polynomials. # # The algorithm does NOT require that K be a field, that is, the defining # polyomials may be reducible. It is desirable, for example, to be able # to work directly with a number ring like K = Q((-1)^(1/3)). In this # example the defining polyomial, m(z) is z^3+1 = (z+1)*(z^2-z+1). # The algorithm will either output the gcd G~ or it will find a zero divisor, # d(z). If this happens an error is generated of the form # # ERROR( "zero divsior found", d(z) ) # # References # # [1] W. S. Brown. On Euclid's Algorithm and the # Computation of Polynomial Greatest Common Divisors, # Journal of the ACM, 18 (1971), pp. 476-504. # # [2] M. J. Encarnacion, # Computing GCDs of Polynomials over Algebraic Number Fields, # J. Symbolic Computation, 20 (1995), pp. 299-313. # # [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, 2000. # # [4] M. van Hoeij, M. Monagan. A Modular GCD algorithm over # Number Fields presented with Multiple Extensions, # Proceedings of ISSAC '2002, ACM Press, pp. 109-116, 2002. # NEXTPRIME := proc(p) option inline; modp1(Prime(p)) end: NGCD := proc(A::POLYNOMIAL,B::POLYNOMIAL) local R,F,K,M,X,E,N,a,b,gamma,d,m,gm,p,ap,bp,gp,h, G,t,P,EE,GG,DD,PP,NN,RR,MP, zerodeg,deggp,i,k,g,terminating,minpolys,splitpolys,x,y,format,defect, zerodivisors, zds, n, mzd, z, pt, gt, dt, rt, tt, st; pt := 0; gt := 0; dt := 0; rt := 0; tt := time(); if getring(A)<>getring(B) then ERROR("inconsistent rings") fi; if iszerorpoly(A) then RETURN( ipprpoly(pprpoly(B)) ) fi; if iszerorpoly(B) then RETURN( ipprpoly(pprpoly(A)) ) fi; R := getring(A); M,X,E := op(R); ASSERT(M=0); # if M<>0 then PGCD is the procedure to use M := nops(E); # the number of field extensions N := nops(X)-M; # the number of polynomial variables K := getcofring(R); # K = Q(alpha1,...,alphan) # Do not invert the lc(a) or lc(b) - this may cause a blowup a := ipprpoly(A); b := ipprpoly(B); zerodeg := [seq(0,i=1..N)]; MP := [seq(getalgext(R,-i),i=1..M)]; minpolys := map(convert,MP,POLYNOMIAL); format := cat("NGCD: GCD in Q%a%a mod <", "%a, "$(M-1), "%a>\n"); printf(format,X[N+1..N+M],X[1..N],seq(minpolys[i],i=1..M)); # lc-bad test: p must not divide gamma gamma := ilcm( seq(ilcrpoly(p), p=MP), seq(igcd(icoeffsrpoly(lcrpoly(p))), p=[a,b]) ); #print(GAMMA=normrpoly(lcrpoly(a)),DELTA=defectrpoly(K)); splitpolys,minpolys := whichsplit(minpolys,MP,X); if splitpolys<>[] then format := cat("NGCD: I will split ", "%a, "$(nops(splitpolys)-1), "%a\n"); printf(format,op(splitpolys)); fi; if minpolys<>[] then format := cat("NGCD: I won't split ", "%a, "\$(nops(minpolys)-1), "%a\n"); printf(format,op(minpolys)); fi; m := 1; # the product of the moduli p := 1; for k do p := NEXTPRIME(p); while gamma mod p = 0 or not primesplits(splitpolys,X,p) do printf("%d ... ",p); p := NEXTPRIME(p); od; printf(" NGCD: Prime %d = %d\n", k, p); st := time(); ap, bp := phirpoly(a,p), phirpoly(b,p); pt := pt+time()-st; st := time(); gp := traperror( AGCD(ap,bp) ); gt := gt+time()-st; if gp=lasterror then if op(1,[gp]) = "inverse does not exist" or op(1,[gp]) = inverse does not exist then # test for zero divisor # assemble the zero divisor as a POLYNOMIAL from the error data printf(" NGCD: fail prime %d\n",p); gp := gp[2]; if type(gp,list) then gp := POLYNOMIAL( [gp[1],X[-gp[2]..-1],gp[3]], gp[4] ); fi; if not assigned(zerodivisors) then zerodivisors := [] fi; zerodivisors := [op(zerodivisors),gp]; zds := map( liftrpoly, zerodivisors ); n := nops(getvars(zds[-1])); zds := select( proc(z) nops(getvars(z))=n end, zds ); n := degrpoly(zds[-1]); zds := select( proc(z) degrpoly(z)=n end, zds ); n := nops(getvars(zds[-1])); mzd := MP[n]; # zds[-1] should divide mzd mod p if not divrpoly(phirpoly(mzd,p),zds[-1]) then ERROR("should not happen") fi; if isFibonacci(nops(zds)) then zds := map( retextsrpoly, keeplastFn( zds ) ); zds := irrrpoly( ichremrpoly( zds ) ); if zds <> FAIL then zds := ipprpoly( subsop(1=op(1,mzd),zds) ); if divrpoly(mzd,zds) then #ERROR("zero divisor found",rpoly(mzd),rpoly(zds)); ERROR("zero divisor found",phirpoly(zds,mzd)); fi; fi; fi; else ERROR(lasterror) fi; else gp := retextsrpoly(gp); # Drop extensions for Chinese remaindering deggp := degrpoly(gp,X[1..N]); if deggp=zerodeg then # GCD is one RETURN(rpoly(1,R)) elif m=1 then # first time, don't chrem, just update m,gm,d,h := p,gp,deggp,irrrpoly(gp); 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,h := p,gp,deggp,irrrpoly(gp); else # termination tests if phirpoly(liftrpoly(modsrpoly(gm)),p)=gp then h := liftrpoly(modsrpoly(gm)); g := subsop(1=R,h); printf("NGCD: Trial divisions over Z starting after %d primes\n",k); st := time(); if divrpoly(a,g) and divrpoly(b,g) then dt := dt+time()-st; tt := time()-tt; print(GCD=100*gt/tt,REC=100*rt/tt,DIV=100*dt/tt,PHI=100*pt/tt); RETURN(g) fi; printf("NGCD: Trial divisions failed\n"); fi; if h<>FAIL and phirpoly(h,p)=gp then # begin termination test g := subsop(1=R,ipprpoly(h)); printf("NGCD: Trial divisions over Q starting after %d primes\n",k); st := time(); if divrpoly(a,g) and divrpoly(b,g) then dt := dt+time()-st; tt := time()-tt; print(GCD=100*gt/tt,REC=100*rt/tt,DIV=100*dt/tt,PHI=100*pt/tt); RETURN(g) fi; printf("NGCD: Trial divisions failed\n"); fi; st := time(); gm := ichremrpoly([gm,gp]); # Chinese remaindering h := irrrpoly(gm); # Rational reconstruction m := m*p; # Update the modulus rt := rt+time()-st; fi; fi; od; end: #macro( Phi=NGCDnf/Phi): #Phi := proc() local m,f,T; # T := table(sparse); # for m to 30 do # f := numtheory[cyclotomic](m,_X); # if degree(f,_X)<=6 then T[f] := m fi; # od; # eval(T) # end(): ## Note, the magic constant 6 is explicit in AGCD too primesplits := proc(minpolys,X,p) local a,b,c,d,m,x,g,j,w,R,r,k; if minpolys=[] then RETURN(true) fi; m := minpolys[1]; x := X[-1]; d := degree(m,x); if d=2 then (a,b,c) := coeff(m,x,2),coeff(m,x,1),coeff(m,x,0); # m = a*x^2+b*x+c splits into distinct linear factors iff # (i) p does not divide m otherwise we two equal roots # (ii) the disciminant of m is a perfect square mod p j := numtheory[jacobi](b^2-4*a*c,p); if j=0 or j=-1 then RETURN(false) fi; elif d=3 and coeff(m,x,1)=0 and coeff(m,x,2)=0 and coeff(m,x,3)=1 then # m = x^3+a splits ==> sqrt(-3) \in Zp ==> J(-3/p) = J(p/3) = 1 ==> p == 1 mod 3 j := mods(p,3); if j=0 or j=-1 then RETURN(false) fi; k := iquo(p-1,3); a := coeff(m,x,0); if mods( a &^ k, p ) <> 1 then RETURN(false) fi; # elif Phi[subs(x=_X,m)] > 0 then # if irem(p, Phi[subs(x=_X,m)]) <> 1 then RETURN(false) fi; else if d=3 or d=6 then j := mods(p,3); if j=0 or j=-1 then RETURN(false) fi; fi; # m splits in Zp[x] into distinct linear factors iff # (i) GCD(m,m')=1 i.e., no repeated roots and # (ii) GCD(x^p-x,m)=m upto a scalar) g := Gcd(m,diff(m,x)) mod p; if g<>1 then RETURN(false) fi; w := Powmod(x,p,m,x) mod p; g := Gcd(m,w-x) mod p; if degree(g,x) 0 then # [m1,m2],minpolys[3..-1]; # Prob = 1/2/n where n=deg[m2] elif degree(m2,y)=4 and separable(r2) then if not has(m2,x) and coeff(m2,y,1)=0 and coeff(m2,y,3)=0 and coeff(m2,y,0)=coeff(m2,y,4) then [m1,m2],minpolys[3..-1] # Prob = 1/8 else [m1],minpolys[2..-1]; # Prob = 1/2 fi; else [m1],minpolys[2..-1]; # Prob = 1/2 fi; fi; elif degree(m1,x)=3 and separable(r1) then if nops(minpolys)=1 then [m1],[] # Prob = 1/6 else m2,r2,y := minpolys[2],minrpolys[2],X[-2]; if degree(m2,y)=2 and not has(m2,x) and separable(r2) then [m1,m2],minpolys[3..-1] # Prob = 1/12 else [m1],minpolys[2..-1]; # Prob = 1/6 fi; fi; # elif Phi[ subs(x=_X,m1) ] > 0 then # # n = deg(m1,x) # if nops(minpolys)=1 then [m1],[] # Prob = 1/n # else # m2,y := minpolys[2],X[-2]; # if degree(m2,y)=2 and not has(m2,x) # then [m1,m2],minpolys[3..-1] # Prob = 1/2/n # else [m1],minpolys[2..-1]; # Prob = 1/n # fi; # fi; elif degree(m1,x)=4 and coeff(m1,x,1)=0 and coeff(m1,x,3)=0 and separable(r1) then if nops(minpolys)=1 then [m1],[] # Prob = 1/4 or 1/8 else m2,r2,y := minpolys[2],minrpolys[2],X[-2]; if degree(m2,y)=2 and not has(m2,x) and coeff(m1,x,0)=coeff(m1,x,4) and separable(r2) then [m1,m2],minpolys[3..-1] # Prob = 1/8 else [m1],minpolys[2..-1]; # Prob = 1/4 or 1/8 fi; fi; else [],minpolys fi; end: normrpoly := proc(a::POLYNOMIAL,K::POLYRING) local F,R,b,m,N,r,X,x; if nargs=1 then F := [getchar(a),[],[]]; RETURN( normrpoly(a,F) ) fi; R := getcofring(a); if not issubring(K,R) then ERROR("invalid input") fi; if R[3]=K[3] then RETURN(a) fi; b := liftrpoly(a); m := getalgext(a); m := scarpoly(invrpoly(lcrpoly(m)),m); # makes the minimal polynomial monic X := getpolvars(b); N := nops(X); x := X[-1]; if N > 1 then R := getring(b); R := [R[1],subsop(1=X[N],N=X[1],R[2]),R[3]]; b := convert(b,POLYRING,R); m := convert(m,POLYRING,R); r := mresrpoly(m,b); else r := resrpoly(m,b); fi; normrpoly(r,K); end: pnormrpoly := proc(a::POLYNOMIAL,K::POLYRING) # This norm is a root of the previous one -- often this will be preferable. local F,R,b,m,r; if nargs=1 then F := [getchar(a),[],[]]; RETURN( pnormrpoly(a,F) ) fi; R := getcofring(a); if not issubring(K,R) then ERROR("invalid input") fi; if R[3]=K[3] then RETURN(a) fi; b := liftrpoly(a); m := getalgext(a); if degrpoly(b)=0 then RETURN( pnormrpoly( retvarpoly(b), K ) ) fi; r := resrpoly(m,b); pnormrpoly( r, K ) end: charrpoly := proc(a::POLYNOMIAL,x::name,K::POLYRING) local F,R,b; if nargs=2 then F := [getchar(a),[],[]]; RETURN(charrpoly(a,x,F)) fi; R := getcofring(a); if R <> getring(a) then ERROR("first argument must not be a polynomial") fi; if not issubring(K,R) then ERROR("second argument must be a subring of the first") fi; # Form x-a b := list2rpoly([negrpoly(a),rpoly(1,R)],x); normrpoly(b,K) end: minrpoly := proc(a::POLYNOMIAL,x::name,K::POLYRING) local F,c; if nargs=2 then F := [getchar(a),[],[]]; RETURN(minrpoly(a,x,F)) fi; c := charrpoly(a,x,K); # Take the square-free part of c quorpoly(c,gcdrpoly(c,diffrpoly(c))); end: defectrpoly := proc(F::POLYRING) local M,X,E,l,m,r,K; M,X,E := op(F); if M<>0 then ERROR("defect if only defined for characteristic 0") fi; if nops(X)<>nops(E) then ERROR("input cannot be a polynomial") fi; if nops(E)=0 then RETURN( rpoly(1,F) ) fi; m := getalgext( F, 1 ); K := getcofring( F, 2 ); r := resrpoly( diffrpoly(m), m ); l := lcrpoly(m); if type(l,rational) then l := rpoly(l,K) fi; r := mulrpoly(invrpoly(l),r); mulrpoly( pnormrpoly(r), defectrpoly(K) ); end: isdefectzero := proc(F::POLYRING,p::integer) local K,i,m,l,d; if nargs=3 then d := args[3]; if not getring(d)=F then ERROR("third argument in wrong ring") fi; else d := rpoly(1,F) fi; K := F; for i to nops(F[3]) do m := getalgext(F,-i); l := lcrpoly(m); if type(l,rational) then l := POLYNOMIAL([0,[],[]],l) fi; if iszerorpoly(phirpoly(l,p)) then RETURN(true) fi; od; K := phirpoly(F,p); d := phirpoly(d,p); if iszerorpoly(d) then RETURN(true) fi; while K[3] <> [] do d := liftrpoly(d); m := getalgext(K,1); K := getcofring(K,2); # N(d)*N(Discrim(K)) | N(res(d,m)*res(m',m)) = N(res(d*m',m)) d := traperror( resrpoly(mulrpoly(d,diffrpoly(m)),m) ); if d = inverse does not exist then printf(" cannot compute defect "); RETURN(true) elif d = lasterror then ERROR(lasterror) fi; if iszerorpoly(d) then printf(" defect might be zero "); RETURN(true) fi; od; false; end: lcbad := proc(f,p) local g; if type(f,{list,set}) then for g in f do if lcbad(g,p) then RETURN(true) fi; od; false; else evalb( {icoeffsrpoly(lcrpoly(f))} mod p = {0} ) fi; end: isFibonacci := proc(n::nonnegint) local Fn,Fp; if n=0 then RETURN(true) fi; Fn,Fp := 1,1; while n>Fn do Fn,Fp := Fn+Fp,Fn od; evalb( n=Fn ); end: keeplastFn := proc(L) local n,Fn,Fp,m; # Let n = nops(L) where n > 0. # Assuming nops(L) = Fn for some fibonacci number Fn, # discard the first Fn-2 entries of L. n := nops(L); Fn,Fp := 1,1; while n>Fn do Fn,Fp := Fn+Fp,Fn; od; m := Fn-Fp; L[m+1..n]; end: