
// Magma code accompanying the paper
// "An explict family of cubic number fields with large 2-rank of the class group"
//
// Avinash Kulkarni, 2016-09-09


function GetSmallBasis(M)
	Mat := Matrix(M);
	L := Lattice(Mat);
	StdLat := StandardLattice(NumberOfColumns(Mat));

	return ChangeRing(BasisMatrix(LLL(StdLat meet L)), CoefficientRing(Mat));
end function;

function CheckGeneralPosition(L)
	
	k := Integers();
	_<x,y,z> := PolynomialRing(k,3);
	pts := {ChangeUniverse(pt cat [1],k) : pt in L};

	cons := {};
	
	// No three lie on a line
	for S in Subsets(pts,3) do
		mons := Monomials((x+y+z));
		M := Matrix(k,3,3,[[Evaluate(m,pt) : pt in S] : m in mons]);
		Include(~cons, Determinant(M));
	end for;
	
	// No six lie on a conic
	for S in Subsets(pts,6) do
		mons := Monomials((x+y+z)^2);
		M := Matrix(k,6,6,[[Evaluate(m,pt) : pt in S] : m in mons]);
		Include(~cons, Determinant(M));
	end for;
	
	// Any cubic through all 3 is non-singular at each point
	for s in pts do
		mons := Monomials((x+y+z)^3);
		M := Matrix(k,10,8,[[Evaluate(m,pt) : pt in pts] : m in mons]);
		N := Matrix(k,10,3,[[Evaluate(Derivative(mon,1),s),Evaluate(Derivative(mon,2),s),Evaluate(Derivative(mon,3),s)] : mon in mons]);

		MN := HorizontalJoin(M,N); 
 
		Include(~cons, GCD(Minors(MN,10)));
	end for;
	
	if 0 in cons then
		return false,cons;
	end if;
	return true, cons;
end function;

function b(n)
	return [n^2 - 1, n^3 - 1];
end function;



k:=Rationals();

S<x,y> := PolynomialRing(k,2);

Lmon6 := Monomials((x+y+1)^6);
Lmon3 := Monomials((x+y+1)^3);

Lmon3 := [mon : mon in Lmon3 | TotalDegree(mon) gt 0];
Lmon6Node:= [mon : mon in Lmon6 | TotalDegree(mon) gt 1];

//Note by the choice of monomials we have fixed [0:0:1] to be in our collection of base points. The function b ensures that all 8 points lie on the cuspidal cubic parametrized by t -> (t^2-1,t^3-1).

Lfix := [b(-1),b(-2),b(2),b(3),b(4),b(5),b(7)];
Lpts := Lfix;

boo, strcons := CheckGeneralPosition(Lpts cat [b(1)]);
assert boo;

eqns3 := [[Evaluate(mon,P) : mon in Lmon3] : P in Lpts];

eqns6 := [[Evaluate(mon,P) : mon in Lmon6Node] : P in Lpts]
			cat [[Evaluate(Derivative(mon,1),P) : mon in Lmon6Node] : P in Lpts]
				cat [[Evaluate(Derivative(mon,2),P) : mon in Lmon6Node] : P in Lpts];

				
//The basis for L(-KX) should be computed first.

M3 := Matrix(k,#eqns3, #Lmon3,eqns3);
ker:= NullspaceOfTranspose(M3);
B := GetSmallBasis(Basis(ker));

u:= Polynomial(Eltseq(B[1]), Lmon3);
v:= Polynomial(Eltseq(B[2]), Lmon3);

//Solve eqns6 to get the 4-dimensional space of solutions for L(-2KX)

M6 := Matrix(k,#eqns6, #Lmon6Node,eqns6);
ker6 := NullspaceOfTranspose(M6);
B := GetSmallBasis(Basis(ker6));
assert Dimension(ker6) eq 4;

//Ensure that w is linearly independent from u^2, uv, v^2. 

A := Matrix([[MonomialCoefficient(f,mon): mon in Lmon6Node] : f in [u^2,u*v,v^2]]);
wec := B[1];

for j in [1 .. NumberOfRows(B)] do
	wec := B[j];
	if Rank(VerticalJoin(A,Matrix(wec))) eq 4 then
		break;
	end if;
end for;

w := Polynomial(Eltseq(wec), Lmon6Node);

//We now construct the curve.
R<x,y,z> := PolynomialRing(k,3);
FracR := FieldOfFractions(R);

function Homogenize(f)
	return R ! z^(TotalDegree(f))*Evaluate(f,[x/z,y/z]);
end function;

u := R ! Homogenize(u);
v := R ! Homogenize(v);
w := R ! Homogenize(w);

J := JacobianMatrix([u,v,w]);
Q := Determinant(J);
						
P2 := ProjectiveSpace(k,2);
C := Curve(P2,Q);


// We construct a morphism to an affine model of the form 0 = c0 * w^3 + ... + c6(t). 
A2<t,W> := AffineSpace(k,2);
phi := map<C->A2| [u/v, w/v^2]>;
Caff := Curve(Image(phi));
phi := Restriction(phi,C,Caff);

F := DefiningEquation(Caff);

// A monic version of the defining equation.
F := F/Coefficient(F,W,3);


//************************************************************************************************************

// Computing the bad primes

//************************************************************************************************************

// Every prime for which the reduction of the base points remains in general position is good. Thus the bad primes are those dividing at least one of the determinants in the general position conditions.

BadPrimes := Factorization(LCM(strcons));

P112 := ProjectiveSpace(k,[1,1,2]);
kP112 := FieldOfFractions(CoordinateRing(P112));
Fhom := Numerator(Evaluate(F,[P112.1/P112.2,P112.3/P112.2^2]));
Fhom := LCM([Denominator(coeff) : coeff in Coefficients(Fhom)])*Fhom;
Czz := Curve(P112,Fhom);


// Correctness depends on the reliability of Magma's Regular Model functions
for bp in BadPrimes do
	p := bp[1];
	if p lt 5 then // We avoid 2,3 since the reduction type is likely horrible.
	   continue;
	end if;
	try
	   Cp := RegularModel(Czz,p);
	   CompGrp := ComponentGroup(Cp);
	   if IsOdd(#CompGrp) then
 		Exclude(~BadPrimes,bp);
	   end if;
	catch e;
	end try;
end for;

BadPrimes := [bp[1] : bp in BadPrimes];

// We insert 2,3,5 into the bad primes if they are not already there.

for p in [5,3,2] do
	if Index(p,BadPrimes) eq 0 then
		Insert(~BadPrimes,1,p);
	end if;
end for;

//************************************************************************************************************

// Computing the explict congruence data

//************************************************************************************************************

// We recall our 8 constructed base points and determine the 9th point defining the identity section.

PencilBasePoints := RationalPoints(Scheme(P2,u) meet Scheme(P2,v)); 
DelPezzoBase := {P2 ! [pt[1],pt[2],1] : pt in Lpts cat [[0,0]]};
idE := Random(PencilBasePoints diff DelPezzoBase);  

// Verify that MAGMA's generators for the function field correspond to t and W.
kC:=FunctionField(C);    
assert Evaluate(Q,[kC.1,kC.2,1]) eq 0;


// This function attempts to find a unique effective degree 3 place in the decomposition of a given divisor.
// An error occurs otherwise.
function getRightPlace(D)
    supp := Decomposition(D);  

    // Checks through all possible decompositions.
    Dplus := {dv : dv in supp | dv[2] gt 0 and Degree(dv[1]) lt 4}; 
    summands := &cat [ [ a[1] : i in [1 .. a[2]] ] : a in Dplus ];   
    I := {1 .. #summands};
    singles := { s : s in summands | Degree(s) eq 3};
    pairs := { summands[ss[1]] + summands[ss[2]] where ss := Setseq(s) : s in Subsets(I,2) | 
			Degree(summands[ss[1]] + summands[ss[2]]) eq 3 where ss := Setseq(s)};
    triples := { summands[ss[1]] + summands[ss[2]] + summands[ss[3]] where ss := Setseq(s) : s in Subsets(I,3) | 
			Degree(summands[ss[1]] + summands[ss[2]] + summands[ss[3]]) eq 3 where ss := Setseq(s)};

    candidates := singles join pairs join triples;
    if #candidates ne 1 then
	error "No unique effective degree 3 place. Case not implemented.";
    end if;

   return Random(candidates);
end function;


// Since the identity section corresponds to the vanishing theta characteristic we can easily compute its divisor.
dv2:=Divisor(Caff,1/t);
Osect:= getRightPlace(dv2);

// Given one of the base points on the plane model, compute the Weil pairing function for the corresponding [2]-torsion point.
function makeDescentFunction(bp)
    plcs:=Places(C ! bp);
    plc:=plcs[1];
    uv_kC:=Evaluate((u/v),[kC.1,kC.2,1]);
    val:=Evaluate(uv_kC,plc);
    P1 := getRightPlace(Divisor(Caff,Evaluate(MinimalPolynomial(val),t)));

    bl,f12:=IsPrincipal(2*(P1-Osect));
    return f12;
end function;


function findLocalConditions(f12)	
    // Here, one should compute the Taylor expansions of the numerator and denominator about the cusp point. By design the cusp point is [0,0] so the Taylor expansion is already computed.

    denF12 := Denominator(f12);
    numF12 := Numerator(f12);

    // We need the numerator and denominator to have integer coefficients.
    d1 := LCM([Denominator(x) : x in Coefficients(denF12)]);
    d2 := LCM([Denominator(x) : x in Coefficients(numF12)]);
    d := LCM(d1,d2);

    denF12 := d*denF12;
    numF12 := d*numF12;

    assert CanChangeUniverse(Coefficients(denF12),Integers()) and CanChangeUniverse(Coefficients(numF12),Integers());

    a0 := Evaluate(numF12,[0,0]);
    b0 := Evaluate(denF12,[0,0]);

    // Denote || to be the p-adic absolute value. Since the coefficients of the polynomials both denF12 and numF12 are integers, if l := min{ |a0|,|b0|} and |P-P0| < l, then
    // f12(P) = f12(P0) mod p^l. In particular the valuation of f12(P)/f12(P0) will be 1.

	if Degree(denF12) eq 0 then
	Lambdas := [Valuation(a0,p) + 1 : p in BadPrimes];	
	else
	Lambdas := [Max(Valuation(a0,p),Valuation(b0,p)) + 1 : p in BadPrimes];
	end if;

    // As in the proof of Lemma 3.4, we need to ensure that the coefficients of the defining form of Caff evaluated at p are small enough. We can determine a sufficient condition using only the coefficients of Caff.

    coeff_error := LCM([Denominator(c) : c in Coefficients(F)]);    
    nbhd := coeff_error*&*[BadPrimes[i]^(Lambdas[i]) : i in [1 .. #BadPrimes]];
    return nbhd, Lambdas;
end function;


descFuns := [makeDescentFunction(bp) : bp in DelPezzoBase];
Neighbourhood := LCM([ findLocalConditions(f) : f in descFuns]);







