////////////////////////////////////////////////////////////////////////
// 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;