indred.magma

// Infer the structure of coefficients in the sum by computing P+nP 
// for reduced Weierstrass form

SetAutoColumns(false);
SetColumns(0);

// Formal addition for reduced form
function formalAddReduced(P, Q)
	Rk := Parent(P[1]);
	
	A := Rk.1; B := Rk.2;
	X := Rk.3;

	X1 := Rk!P[1]; Y1 := Rk!P[2]; Z1 := Rk!P[3];
	X2 := Rk!Q[1]; Y2 := Rk!Q[2]; Z2 := Rk!Q[3];

	g1 := X2*Y1 + X1*Y2;
	g2 := Y2*Z1 + Y1*Z2;

	H1 := - X2*Z1*A - X1*Z2*A - 3*Z1*Z2*B + Y1*Y2;
	H2 := Z1*Z2*A^2 - X1*X2*A - 3*X2*Z1*B - 3*X1*Z2*B;
	H3 := Z1*Z2*A + 3*X1*X2;
	H4 := X2*Z1*A + X1*Z2*A + 3*Z1*Z2*B + Y1*Y2;

	X3 := g1*H1 + g2*H2;
	Y3 := H1*H4 - H2*H3;
	Z3 := g1*H3 + g2*H4;
	
	Q := [Rk!X3, Rk!Y3, Rk!Z3];
	
	// Go back to standard form
	while Q[2] ne Rk!1 do
		coeff := Coefficients(Q[2], X);
		for j in [2..#coeff] do
			if coeff[j] ne 0 then
				inv_y := 1 + Rk!X^(j-1) * (-1) * Rk!coeff[j];
				Q[1] *:= inv_y;
				Q[2] *:= inv_y;
				Q[3] *:= inv_y;
				break;
			end if;
		end for;
	end while;
	return Q;
end function;

function fit(ys)
	// ======================
	// ys := y corresponding to [1..#ys]
	//
	// Returns: coeffs of polynomial
	// ======================
	//print "xs:",xs;
	//print "ys:",ys;
	Q := RationalField();
	P<x> := PolynomialRing(Q,1);
	f := Interpolation([Q!i : i in [1..#ys]], ys);
	return f;
end function;

k := 1331;
proof := false; // Compute the induction step

// We pick a generic curve
R<A, B, X, n> := PolynomialRing(Rationals(), 4);

// Z = F(X)
load "zfxred_stored_300.magma";
print "Loaded";

n_coeff := k+1; // Number of coefficients to determine (up to X^n_coeff)

// Closed form for X for the first few coefficients
psi1 := n;
formal_X := psi1*X; 

// Try to obtain the structure of X for each value of k
for k in [2..n_coeff+1] do 
	
	print "=======\nk = ", k, "\n=======";
	if (k-1) mod 2 eq 0 then
		SetOutputFile("psired_stored.magma");
		printf "psi%o := 0;\n", k-1;
		UnsetOutputFile();
		print "Even k --> skip";
		continue;
	end if;
	I := ideal<R | R.3^k>; // R.6 := X
	Rk := R/I;

	// Casting of variables
	A := Rk!R.1;
	B := Rk!R.2;
	X := Rk!R.3;
	n := Rk!R.4;

	// Casting of F
	Fk := Rk!F;
	
	// Points
	// P : (X:1:f(X))
	P := [Rk!X, Rk!1, Fk];
	Q := P;
	
	residuals := [Rk!0];

	for i in [2..k+10] do 		
		Q := formalAddReduced(P, Q);	
		//print "Q: ", Q[1];
		res := Q[1] - Rk!Evaluate(formal_X, n, i); // nP - n(formalP) where formalP is built with psi_i up to k-1
		Append(~residuals, res);
		
		if res eq 0 then continue; end if;

		cff := Coefficients(res, X);
		assert IsDivisibleBy(res, Rk!X^(k-1)); // Check correctness of the previous psi_i
	end for;
	
	if residuals[#residuals] eq Rk!0 then
		print "Empty residuals --> skip";
		SetOutputFile("psired_stored.magma");
		printf "psi%o := %o;\n", k-1, Coefficient(formal_X, X, k-1);
		UnsetOutputFile();
		continue;
	end if;

	residuals := [Coefficient(r, Rk!X, k-1) : r in residuals];
	//print residuals;

	all_monomials := Monomials(residuals[#residuals]);
	final_coeff := 0;
	for m in all_monomials do
		//printf "-----\nFitting poly for %o...\n", m;
		
		ys := [MonomialCoefficient(residuals[i], m) : i in [1..k+5]];
		
		// Check on last values
		fitted_poly := fit(ys);
		for i in [k+5..k+10] do
			assert  Evaluate(fitted_poly, i) eq MonomialCoefficient(residuals[i], m);
		end for;
		//print "Check passed!";
		
		// Cast into R (R.4 = n)
		fitted_poly := MultivariatePolynomial(R, fitted_poly, R.4);
		//print "Fitted: ",fitted_poly;
		final_coeff +:= R!m * fitted_poly;

	end for;
	SetOutputFile("psired_stored.magma");
	printf "psi%o := %o;\n", k-1, final_coeff;
	//printf "final_coeff := %o\n", final_coeff;
	UnsetOutputFile();
	printf "Computed psi%o\n", k-1;

	//printf "##################################\nFinal coeff for X^%o: %o\n##################################\n", (k-1), final_coeff;

	// R.3 = X
	formal_X +:= final_coeff*(R.3)^(k-1);

	// Check induction
	if proof then
		// Basis Point
		P_1 := [Rk!Evaluate(formal_X, n, 1), Rk!1, Evaluate(Fk, X, Rk!Evaluate(formal_X, n, 1))];
		assert P_1 eq [Rk!X, Rk!1, Evaluate(F, X, Rk!X)];
		print "Induction base ok...";

		// Induction Step
		P_n := [Rk!formal_X, Rk!1, Rk!Evaluate(Fk, X, formal_X)];
		P_n1 := [Rk!Evaluate(formal_X, n, Rk!(n+1)), Rk!1, Rk!Evaluate(F, X, Rk!Evaluate(formal_X, n, Rk!(n+1)))];

		assert formalAddReduced(P_n, P_1) eq P_n1;
		print "Induction step ok!";
	end if;
	
end for;