A := SymmetricMatrix([3,1,4,1,5,9,2,6,5,3]);
B := SymmetricMatrix([2,7,1,8,2,8,1,8,2,8]);

P<x> := PolynomialRing(Integers());

function LS2Basis(S)
	// The interesting case is where the base field doesn't
	// have class number 1. But for the example in the MSS
	// paper it does have h = 1.

	Senum:=[t:t in S];
	K := FieldOfFractions(Order(Senum[1]));
	hK := ClassNumber(K);
	assert(hK eq 1);

	// printf "Base field is %o of class number %o\n",K,hK;

	// In the class-number-1 case, K(S,2) is generated by
	// the generators of the units of K and the elements
	// Ti such that Si = Ti OK -- it's literally just the
	// units union the primes.

	Bug,Bugh := UnitGroup(K);

	LS2basis := [Bugh(Bug.m) : m in [1 .. NumberOfGenerators(Bug)]];

	for t in S do
		a,b:=IsPrincipal(t);
		Append(~LS2basis,b);
	end for;
	
//	printf "LS2basis comes out as %o",[<t,Norm(t)> : t in LS2basis];
	return LS2basis;
end function;

function ProductWithSquarefreeNorm(RB,a)
	N := #RB;p := 2^N;
	K := Parent(RB[1]);

	CorrectlyNormed := [];

	for t in [0 .. p-1] do
		tt:=t;
		r:=K!1;
		for p in [1 .. N] do
			if (IsOdd(tt)) then
				r*:=RB[p];
			end if;
			tt div:= 2;
		end for;
		sfp:=SquareFree(Norm(r));
		if (sfp eq a) then
			Append(~CorrectlyNormed, r);
		end if;
	end for; 

	return CorrectlyNormed;
end function;

function IsFourByFourMatrix(M)
	return (Type(M) eq AlgMatElt) and (NumberOfRows(M) eq 4) and (NumberOfColumns(M) eq 4);
end function;

function QFFromMatrix(Coords,M)
	assert IsFourByFourMatrix(M);
	return &+[Coords[i]*Coords[j]*M[i,j] : i in [1..4], j in [1..4]];
end function;

function MatrixFromQF(Q)
	// Construct a 4x4 symmetric matrix from Q
	assert Type(Q) eq RngMPolElt;
	assert Rank(Parent(Q)) eq 4;
	assert TotalDegree(Q) eq 2;

	// We would like to *assert* that Parent(Q) is lexicographically
	// ordered ... not sure how to do that, though

	T:=Coefficients(Q);
	for i in [2,3,4,6,7,9] do
		T[i] /:= 2;
	end for;
	Reordering:=[1,2,5,3,6,8,4,7,9,10];
	T2:=[T[i] : i in Reordering];

	return SymmetricMatrix(T2);
end function;

function FourDescendentDiscriminant(A,B)
	assert IsFourByFourMatrix(A);
	assert IsFourByFourMatrix(B);
	return Discriminant(Determinant(A*x-B));
end function;

function FourDescendentSigmaInvariants(A,B)
// See MSS: these are the invariants [also the discriminant]

	assert IsFourByFourMatrix(A);
	assert IsFourByFourMatrix(B);
	DD := Determinant(A+x*B);
	return [Coefficient(DD,t) : t in [0..4]];
end function;

function FourDescendentCovariants(A,B)
// See MSS: these are the covariants

	assert IsFourByFourMatrix(A);
	assert IsFourByFourMatrix(B);

	S:=FourDescendentSigmaInvariants(A,B);

	adjs := Adjoint(x*Adjoint(A)+Adjoint(B));

	d1List := [Coefficient(adjs[i,j],2) / S[1] : i in [1..4], j in [1..4]];
	d2List := [Coefficient(adjs[i,j],1) / S[5] : i in [1..4], j in [1..4]];

	d1 := Matrix(Integers(),4,d1List);
	d2 := Matrix(Integers(),4,d2List);

	F1:=QFFromMatrix(d1);F2:=QFFromMatrix(d2);
	Q1:=QFFromMatrix(A);Q2:=QFFromMatrix(B);

	G:=Determinant(JacobianMatrix([F1,F2,Q1,Q2])) div 16;

	return [F1,F2,Q1,Q2,G];
end function;

function MerrimanQuarticForm(p)

	// p is a non-monic quartic polynomial
	// we want to return A and a monic quartic polynomial
	// such that Ay^2 = x^4 + ... is isomorphic to
	// y^2 = p. Note that A needs to be squarefree for the
	// later norm-comparison step

	// No need to do careful reduction, here [unless I want to for
	// fun] since we use the result only as a number field, and
	// we run OptimizedRepresentation() over it first.

	assert Type(p) eq RngUPolElt;
	assert Degree(p) eq 4;

	A := LeadingCoefficient(p);
	mc := A^3;
	zz := &+[(Coefficient(p,j)*mc*x^j) div A^j : j in [0..4]];
	A := SquareFree(A);

	return A,zz;
end function;

function MerrimanNiceNumberField(p)
	assert Type(p) eq RngUPolElt;
	assert Degree(p) eq 4;
	a,T := MerrimanQuarticForm(p);
	NastyNF<nnfv> := NumberField(T);
	NiceNF<nfv> := OptimizedRepresentation(NastyNF);
	return NastyNF, NiceNF,NiceNF!nnfv;
end function;

function MerrimanDescendent(epsilon)
	// The descendents corresponding to a particular epsilon
	// are the 4x4 matrices representing the quadratic forms
	// corresponding to the u3 and u4 coefficients in
	// epsilon*(u1 + u2 R + u3 R^2 + u4 R^3)

	K:=Parent(epsilon);
	Q<u1,u2,u3,u4>:=PolynomialRing(K,4);
	P:=PolynomialRing(Rationals());
	thugp := hom<K->P | P.1>;
	T := u1 + u2*K.1 + u3*K.1^2 + u4*K.1^3;
	ps := epsilon*T^2;
	R := Monomials(ps); S:=Coefficients(ps);
"thugp=",thugp;
// thugp(3*K.1^3 + K.1^7);
"K.1=",K.1;
"S[1]=",S[1];
"K!S[1]=",K!S[1];
"Type(S[1])=",Type(S[1]);
Parent(S[1]) eq K;
thugp(K!S[1]);  

       return 0; // that's all I care about for debugging...

end function;

function MerrimanDescent(p)
	a,T := MerrimanQuarticForm(p);
	Nasty, Nice, Theta := MerrimanNiceNumberField(p);
	DeltaG := Discriminant(p);

	// Construct the set of prime ideals in K which divide a or Delta(G);

	ONice := Integers(Nice);
	OA := Factorisation(ideal <ONice | a>);
	ODeltaG := Factorisation(ideal <ONice | DeltaG>);

	// Kill off the multiplicities and merge the sets, then compute L(S,2)

	SElts := {t[1] : t in OA} join {t[1] : t in ODeltaG};
	RB := LS2Basis(SElts);

	// Pull out the elements with the correct square-free norm
	
	RX := ProductWithSquarefreeNorm(RB,a);

	if (#RX eq 0) then
		print "This quartic has no descendents, by the norm condition, and so must represent an element of Sha";
	else
		printf "This quartic has %o descendents",#RX;
	end if;
	
	for t in RX do
		print MerrimanDescendent(Nasty!t);
	end for;
	return #RX;
end function;

r := 2*x^4 - 2*x^3 + 8*x - 12;

print MerrimanDescent(r);
