/***
* This script was uploaded to https://software.mis.mpg.de/
* For an explanation of the functions and their usage, please check the website
* In particular, see the section on experimental results for various examples
**/
K<I> := CyclotomicField(4);
P2K  := ProjectiveSpace(K,2);
k    := Rationals();
P2<x,y,z> := ProjectiveSpace(k,2);
R := CoordinateRing(P2);
P3<x0,x1,x2,x3> := ProjectiveSpace(k,3);
/*
// A point configuration with one complex conjugated pair:

		Lpts := [ Eltseq(Vector(K, [-256545772*I + 313562678 , 354858009*I + 386629561 , 1])),
		          Eltseq(Vector(K, [256545772*I + 313562678 , -354858009*I + 386629561 , 1])),
  	    		  Eltseq(Vector(K, [-336705375 , -207070120 , 1])),
			        Eltseq(Vector(K, [301660661 , 363382661 , 1])),
			        Eltseq(Vector(K, [387324545 , -380342600 , 1])),
			        Eltseq(Vector(K, [-513184062 , 131906892 , 1])),
			        Eltseq(Vector(K, [-241165751 , 143263877 , 1])),
			        Eltseq(Vector(K, [-210847615 , 383053860 , 1]))];

		Lpts := [P2K ! pt : pt in Lpts];

		Lrealpts :=[ Eltseq(Vector(K, [-336705375 , -207070120 , 1])),
			     Eltseq(Vector(K, [301660661 , 363382661 , 1])),
			     Eltseq(Vector(K, [387324545 , -380342600 , 1])),
			     Eltseq(Vector(K, [-513184062 , 131906892 , 1])),
			     Eltseq(Vector(K, [-241165751 , 143263877 , 1])),
			     Eltseq(Vector(K, [-210847615 , 383053860 , 1]))];
*/

Lpts := [ [1, 0 , 0], [0 , 1 , 0], [0 , 0 , 1], [1 , 1 , 1], [10 , 11 , 1], [27 , 2 , 17], [-19 , 11 , -12], [-15 , -19 , 20] ];
Lpts := [P2 ! pt : pt in Lpts];

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

TRITGT_NOT_RAT := "Tritangent not defined over Rational Field.";


function CheckGeneralPosition(L)

	// We assume that L is a non-empty list of points in the projective plane.
	assert #L ne 0;
	assert Type(L[1]) eq Pt;
	P2 := Scheme(L[1]);
	assert IsProjective(P2) and IsAmbient(P2) and Dimension(P2) eq 2;

	// Check that the points are in general position
	k := BaseRing(P2);
	_<x,y,z> := PolynomialRing(k,3);
	pts := {ChangeUniverse(Eltseq(pt), k) : pt in L};

	// No coincident points
	if #L ne #pts then return false; end if;

	// 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]);
		if Determinant(M) eq 0 then return false; end if;
	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]);
		if Determinant(M) eq 0 then return false; end if;
	end for;

	// Any cubic through all 8 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);
		if Rank(MN) ne 10 then return false; end if;
	end for;

	return true;
end function;

// Given a projective space, n, and a list [<Pi, mi>] of points with multiplicites.
// Return the linear system of degree n homogenous polynomials vanishing at each Pi with multiplicity mi

function LinearSystemFromPoints(X,n,S)
	L := LinearSystem(X,n);
	for elt in S do
		L := LinearSystem(L,elt[1],elt[2]);
	end for;

	return L;
end function;

// INPUTS:
// Q     -- defining polynomial of the curve on the projective plane.
// u,v,w -- The polynomials used in constructing the curve.
function GetConeModel(Q,u,v,w)

	R<x,y,z> := Parent(u);
	k := BaseRing(Parent(u));
	A<T> := PolynomialRing(k);
	_<W> := PolynomialRing(A);

	L6Sections := [ u^i*v^(6-i) : i in [0 .. 6]] cat [ w*u^i*v^(4-i) : i in [0 .. 4]] cat [w^2*u^i*v^(2-i) : i in [0 .. 2]] cat [w^3];
	image_basis:= [ T^i : i in [0 .. 6]] cat [ W*T^i : i in [0 .. 4]] cat [W^2*T^i : i in [0 .. 2]] cat [W^3];

	mons := Monomials((x+y+z)^18);
	coeff_matrix := Matrix(k, [[MonomialCoefficient( s, mon) : mon in mons] : s in L6Sections]);

	Qsq := Q^2;
	vecQsq := Vector(k, [MonomialCoefficient( Qsq, mon) : mon in mons]);

	img_coeffs, N := Solution(coeff_matrix, vecQsq);
	assert Dimension(N) eq 0;

	return &+[ img_coeffs[i] * image_basis[i] : i in [1 .. #image_basis]];
end function;

// INPUTS:
// P2 -- ambient projective space. Must be P2 over a field K.
// Lpts -- Either
//		 i)  a list of eight points in P2(K).
//		 ii) a list of eight points in P2(K^sep). The set of eight points should
//		     define a cluster defined over K.
//
// OUTPUTS:
// C     -- A singular model of the branch curve in P2
// Q     -- The Defining Equation of C
// Qcone -- The defining equation of the branch curve in a cone P(1:1:2)
// u,v   -- The cubics generating the linear system of cubics through the eight points
// w     -- A sextic such that <u^2, uv, v^2, w> generates the sextics vanishing doubly
//	    at each point
function ConstructBranchCurve(P2, Lpts);

	assert IsProjective(P2) and IsAmbient(P2) and Dimension(P2) eq 2 and IsField(BaseRing(P2));
	R := CoordinateRing(P2);

	if Type(Lpts) eq SeqEnum then

		assert #Lpts eq 8 and &and[Type(pt) eq Pt : pt in Lpts];
		if not CheckGeneralPosition(Lpts) then error "Points not in general position."; end if;

		if &and[Scheme(pt) eq P2 : pt in Lpts] then
			L3 := LinearSystemFromPoints(P2, 3, [< pt, 1> : pt in Lpts]);
			L6 := LinearSystemFromPoints(P2, 6, [< pt, 2> : pt in Lpts]);
		else
			// This is the case our points define a cluster over K.
			calP := &join [Cluster(pt) : pt in Lpts];
			calP := Scheme(P2, [R ! f : f in DefiningEquations(calP)]);
			calPtimes2 := Scheme(P2, DefiningIdeal(calP)^2);

			L3 := LinearSystem(P2, 3);
			L3 := LinearSystem(L3, calP);
			L6 := LinearSystem(P2, 6);
			L6 := LinearSystem(L6, calPtimes2);
		end if;
	else
		error "ConstructBranchCurve requires a list of points as input.";
	end if;

	u := Sections(L3)[1];
	v := Sections(L3)[2];

	L3Squared := LinearSystem(P2, [u^2,u*v,v^2]);
	w := Sections(Complement(L6,L3Squared))[1];

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

	Qcone := GetConeModel(Q,u,v,w);

	return C, Q, Qcone, u,v,w;
end function;

// Function to recover the canonical model of the curve defined by
// Qcone in P3.
function GetCanonicalModelFromCone(P3, Qcone)
    k := BaseRing(P3);
    _<x0,x1,x2,x3> := CoordinateRing(P3);
    R<t> := PolynomialRing(k);
    coeffs := Coefficients(Qcone);
    mons := [R ! m : m in Monomials(Qcone)];
    q := x0*x2-x1^2;

    F := CoordinateRing(P3) ! (x1^6*(&+[Evaluate(coeffs[i], x0/x1)
					* Evaluate(mons[i], x0*x3/x1^2)
					: i in [1 .. #coeffs]]));
    X := Curve(Scheme(P3, [q, F]));

    // X is not defined properly since the "inverse" of P112 --> P3 is a bit delicate.

    comps := IrreducibleComponents(X);
    C := [co : co in comps | IsReduced(co)][1];
    return Curve(C);
end function;


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

// The tritangent planes of a space sextic.



function MakeMultiplicityMap(mults)
	multperms := { [mults[f[i]] : i in [1 .. #mults]] : f in Permutations({1 .. #mults}) };

	// Remove the multiplicites corresponding to the image under Bertini.
	for x in multperms do
		minusX := [2-y : y in x];
		// Note that destructive operations do not apply to the for loop condition, but they
		// do apply to the variable `multperms`. This is why we need the if statement.
		if x in multperms then Exclude(~multperms, minusX); end if;
	end for;

	function MultMap(Lpts)
		assert #Lpts eq #mults;
		return { [<Lpts[i], mp[i]> : i in [1 .. #Lpts]] : mp in multperms};
	end function;
	return MultMap;
end function;

multmap0 := MakeMultiplicityMap([1,0,0,0,0,0,0,0]);
multmap1 := MakeMultiplicityMap([1,1,0,0,0,0,0,0]);
multmap2 := MakeMultiplicityMap([1,1,1,1,1,0,0,0]);
multmap3 := MakeMultiplicityMap([2,1,1,1,1,1,1,0]);


// INPUTS:
// F     -- A polynomial whose zero locus in P2 is the union of an exceptional curve and its
//	    image under the Bertini involution.
// u,v,w -- The polynomials used in constructing the curve.
// A     -- The ring of definition.
function GetConeTritangent(F, u,v,w, A)
	R<x,y,z> := Parent(u);
	k := BaseRing(Parent(u));
	_<W> := PolynomialRing(A);
	T := A.1;

	L2Sections := [ u^i*v^(2-i) : i in [0 .. 2]] cat [w];
	image_basis:= [ T^i : i in [0 .. 2]] cat [ W ] ;

	mons := Monomials((x+y+z)^6);
	coeff_matrix := Matrix(k, [[MonomialCoefficient( s, mon) : mon in mons] : s in L2Sections]);

	vecf := Vector(k, [MonomialCoefficient( F, mon) : mon in mons]);

	img_coeffs, N := Solution(coeff_matrix, vecf);
	assert Dimension(N) eq 0;

	return &+[ img_coeffs[i] * image_basis[i] : i in [1 .. #image_basis]];
end function;





function GetConeCurve(P2, n, Qcone, S, u, v, w : descendToRationals:=false)
  if n ne 0 then
	  L := LinearSystemFromPoints(P2,n,S);
	  assert #Sections(L) eq 1;
	  f := Sections(L)[1];
	  minusS := [ <a[1], 2-a[2]> : a in S];
  else
    f := 1;
	  minusS := [ <a[1], 2+a[2]> : a in S];
  end if;

	L := LinearSystemFromPoints(P2, 6-n,minusS);
	assert #Sections(L) eq 1;
	g := Sections(L)[1];

	// The polynomial defining the exceptional curve times its Bertini conjugate.
	// We normalize to ensure that at least one coefficient of F is rational.
	F := f*g; F := F/LeadingCoefficient(F);

	if descendToRationals then
		k  := Rationals();
		P2 := ChangeRing(P2, k);

		// TODO: Enforce that the representative inputs are defined over k.

		try
			F  := Parent(u) ! ChangeRing(F, k);
		catch e
			// In the case that F is not defined over Q, we know that the tritangent is not
			// defined over R. In this case, we return "false" along with a dummy
			// discriminant value.
			return false, TRITGT_NOT_RAT;
		end try;
	end if;

	return GetConeTritangent(F, u, v, w, BaseRing(Parent(Qcone)));
end function;


function AllTritangentsInP3()

    C, Q, Qcone, u,v,w := ConstructBranchCurve(P2,Lpts);

    Tritangents8 :=    [ GetConeCurve(P2, 0, Qcone, S, u, v, w) : S in multmap0(Lpts)];
    Tritangents28 :=   [ GetConeCurve(P2, 1, Qcone, S, u, v, w) : S in multmap1(Lpts)];
    Tritangents56 :=   [ GetConeCurve(P2, 2, Qcone, S, u, v, w) : S in multmap2(Lpts)];
    Tritangents56_2 := [ GetConeCurve(P2, 3, Qcone, S, u, v, w) : S in multmap3(Lpts)];
    Tgt_list := Tritangents8 cat Tritangents28 cat Tritangents56 cat Tritangents56_2;


    // Put tritangent equations and the curve in P3
    P3<x0,x1,x2,x3> := ProjectiveSpace(k,3);
    Tritangents := [];
    for tt in Tgt_list do
    cox3 := Coefficient(Coefficient(tt,1),0);
    cox2 := Coefficient(Coefficient(tt,0),0);
    cox1 := Coefficient(Coefficient(tt,0),1);
    cox0 := Coefficient(Coefficient(tt,0),2);
    ttcoeffs := [cox0,cox1,cox2,cox3];
    tgteq := &+[P3.i*ttcoeffs[i] : i in [1 .. 4]];
    Include(~Tritangents, tgteq);
  end for;
  return Tritangents;
end function;

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

//Totally real check:

// INPUTS:
// P2    -- ambient projective space. Must be P2
// Qcone     -- defining polynomial of the curve on the quadric cone.
// F     -- A polynomial whose zero locus in P2 is the union of an exceptional curve and its
//	    image under the Bertini involution.
// u,v,w -- The polynomials used in constructing the curve.
function OptimizedCheckCone(P2, Qcone, F , u, v, w )

	ell := GetConeTritangent(F, u, v, w, BaseRing(Parent(Qcone)));
	assert Degree(ell) eq 1;

	g := BaseRing(Parent(Qcone)) ! Evaluate(Qcone, -Coefficient(ell,0)/Coefficient(ell,1));

	// g should be a square up to leading coefficient, so we want to take the square root.
	h := GCD(g, Derivative(g));
	assert Degree(ExactQuotient(g, h^2)) eq 0;
	assert LeadingCoefficient(h) eq 1;

	if Degree(h) lt 1 then return true; end if;
	assert Degree(h) lt 4;
return Discriminant(h) ge 0, Discriminant(h), ell;
end function;


// INPUTS:
// X -- ambient projective space. Must be P2
// n -- degree of plane curve corresponding to exceptional curve
// Qcone -- defining polynomial of the curve on the quadric cone.
// S -- List of points with multiplicity defining the exceptional curve
// u,v,w -- The polynomials used in constructing the curve.
// PARAMETERS:
// descendToRationals -- determines whether to descend objects to QQ.

function CheckConeCurve(P2, n, Qcone, S, u, v, w : descendToRationals:=false)
	L := LinearSystemFromPoints(P2,n,S);
	assert #Sections(L) eq 1;
	f := Sections(L)[1];

	minusS := [ <a[1], 2-a[2]> : a in S];
	L := LinearSystemFromPoints(P2, 6-n,minusS);
	assert #Sections(L) eq 1;
	g := Sections(L)[1];

	// The polynomial defining the exceptional curve times its Bertini conjugate.
	// We normalize to ensure that at least one coefficient of F is rational.
	F := f*g; F := F/LeadingCoefficient(F);

	if descendToRationals then
		k  := Rationals();
		P2 := ChangeRing(P2, k);

		// TODO: Enforce that the representative inputs are defined over k.

		try
			F  := Parent(u) ! ChangeRing(F, k);
		catch e
			// In the case that F is not defined over Q, we know that the tritangent is not
			// defined over R. In this case, we return "false" along with a dummy
			// discriminant value.
			return false, TRITGT_NOT_RAT;
		end try;
	end if;

	return OptimizedCheckCone(P2, Qcone, F, u, v, w );
end function;

// INPUTS:
// C     -- The model of the curve on the projective plane.
// P     -- A point on the plane model. One of the eight original points chosen in
//	    the construction
// PARAMETERS:
// assertReal -- determines whether to check if the base ring is a finite extension of Q.

function IsTotallyRealTangentCone(C,P : assertReal:=true)
	// We assume that P does not lie on the hyperplane at infinity.
//	assert P[3] ne 0;
	assert Scheme(P) eq C;

	P2<x,y,z> := AmbientSpace(C);
	k := BaseRing(P2);
	tc := TangentCone(C,P);

	assert Degree(tc) eq 3;

	if assertReal then
		assert (k eq Rationals() or IsNumberField(k));
	end if;

	if not IsOrdinarySingularity(P) then
		// The only way this happens is if some of the contact points with C coincide.
		// In particular, they must all be real
		return true, 0;
	end if;

	F := Evaluate(DefiningEquation(tc), [x,y,0]);

	// Define the number field over which the tangent directions split.
	R<t> := PolynomialRing(BaseRing(P2));

	if Coefficient(F, x , 3) ne 0 then
		f := Evaluate(F, [t,1,0]);
	else
		f := Evaluate(F, [1,t,0]);
	end if;

	return (Discriminant(f) ge 0), Discriminant(f);
end function;

function TotRealTritCount(P2,P2K,Lpts,Lrealpts)

    C, Q, Qcone, u,v,w := ConstructBranchCurve(P2,Lpts);

    //check if each tritangent corresponding to blow ups at points is real
    if Lrealpts ne [] then
	Lrealpts := [P2 ! pt : pt in Lrealpts];
	pointconds:= [IsTotallyRealTangentCone(C, C ! pt) : pt in Lrealpts];
    else
        pointconds := [];
    end if;


    // Check if each of the lines gives rise to totally real tritangent
    lineconds := [CheckConeCurve(P2K, 1, Qcone, S, u, v, w :  descendToRationals:=true) : S in multmap1(Lpts)];
    //    --    each of the conics gives rise to totally real tritangent

    conicconds:= [CheckConeCurve(P2K, 2, Qcone, S, u, v, w :  descendToRationals:=true) : S in multmap2(Lpts)];
    //    --    each of the cubics gives rise to totally real tritangent

    cubicconds:= [CheckConeCurve(P2K, 3, Qcone, S, u, v, w :  descendToRationals:=true) : S in multmap3(Lpts)];

    // Count the number of totally real tritangents

    numtotsreal :=   #[l : l in pointconds | l eq true]
		     + #[l : l in lineconds | l eq true]
		     + #[l : l in conicconds | l eq true]
		     + #[l : l in cubicconds | l eq true];

    return numtotsreal;

end function;
