# To do: # -Simple optimizations in the modular portion of cgcd, such as terminating # when gcd = 1 modulo a prime. # RECDEN := Algebraic:-RecursiveDensePolynomials; randtriset := proc(X, d) # Create a random triangular set: # T[1] will be in X[-1] with degree d[1]. # T[2] will be in X[-1],X[-2] with degree d[2]. # etc. # # All polies will be monic in their main variable. local n,T,R,i; n := nops(d); T := [seq(1..n)]; C := rand(1..99); T[1] := X[n]^d[1] + randpoly(X[n], degree=d[1]-1) + C(); for i from 2 to n do T[i] := X[n-i+1]^d[i] + randpoly(X[n-i+1..n], degree=d[i]-1) + C(); od; return construct_triangular_set(T,X); end proc: randpolyn := proc(X,D) local i,C,f,n; n := nops(X); if n=1 then C := rand(1..99); f := add(C()*X[1]^i,i=0..D[1]-1); else f := add( randpolyn(X[2..-1],D[2..-1])*X[1]^i,i=0..D[1]-1); fi; end: randtriset := proc(X,D) T := []; n := nops(X); Y := [seq(X[n-i],i=0..n-1)]; for i from 1 to n do x := Y[i]; d := D[i]; T := [op(T),x^d+(randpolyn(Y[1..i],D[1..i]))]; od; return construct_triangular_set(T,X); end: reduce_by_tset := proc(f,rT) # Input :: f := a polynomial in the recden format. # T := a triangular set in the recden format. # # This method reduces f modulo T. It will work over Q or Z/p. local F,t; F := f; for t in rT do F := phirpoly(F,t); od; return F; end proc: mvarrpoly := proc(a::POLYNOMIAL) local i,n; # if tdegrpoly(a) <= 0 then # return 0; # fi; n := nops(getvars(a)); i := 1; while degrpoly(a,i) = 0 do i := i + 1; od; return i; end proc: construct_triangular_set := proc(T, Z) # Input :: T := A triangular set with variables in Z. # Z := A list of variables. # # It is assumed that T is of the form t1(Z[-1]), t2(Z[-1],Z[-2]), etc, and that # the polynomials in T are in Maple's format. # # This method also turns T into a reduced triangular set. local rT,i,n; n := nops(Z); rT := [rpoly(T[1],Z[-1]), seq(1..n-1)]; for i from 2 to nops(Z) do rT[i] := rpoly(T[i], [seq(Z[j],j=n-i+1..n,1)], [seq(T[j],j=1..i-1,1)]); od; return rT; end proc: construct_triangular_set_modp := proc(T, Z, p) # Input :: T := A triangular set with variables in Z. # Z := A list of variables. # p := a prime number. # # It is assumed that T is of the form t1(Z[-1]), t2(Z[-1],Z[-2]), etc, and that # the polynomials in T are in Maple's format. # # This method also turns T into a reduced triangular set. local rT,i,n; n := nops(Z); rT := [rpoly(T[1],Z[-1],p), seq(1..n-1)]; for i from 2 to nops(Z) do rT[i] := rpoly(T[i], [seq(Z[j],j=n-i+1..n,1)], p, [seq(T[j],j=1..i-1,1)]); od; return rT; end proc: is_radical := proc(T, p) # It is assumed that T is of the form t1(Z[-1]), t2(Z[-1],Z[-2]), etc, and that # the polynomials in T are in the recden format and not reduced modulo p. # # If a zero-divisor is encountered, pass the error through. local g,t,tp; for t in T do tp := phirpoly(t,p); g := RECDEN:-gcdrpoly(tp, diffrpoly(tp, 1)); if degrpoly(g,1) <> 0 then return false; fi; od; return true; end proc: isBad := proc(p,a::POLYNOMIAL,b::POLYNOMIAL, T) local t; # Test if p is a bad prime: see Definition 4 in paper for t in T do if iszerorpoly( phirpoly(lcrpoly(t),p) ) then return true fi; od; if idenomrpoly(a) mod p = 0 then return true; fi; if idenomrpoly(b) mod p = 0 then return true; fi; if iszerorpoly( phirpoly(lcrpoly(a),p) ) then return true fi; if iszerorpoly( phirpoly(lcrpoly(b),p) ) then return true fi; return false; end: makereduced := proc(T) local n,T2,i; n := nops(T); T2 := [seq(retextsrpoly(T[i]), i=1..n)]; for i from 2 to n do T2[i] := reduce_by_tset(T2[i], T2[1..i-1]); od; return T2; end proc: cgcd := proc(a,b,initial_prime:=11) local T,p,P,a_p,b_p,G_p,G,H,k,R,u,v,a0,b0,a_1,a_2,b_1,b_2,j,Z,n,M,T1,T2,i,N,counter,powoftwo,t,lc_u,st; global NPRIMES,DIVTIME; Z := getvars(a); T := getalgexts(a); # T is the algebraic extensions; polynomials are # primitive over Z n := nops(T); a0 := retextsrpoly(a); # a0 and b0 are the base polynomials of a and b. b0 := retextsrpoly(b); p := initial_prime - 1; G := rpoly(0,Z); P := 1; H := FAIL; k := 1; counter := 1; powoftwo := 1; while true do # Find a good, nonradical prime. N := false; while not N = true do p := nextprime(p); while isBad(p,a0,b0,T) do printf(" cgcd: %d is bad\n",p); p := nextprime(p); od; # Decide if p is a radical prime. If not, go to next prime. N := traperror(is_radical(T,p)); if N[1] = "inverse does not exist" then while N[1] = "inverse does not exist" do # See if the zero-divisor lifts. u := liftrpoly(N[2]); u := scarpoly(1/ilcrpoly(u),u); # make u monic j := nops(getvars(u)); t := T[j]; t := scarpoly(1/ilcrpoly(t),t); # make t monic v := quorpoly(phirpoly(t,p), u); M := traperror(hensellift_univar(u, v, t, p)); if M = FAIL or M[1] = "inverse does not exist" then N := M; # Pick a new prime if lifting fails, or try # lifting a new zero-divisor. else # Recursively call with new triangular sets. M := [M]; T1 := T; T2 := T; M[1] := scarpoly(idenomrpoly(M[1]),M[1]); M[2] := scarpoly(idenomrpoly(M[2]),M[2]); T1[j] := M[1]; T2[j] := M[2]; T1 := makereduced(T1); T2 := makereduced(T2); a_1 := reduce_by_tset(a0, T1); a_2 := reduce_by_tset(a0, T2); b_1 := reduce_by_tset(b0, T1); b_2 := reduce_by_tset(b0, T2); return cgcd(a_1, b_1, initial_prime), cgcd(a_2, b_2, initial_prime); fi; od; elif N = false then printf(" cgcd: %d is nonradical\n",p); # Pick a new prime. fi; od; #printf(" cgcd: trying prime p=%d\n",p); a_p := phirpoly(a,p); b_p := phirpoly(b,p); G_p := traperror(RECDEN:-gcdrpoly(a_p,b_p)); if G_p[1] = "inverse does not exist" then while G_p[1] = "inverse does not exist" do # See if the zero-divisor lifts. u := liftrpoly(G_p[2]); u := scarpoly(1/ilcrpoly(u),u); # make u monic j := nops(getvars(u)); t := T[j]; t := scarpoly(1/ilcrpoly(t),t); # make t monic v := quorpoly(phirpoly(t,p), u); M := traperror(hensellift_univar(u,v,t,p)); if M = FAIL or M[1] = "inverse does not exist" then G_p := M; # Pick a new prime. else # Recursively call with new triangular sets. M := [M]; M[1] := scarpoly(1/idenomrpoly(M[1]),M[1]); M[2] := scarpoly(1/idenomrpoly(M[2]),M[2]); T1 := T; T2 := T; T1[j] := M[1]; T2[j] := M[2]; T1 := makereduced(T1); T2 := makereduced(T2); a_1 := reduce_by_tset(a0, T1); a_2 := reduce_by_tset(a0, T2); b_1 := reduce_by_tset(b0, T1); b_2 := reduce_by_tset(b0, T2); return cgcd(a_1, b_1, initial_prime), cgcd(a_2, b_2, initial_prime); fi; od; else # Lift G_p over entire triangular set. G_p := retextsrpoly(G_p); if iszerorpoly(G) then G := G_p; # Initialization of gcd(a,b). NPRIMES := 1; else if degrpoly(G_p,1) = degrpoly(G,1) then # Current and past primes were lucky. G := modsrpoly(ichremrpoly([G,G_p])); P := P*p; NPRIMES := NPRIMES + 1; # Attempt trial division if its a power of 2. if k = powoftwo then H := irrrpoly(G); if H <> FAIL then H := reduce_by_tset(H, T); st := time(); if divrpoly(a,H) and divrpoly(b,H) then DIVTIME := time()-st; return H; fi; fi; powoftwo := powoftwo*2; fi; k := k + 1; elif degrpoly(G_p,1) < degrpoly(G,1) then printf(" cgcd: all past primes were unlucky\n"); G := modsrpoly(G_p); NPRIMES := 1; P := p; # Attempt trial division. k := 1; powoftwo := 1; H := irrrpoly(G); if H <> FAIL then H := reduce_by_tset(H, T); st := time(); if divrpoly(a,H) and divrpoly(b,H) then DIVTIME := time()-st; return H; fi; fi; else printf(" cgcd: %d is unlucky\n",p); fi; counter := counter + 1; # It shouldn't last more than 128 primes. This conditional is only used # while debugging. if counter > 256 then printf("I end up in here.\n"); return FAIL; fi; fi; fi; od; end proc: