////////////////////////////////////////////////////////////////////////
// Couveignes's zone
// Based on his article in the most famous LNM
// TODO: make it faster?
////////////////////////////////////////////////////////////////////////
// Assumes deg(f) odd and that prodalg is the square of an element of
// Z[X]/(f(X)).
// max(|a|, |b|) <= U.
JMCBound := function(f, m, U, cardS)
local d, i, Nf;
d:=Degree(f);
Nf:=Sqrt(&+ [ Coefficient(f, i)^2 : i in [0..d] ]);
bd:=&+ [ (m/Nf)^i : i in [0..d] ];
return d^(3/2)*Nf^d*(2*U*Nf)^(cardS div 2)*bd;
end function;
JMCUseComb := function(N, f, m, alg, maxp, prodrat, U, cardS, cardS0, bound)
d:=Degree(f);
cd:=Coefficient(f, d);
x:=Parent(f).1;
g:=PolynomialRing(Integers()) ! Eltseq(cd^(d-1)*Evaluate(f, x/cd));
printf "g = %o\n", g;
cd:=Integers() ! cd;
prodmi:=1; Bmmi:=0;
mi:=Max(maxp, 10^3);
NFSPRINTLEVEL:=1;
tmi:=[];
txi:=[];
nch:=0;
while true do
nch:=nch+1;
// cd <> 0 mod p??
repeat
mi:=NextPrime(mi);
Fmi<Z>:=PolynomialRing(GF(mi));
until IsIrreducible(Fmi ! Eltseq(f));
if NFSPRINTLEVEL ge 1 then
if (nch mod 100) eq 0 then
printf "ch[%o] = %o at %o (%o)\n",nch,mi,Cputime(),1.0*prodmi;
end if;
end if;
tp:=Cputime();
// prodalgmi <- g'(omega)^2 mod mi
// prodalgmi:=Rem(diff(g, X)^2, g, X) mod mi;
gmi:=Fmi ! Eltseq(g);
Fq<W>:=ext<GF(mi) | gmi>;
FqV<V>:=PolynomialRing(Fq);
prodalgmi:=Fq ! Eltseq((Derivative(g)^2 mod g));
Nbeta:=1;
for i:=1 to #alg do
ab:=Split(alg[i]);
ab:=StringToIntegerSequence(ab[1]);
if (ab[1] eq 0) and (ab[2] eq 0) then ind:=i+1; break; end if;
prodalgmi:=prodalgmi*(cd*ab[1]-ab[2]*W);
end for;
if NFSPRINTLEVEL ge 2 then printf "prod = %o\n", prodalgmi; end if;
// Nbeta <- +/- N(g'(omega)) * cd^((d-1)*cardS/2) * prodlpalpg
Nbeta:=Nbeta * Resultant(Derivative(gmi), gmi);
Nbeta:=Nbeta * Modexp(cd, (d-1)*(cardS div 2), mi);
// take care to free relations
Nbeta:=Nbeta * Modexp(cd, cardS0 div 2, mi);
for i:=ind to #alg do
ab:=Split(alg[i]);
ab:=StringToIntegerSequence(ab[1]);
if (ab[1] eq 0) and (ab[2] eq 0) then break; end if;
// ab = [p, r, e]
Nbeta:=Nbeta * Modexp(Abs(ab[1]), ab[3], mi);
end for;
// do not forget the sign problem (really???)
// if (ex[1+kchar+1+krat+1]/2) mod 2 = 1 then
// Nbeta:=-Nbeta mod mi;
// end if;
betami:=Roots(V^2-prodalgmi);
betami:=betami[1][1]; // in Fq
// could be Resultant(betami, g, X) mod mi;
// but take care to the fact that if deg(betami) <> d-1
// there is a sign problem... Try "Ex2();" to be convinced.
Nbetami:=betami ^ ((mi^d-1) div (mi-1));
if Nbeta ne Nbetami then
betami:=-betami;
Nbetami:=-Nbetami;
end if;
if NFSPRINTLEVEL ge 2 then
printf "B(X) = %o mod %o\n", betami, mi;
end if;
if Nbeta ne Nbetami then
error("Nbetami <> +/- Nbeta mod mi");
end if;
// Bmmi:=numtheory[mcombine](prodmi,Bmmi,mi,subs(X=cd*m, betami) mod mi);
tmp:=PolynomialRing(Integers()) ! Eltseq(betami);
tmp:=Evaluate(tmp, cd*m);
Bmmi:=CRT([Bmmi, tmp], [prodmi, mi]);
prodmi:=prodmi * mi;
if NFSPRINTLEVEL ge 2 then
printf "gamma mod mi = %o => beta mod mi = %o",prodalgmi,betami;
printf " (computed in %o)\n", Cputime(tp);
printf "prodmi = %o, Bm = %o\n", prodmi, Bmmi;
end if;
// tmi[nch]:=mi;
// txi[nch]:=subs(X=cd*m, betami) mod mi;
//
// phisqrtalg:=mods(Bmmi, prodmi) mod N;
if 2*Bmmi gt prodmi then
phisqrtalg:=-Bmmi mod N;
else
phisqrtalg:=Bmmi mod N;
end if;
// we should have phisqrtalg^2 = prodrat^2 mod N
if (phisqrtalg^2-prodrat^2) mod N ne 0 then
/// printf "phi(sqrtalg)^2 != prodrat^2 mod N\n";
if prodmi gt bound then
error("?? prodmi > bound\n");
end if;
else
printf "// nch = %o\n", nch;
fact:=GCD(phisqrtalg-prodrat, N);
printf "--> Igcd = [%o, %o]\n", fact, N div fact;
printf "[%o, %o]\n", IsProbablePrime(fact), IsProbablePrime(N div fact);
// Bm:=JMCCRT(N, tmi, txi, nch);
// if Bm <> phisqrtalg then
// ERROR("Bm <> phisqrtalg");
// end if;
break;
end if;
end while;
return fact;
end function;
function finishWithJMC(N, f, m, ratdep, algdep)
d:=Degree(f);
cd:=Integers() ! Coefficient(f, d);
alg:=Split(Read(algdep));
// print alg;
cardS:=0;
cardS0:=0;
U:=0;
for i:=1 to #alg do
ab:=Split(alg[i]);
ab:=StringToIntegerSequence(ab[1]);
if (ab[1] eq 0) and (ab[2] eq 0) then
ind:=i+1;
break;
end if;
cardS:=cardS+1;
if ab[2] eq 0 then
// free relation (a, 0)
cardS0:=cardS0+1;
end if;
if Abs(ab[1]) gt U then U:=Abs(ab[1]); end if;
if Abs(ab[2]) gt U then U:=Abs(ab[2]); end if;
end for;
maxp:=0;
for i:=ind to #alg do
ab:=Split(alg[i]);
ab:=StringToIntegerSequence(ab[1]);
if ab[1] gt maxp then maxp:=ab[1]; end if;
end for;
printf "cardS = %o, cardS0 = %o, U = %o", cardS, cardS0, U;
printf ", maxp = %o\n", maxp;
if (cardS mod 2) eq 1 then
error("cardS is odd");
end if;
if (cardS0 mod 2) eq 1 then
error("cardS0 is odd");
end if;
prodrat:=Split(Read(ratdep));
prodrat:=StringToInteger(prodrat[1]);
df:=Derivative(f);
prodrat:=(Integers()!(prodrat * Evaluate(df, m))) mod N;
prodrat:=(prodrat * (Modexp(cd, d-2+(cardS div 2), N))) mod N;
// print prodrat;
bound:=JMCBound(f, m, U, cardS);
printf "Bound = %o\n", bound;
return JMCUseComb(N, f, m, alg, maxp, prodrat, U, cardS, cardS0,bound);
end function;
Qx<x>:=PolynomialRing(Rationals());
Try := procedure(N, f, m, str)
rat:="c20.dep.rat." * str;
alg:="c20.dep.alg." * str;
finishWithJMC(N, f, m, rat, alg);
end procedure;