
/////////////////////////////////////////////////////////////////
//                                                             //
//           Explicit n-descent on elliptic curves             //
//                   Paper III. Algorithms                     // 
//   J.E. Cremona, T.A. Fisher, C. O'Neil, D.Simon, M.Stoll    //
//                                                             //
/////////////////////////////////////////////////////////////////

// This is a MAGMA file that checks the examples in Section 8.
// (It is NOT a general implementation of the algorithms in the paper.)
// 
// This version: T. Fisher, 31/03/2011, MAGMA V2.17-3
// 07/06/2011:  Some edits to match the paper.

printf "\nSection 8.1. An element of III[3]\n";

a4 := -1496259;
a6 := -693495810;
E := EllipticCurve([a4,a6]);
E0 := MinimalModel(E);
assert CremonaReference(E0) eq "681b1";
print "E =",CremonaReference(E0);
assert ThreeTorsionType(E) eq "Generic";

// A 3-torsion point on E
Q := Rationals();
P<X> := PolynomialRing(Q);
f := X^8 - 6*X^4 + 235*X^2 - 3;
L<u> := NumberField(f);
OL := RingOfIntegers(L);
xT := 12*u^6 -36*u^2 + 2115;
yT := -2820*u^7 - 144*u^5 + 16920*u^3 - 662268*u;
T := E(L)![xT,yT];
assert 3*T eq E!0;
lambdaT := (3*xT^2 + a4)/(2*yT); // slope of tangent line
assert lambdaT eq -3*u^7 + 15*u^3 - 705*u;

// An algebraic representative for a 3-Selmer element
alpha := 1/18*(u^6 - u^4 - 9*u^3 - 5*u^2 - 27*u - 3);

// Setting up fields L, L^+, M, M^+ and the maps between them
Lplus := sub<L|u^2>;
P<Y> := PolynomialRing(L);
g := ExactQuotient(f,Y^2-u^2);
assert g eq Y^6 + u^2*Y^4 + (u^4 - 6)*Y^2 + u^6 - 6*u^2 + 235;
M<v> := NumberField(g);
Mplus := sub<M|v^2>;
sigma := hom<L->L|-u>;
tau := hom<M->M|-v>;
u10 := 1/110322*(235*u^6 + 8*u^4 - 1175*u^2 + 24)*v^5 
  + 1/110322*(-235*u^7 - 8*u^5 + 1175*u^3 - 55185*u)*v^4 
  + 1/110322*(16*u^6 + 470*u^4 - 55241*u^2 + 1410)*v^3 
  + 1/55161*(-8*u^7 - 235*u^5 + 40*u^3 - 705*u)*v^2
  + 1/36774*(-235*u^6 - 18395*u^4 + 1175*u^2 + 18363)*v 
  + 1/36774*(235*u^7 + 8*u^5 - 1175*u^3 + 73572*u);
u01 := tau(u10);
assert {r[1]: r in Roots(f,M)} eq {u,-u,v,-v,u10,-u10,u01,-u01};
iota10 := hom<L->M|u10>;
iota01 := hom<L->M|u01>;
T10 := E(M)![iota10(x): x in Eltseq(T)];
T01 := E(M)![iota01(x): x in Eltseq(T)];
assert T01 + T10 eq T;

// The etale algebras R and (R tensor R)
R := car<Q,L>;
RR := car<Q,L,L,L,L,M>;

// Converting alpha to rho 
s := -u^2;
t := 1/486*(u^7 - u^5 - 5*u^3 - 27*u^2 + 240*u - 27)*v^4 
   + 1/486*(-2*u^7 + 2*u^5 - 27*u^4 + 10*u^3 - 237*u + 27)*v^2 
   + 1/162*(-u^7 - 9*u^6 + u^5 + 86*u^3 + 54*u^2 - 240*u + 45);
assert s^3 eq alpha*sigma(alpha);
assert t^3 eq sigma(alpha)*iota10(alpha)*iota01(alpha);
rho := RR!<1, 1, 1, sigma(alpha)/s,s,t/s>;
 
// The Weil pairing as an element of (R tensor R)
zeta3 := 1/110322*(-8*u^7 - 235*u^5 + 40*u^3 - 705*u)*v^5 
   + 1/110322*(-470*u^7 - 16*u^5 + 2350*u^3 - 55209*u)*v^3  
   + 1/36774*(8*u^7 + 235*u^5 + 18347*u^3 + 705*u)*v - 1/2;
assert zeta3^3 eq 1;
eps := RR!<1, 1, 1, 1, 1, zeta3 >;

// Computing a good basis for R (depending on alpha)
alphaconj := Conjugates(alpha);
aa := ideal<OL|alpha>;
ff := Factorization(aa);
FracOL := FieldOfFractions(OL);
ccinv := &*[PowerIdeal(FracOL)| f[1]^(-Floor(f[2]/3)): f in ff];
basL := [L!x : x in Basis(ccinv)];
for ctr in [1..1] do  // perhaps iterate a few times?
  MM := Matrix(8,8,[Conjugates(x): x in basL]);
  MMconj := Matrix(8,8,[ComplexConjugate(x): x in Eltseq(MM)]);
  D := DiagonalMatrix([Abs(x)^(2/3): x in alphaconj]);
  A := MM*D*Transpose(MMconj);
  A := Matrix(8,8,[Re(A[i,j]):i,j in [1..8]]);
  A := A + Transpose(A);
  _,cob := LLLGram(A);
  basL := [&+[cob[i,j]*basL[j]: j in [1..8]]: i in [1..8]];
end for;
basR := [R!<1,0>] cat [R|<0,x>: x in basL];
assert basR eq [ R | 
    <1, 0>,
    <0, 1>,
    <0, 1/3*(-u^7 + 6*u^3 - 235*u)>,
    <0, 1/18*(80*u^7 - 9*u^6 + u^5 - 481*u^3 + 54*u^2 + 18795*u - 2124)>,
    <0, 1/54*(9*u^7 - u^6 + u^4 - 54*u^3 + 5*u^2 + 2142*u - 240)>,
    <0, 1/54*(9*u^7 - u^6 + u^4 - 54*u^3 + 5*u^2 + 2088*u - 240)>,
    <0, 1/6*(-7*u^7 + u^6 + 42*u^3 - 7*u^2 - 1645*u + 237)>,
    <0, 1/54*(97*u^7 - 9*u^6 + 2*u^5 - 584*u^3 + 63*u^2 + 22785*u - 2133)>,
    <0, 1/54*(462*u^7 - 53*u^6 + 6*u^5 - u^4 - 2769*u^3 + 319*u^2 + 108549*u - 
        12423)>
];

// Structure constants for the obstruction algebra (checked below)
strconsts := [
    [ 1, 0, 0, 0, 0, 0, 0, 0, 0 ],
    [ 0, 1, 0, 0, 0, 0, 0, 0, 0 ],
    [ 0, 0, 1, 0, 0, 0, 0, 0, 0 ],
    [ 0, 0, 0, 1, 0, 0, 0, 0, 0 ],
    [ 0, 0, 0, 0, 1, 0, 0, 0, 0 ],
    [ 0, 0, 0, 0, 0, 1, 0, 0, 0 ],
    [ 0, 0, 0, 0, 0, 0, 1, 0, 0 ],
    [ 0, 0, 0, 0, 0, 0, 0, 1, 0 ],
    [ 0, 0, 0, 0, 0, 0, 0, 0, 1 ],
    [ 0, 1, 0, 0, 0, 0, 0, 0, 0 ],
    [ 0, 0, -3, 0, -3, -3, 0, 0, 0 ],
    [ 0, 2, 3, 0, 0, 0, 0, 3, 0 ],
    [ 12, 0, -3, 0, 3, 3, 0, 0, 0 ],
    [ 27, 2, -3, 0, 0, 0, -3, 0, 0 ],
    [ 27, -4, 0, 0, 0, 0, 3, -3, 0 ],
    [ 0, -2, -3, 3, 0, -6, 0, -3, -6 ],
    [ 0, 7, -3, 9, 0, 0, 0, -3, 0 ],
    [ -3, 4, 3, 0, -3, 6, -3, 3, 0 ],
    [ 0, 0, 1, 0, 0, 0, 0, 0, 0 ],
    [ 0, -4, 0, 0, 0, 0, 3, -3, 0 ],
    [ 8, 0, 1, -3, 3, 0, 0, 0, -3 ],
    [ -27, 0, 3, 2, 0, 0, 0, 3, 0 ],
    [ -4, -9, 2, 0, 2, 0, 0, 0, 0 ],
    [ -4, 9, -1, 6, -3, -4, 0, 0, -3 ],
    [ 28, 0, 1, 3, -3, 9, 2, 0, 3 ],
    [ 10, 0, 1, 3, 6, 9, 0, 2, 3 ],
    [ -50, -9, -2, -3, 3, 0, 0, -3, -1 ],
    [ 0, 0, 0, 1, 0, 0, 0, 0, 0 ],
    [ 12, 0, 6, 0, 6, 0, 0, 0, 0 ],
    [ -27, 0, -3, -4, 0, 0, -3, 0, 0 ],
    [ 0, 0, 0, 0, 0, 0, 0, 0, 0 ],
    [ 27, 0, 3, 2, 0, 0, 3, 0, 0 ],
    [ -54, 0, -6, 2, 0, 0, -6, 0, 0 ],
    [ 27, 0, 3, -5, 0, 0, 3, 0, 0 ],
    [ 27, 0, 3, -14, 0, 0, 3, 0, 0 ],
    [ 18, 0, 9, -2, 9, 0, 0, 0, 0 ],
    [ 0, 0, 0, 0, 1, 0, 0, 0, 0 ],
    [ 27, 2, 3, 0, 0, 0, 0, 3, 0 ],
    [ -4, 0, -1, -6, -1, 0, 0, 0, 3 ],
    [ 27, 0, -3, 2, 0, 0, 0, -3, 0 ],
    [ 8, 9, 0, 9, -2, 0, 0, 0, 0 ],
    [ -4, -9, -3, 3, -1, 2, 0, 0, 3 ],
    [ 4, 0, -6, 6, -8, -9, 2, 0, -3 ],
    [ 4, 0, 6, 6, 1, -9, 0, -4, -3 ],
    [ -23, 9, -3, -3, 1, 0, -9, 3, -1 ],
    [ 0, 0, 0, 0, 0, 1, 0, 0, 0 ],
    [ 27, 2, -3, 0, 0, 0, -3, 0, 0 ],
    [ -4, -9, 2, 0, 0, 2, 0, 0, 0 ],
    [ -54, 0, 6, -4, 0, 0, 0, 6, 0 ],
    [ -4, -9, 0, -6, 2, -4, 0, 0, -6 ],
    [ 8, 0, 0, -3, 0, 4, 0, 0, 6 ],
    [ 4, 9, 12, 0, 0, 7, -4, 0, 0 ],
    [ -50, 9, 9, 0, -9, -2, 0, 2, 0 ],
    [ 46, 0, 3, 0, 0, -2, 0, 3, -4 ],
    [ 0, 0, 0, 0, 0, 0, 1, 0, 0 ],
    [ 0, -5, 0, -9, 0, 0, -3, 3, 0 ],
    [ 28, 0, 1, 3, 6, 0, 2, 0, 3 ],
    [ 27, 0, -3, 16, 0, 0, 0, -3, 0 ],
    [ 4, 9, 12, 0, 7, 0, -4, 0, 0 ],
    [ 4, -9, -6, -6, -6, -5, 2, 0, 3 ],
    [ 98, 0, 6, -3, -6, -9, 5, 0, -3 ],
    [ 35, 0, -3, 15, -15, -9, 7, -11, -3 ],
    [ 50, 9, 9, 3, 6, 0, -2, 3, -8 ],
    [ 0, 0, 0, 0, 0, 0, 0, 1, 0 ],
    [ 0, -5, 6, -6, 6, 0, -3, 3, -6 ],
    [ 10, 0, -5, 3, -3, -9, 0, -4, 3 ],
    [ 27, 0, -3, 7, 0, 0, 0, -3, 0 ],
    [ 4, 9, 9, 0, 7, 9, 0, 2, 0 ],
    [ -50, -9, -15, -6, 3, -5, 0, 2, 3 ],
    [ 35, 0, 3, -3, 3, 0, 7, -5, -3 ],
    [ 134, 0, -6, -3, -6, 0, 0, -7, -3 ],
    [ 68, -18, 9, -6, 6, 0, 0, 1, 1 ],
    [ 0, 0, 0, 0, 0, 0, 0, 0, 1 ],
    [ -3, -2, -3, 0, 3, 0, 3, 0, 0 ],
    [ -50, 0, 1, -3, 0, 0, 3, 0, 2 ],
    [ 18, 0, 3, 4, -9, -9, 0, -6, 0 ],
    [ -23, 0, -3, 9, -2, 0, 6, -9, 2 ],
    [ 46, 0, -3, 3, 0, -2, -3, 0, -4 ],
    [ 50, -27, -15, -6, 0, 0, 1, 0, 7 ],
    [ 68, 0, -3, 3, 0, 0, -3, -2, 7 ],
    [ -82, 0, 6, 0, 0, 9, 0, -3, 2 ]
];

function AlgebraToVectorSpace(basR)
  V := VectorSpace(Q,9);
  bb := [&cat[Eltseq(x): x in b]: b in basR];
  psi := hom<V->V|Matrix(9,9,bb)>;
  mapRtoV := func<r|V!(&cat[Eltseq(x): x in r]) @@ psi>;
  return V,mapRtoV;
end function;

function TraceMap(rr)
  x1 := rr[1] + Trace(rr[5]);
  x2 := rr[2] + rr[3] + rr[4] + Trace(rr[6]);
  return R!<x1,x2>;
end function;

function TensorMap(r1,r2)
  a1,b1 := Explode(r1);
  a2,b2 := Explode(r2);
  bb1 := sigma(b1);
  bb2 := sigma(b2);
  c := iota10(b1)*iota01(b2);
  return elt<RR|a1*a2,b1*a2,a1*b2,bb1*bb2,b1*bb2,c>;
end function;

// Computing the structure constants
BB := [TensorMap(basR[i],basR[j]):i,j in [1..9]];
V,mapRtoV := AlgebraToVectorSpace(basR);
cc := [ mapRtoV(TraceMap(<eps[i]*rho[i]*bb[i] : i in [1..6]>)): bb in BB];
// Creating the obstruction algebra
A := AssociativeAlgebra< V | cc : Check := true >;
// Checking structure contstants are as above.
assert strconsts eq [Eltseq(A.i*A.j): i,j in [1..9]];

assert Identity(A) eq A.1;
assert A.2^2 eq -3*A.3 - 3*A.5 - 3*A.6;
assert A.2*A.3 eq 2*A.2 + 3*A.3 + 3*A.8; 
assert A.2*A.4 eq 12*A.1 - 3*A.3 + 3*A.5 + 3*A.6; 
assert A.2*A.5 eq 27*A.1 + 2*A.2 - 3*A.3 - 3*A.7; 
assert A.2*A.6 eq 27*A.1 - 4*A.2 + 3*A.7 - 3*A.8; 
assert A.2*A.7 eq -2*A.2 - 3*A.3 + 3*A.4 - 6*A.6 - 3*A.8 - 6*A.9; 
assert A.2*A.8 eq 7*A.2 - 3*A.3 + 9*A.4 - 3*A.8; 
assert A.2*A.9 eq -3*A.1 + 4*A.2 + 3*A.3 - 3*A.5 + 6*A.6 - 3*A.7 + 3*A.8;

// Printing out some structure constants
// S<a1,a2,a3,a4,a5,a6,a7,a8,a9> := PolynomialRing(Q,9);
// for i in [2..2] do
// for j in [2..9] do
//   printf "a%o*a%o = %o \\\\ \n",i,j,&+[Eltseq(A.i*A.j)[k]*S.k: k in [1..9]];
// end for; end for;
// print [MinimalPolynomial(x): x in Basis(A)];

assert [MinimalPolynomial(x): x in Basis(A)] eq [ X - 1, X^3 + 162, 
  X^3 - 12*X - 227, X^2, X^3 - 12*X - 470, X^3 - 12*X - 470, 
  X^3 - 147*X - 367, X^3 - 201*X + 1307, X^3 + 123*X + 254 ];

// Computing the discriminant
function MyDiscriminant(A:bas := Basis(A))
  _, mm := MatrixAlgebra(A);
  trd := func<x|(1/3)*Trace(mm(x))>;
  trmat := Matrix(9,9,[trd(bas[i]*bas[j]): i,j in [1..9]]);
//  for i in [1..9] do
//    for j in [1..9] do
//      printf "%o & ",trmat[i,j]; 
//    end for;
//    print " &";
//  end for;
//  print trmat;
  return Determinant(trmat);
end function;

assert Abs(Discriminant(OL)) eq 3^11*227^4;
assert Abs(MyDiscriminant(A)) eq 3^20*227^4;
assert IsPower(alpha*OL,3);

// A basis for a maximal order (checked below)
mat0 := Matrix(9,9,[ 
    [ 1/3, 0, 0, 0, 0, 56/3, 42, 101/3, 2438/3 ],
    [ 0, 1/681, 0, 38/2043, 157/681, 95432/2043, 5561/227, 
           25292/681, 1408079/2043 ],
    [ 0, 0, 1/2043, 167/2043, 181/681, 175106/2043, 254/9, 
           29086/681, 8248/9 ],
    [ 0, 0, 0, 1/9, 0, 529/9, 0, 0, 2041/9 ],
    [ 0, 0, 0, 0, 1/3, 53, 35, 160/3, 934 ],
    [ 0, 0, 0, 0, 0, 1/3, 0, 0, 1/3 ],
    [ 0, 0, 0, 0, 0, 0, 1/3, 0, 8/3 ],
    [ 0, 0, 0, 0, 0, 0, 0, 1/3, 8/3 ],
    [ 0, 0, 0, 0, 0, 0, 0, 0, 1 ]
]); 

bb0 := [&+[mat0[i,j]*A.j: j in [1..9]]: i in [1..9]];
assert bb0 eq [
  (1/3)*(A.1 + 56*A.6 + 126*A.7 + 101*A.8 + 2438*A.9),
  (1/2043)*(3*A.2 + 38*A.4 + 471*A.5 + 95432*A.6 + 50049*A.7 
      + 75876*A.8 + 1408079*A.9),
  (1/2043)*(A.3 + 167*A.4 + 543*A.5 + 175106*A.6 + 57658*A.7 
      + 87258*A.8 + 1872296*A.9),
  (1/9)*(A.4 + 529*A.6 + 2041*A.9),
  (1/3)*(A.5 + 159*A.6 + 105*A.7 + 160*A.8 + 2802*A.9),
  (1/3)*(A.6 + A.9),(1/3)*(A.7 + 8*A.9),(1/3)*(A.8 + 8*A.9),A.9];

// for i in [1..9] do
//   dd := RationalGCD(Eltseq(mat1[i]));
//   printf "b%o = %o(%o) \\\\ \n",i,dd,(1/dd)*&+[mat1[i,j]*S.j: j in [1..9]];
// end for;

// Computing maximal order (and repeating until the answer agrees 
// with the one I got earlier)
while true do 
  bb := Basis(MaximalOrder(A));  
  if bb eq bb0 then break; end if;
  delete A`MaximalOrder; // stop answer being cached
end while;

assert Abs(MyDiscriminant(A:bas:=bb)) eq 1;
assert [MinimalPolynomial(b): b in bb] eq [ 
   X^3 - X^2 + 67882988*X + 153570178243, 
   X^3 + 46000395*X + 93752525874, 
   X^3 + 80434914*X + 198363227932, 
   X^3 + 4444433*X + 1099577331, X^3 + 84844655*X + 243745052250, 
   X^3 - 3*X, X^3 + 725*X + 3507, X^3 + 671*X + 9393, 
   X^3 + 123*X + 254 ];

// Trivialising the obstruction algebra
mat := LLL(Matrix(Q,9,9,[Eltseq(b): b in bb]));
bb1 := [A!mat[i]: i in [1..9]];

assert bb1 eq [
  (1/2043)*(-12*A.2 - 63*A.3 - 4*A.4 - 73*A.6 + 18*A.7 + 8*A.9),
  (1/2043)*(-81*A.2 - 28*A.3 - 27*A.4 + 18*A.6 + 8*A.7 + 54*A.9),
  (1/681)*(-15*A.3 - 8*A.4 + 27*A.5 + 27*A.6 + 12*A.8),
  (1/2043)*(-27*A.2 + 28*A.3 + 21*A.4 + 69*A.5 + 75*A.6 
              + 154*A.7 - 45*A.8 + 18*A.9), 
  (1/2043)*(-15*A.2 - 91*A.3 + 49*A.4 - 12*A.5 + 67*A.6 - 91*A.7 
              - 81*A.8 + 10*A.9), 
  (1/2043)*(138*A.2 - 55*A.3 + 54*A.4 - 27*A.5 + 18*A.6 + 20*A.7 
              - 12*A.8 + 135*A.9), 
  (1/2043)*(-27*A.2 + 83*A.3 + 126*A.4 - 30*A.5 - 24*A.6 - 73*A.7 
              + 138*A.8 + 18*A.9),
  (1/2043)*(-36*A.2 + 37*A.3 + 48*A.4 + 138*A.5 - 81*A.6 - 173*A.7 
	     - 90*A.8 + 24*A.9),
  (1/2043)*(-681*A.1 - 27*A.2 - 85*A.3 - 9*A.4 + 6*A.6 - 73*A.7 + 18*A.9)
];

// for i in [1..9] do
//  dd := RationalGCD(Eltseq(mat[i]));
//    printf "b%o = %o(%o) \\\\ \n",i,dd,(1/dd)*&+[mat[i,j]*S.j: j in [1..9]];
// end for;

assert Abs(MyDiscriminant(A:bas:=bb1)) eq 1;
assert [MinimalPolynomial(b): b in bb1] eq
  [ X^2, X^2, X^2, X^3 - X, X^3 - X, X^3, X^3 - X, X^3 - X, X^2 + X ];

trivdata := [
    [ 1, 0, 0, 0, 1, 0, 0, 0, 1 ],
    [ 6, -6, 3, 6, -6, 0, 0, -9, 0 ],
    [ -4, -3, -3, 0, 2, -3, 9, -9, 2 ],
    [ 0, 0, 0, 0, 0, 0, 12, 15, 0 ],
    [ 2, 3, 3, 0, -4, 3, 18, -18, 2 ],
    [ 2, 9, -6, 0, 2, -6, -9, 9, -4 ],
    [ -5, 3, 3, 0, -11, 3, -9, -9, 16 ],
    [ 1, -12, 3, -12, -8, 3, -9, 9, 7 ],
    [ 7, -9, -3, 9, -11, 6, 15, -15, 4 ]
];
MM := [Matrix(Q,3,3,x) : x in trivdata];

procedure printmat(M)
  print "\\begin{pmatrix}";
  for i in [1..3] do
    printf "%o & %o & %o ",M[i,1],M[i,2],M[i,3];
    if i lt 3 then print "\\\\"; else print ""; end if;
  end for;
  print "\\end{pmatrix} &";
end procedure;

// for i in [1..9] do 
//  printf "\\av_%o \\mapsto \n",i; 
//  printmat(MM[i]); 
// end for;

// Checking the trivialisation
assert forall{[i,j] : i,j in [1..9] |
  MM[i]*MM[j] eq &+[(A.i*A.j)[k]*MM[k]: k in [1..9]]};

// Computing the 27 quadrics
PP<[X]> := PolynomialRing(Q,8);
MC := MonomialCoefficient;
MD := MonomialsOfDegree;
mons := &join[MD(PP,d) : d in [0..2]];
// The type I quadrics
  LX := PolynomialRing(L,8);
  zT := &+[LX| basL[i]*LX.i : i in [1..8]];
  zmT := &+[LX| sigma(basL[i])*LX.i : i in [1..8]];
  Q1 := xT + rho[5]*zT*zmT;
  Q1 := PolynomialRing(Lplus,8)!Q1;
  coeffs := [Eltseq(MC(Q1,Parent(Q1)!mon)): mon in mons];
  quads1 := [&+[coeffs[i][j]*mons[i]: i in [1..#mons]]: j in [2..4]];
// The type II quadrics
  MX := PolynomialRing(M,8);
  zT1 := &+[MX| iota10(basL[i])*MX.i : i in [1..8]];
  zT2 := &+[MX| iota01(basL[i])*MX.i : i in [1..8]];
  zT := &+[MX| basL[i]*MX.i : i in [1..8]];
  zmT := &+[MX| sigma(basL[i])*MX.i : i in [1..8]];
  c := (1/3)*(iota10(lambdaT) + iota01(lambdaT) + 2*lambdaT);
  Q2 := c*zT - rho[4]*zmT^2 + rho[6]*zT1*zT2;
  Q2 := PolynomialRing(Mplus,8)!Q2;
  AbsEltseq := func<x|&cat[Eltseq(y): y in Eltseq(x)]>;
  coeffs := [AbsEltseq(MC(Q2,Parent(Q2)!mon)) : mon in mons];
  quads2 := [&+[coeffs[i][j]*mons[i]: i in [1..#mons]]: j in [1..24]];

function TidyUpPolynomials(polys,mons)
  mat := Matrix(#polys,#mons,[MC(f,mon): mon in mons,f in polys]);
  mat := LLL(Matrix(Basis(PureLattice(Lattice(mat)))));
  return [&+[mat[i,j]*mons[j]: j in [1..#mons]]: i in [1..#polys]];
end function;

quads := TidyUpPolynomials(quads1 cat quads2,mons);

// Projecting away from the "T = 0" component
mons1 := [1] cat [PP.i: i in [1..8]];
mat := Matrix(Q,27,9,[MC(q,mon): mon in mons1,q in quads]);
km := KernelMatrix(mat);
quads := [&+[km[i,j]*quads[j]: j in [1..27]]: i in [1..18]];
quads := TidyUpPolynomials(quads,mons);

printed_quads := [
    X[4]*X[6] - X[5]*X[6] - X[5]*X[7] + X[6]*X[8] + X[7]*X[8],
    2*X[1]*X[6] - X[3]*X[6] + X[4]*X[5] + X[5]^2 - X[5]*X[8] - 2*X[6]*X[7]
        + X[6]*X[8],
    -X[1]*X[6] - X[3]*X[6] - 2*X[4]*X[5] + X[4]*X[8] + X[5]^2 - X[5]*X[8] 
        + X[6]*X[8],
    -X[1]*X[6] + X[2]*X[4] - X[3]*X[6] + X[4]*X[5] - X[4]*X[7] - X[5]^2 - 
        X[5]*X[6] - X[5]*X[7] + X[5]*X[8] - X[6]*X[8] + X[7]*X[8],
    -X[2]*X[4] + X[2]*X[5] - X[4]*X[5] - X[4]*X[6] + X[4]*X[7] + X[5]*X[6]
        + X[5]*X[7] + X[6]*X[8] + X[7]*X[8],
    -X[1]*X[8] - X[2]*X[6] + X[2]*X[7] - X[3]*X[8] + X[5]*X[6] - X[5]*X[7]
        + X[6]^2 - X[6]*X[8] - X[7]^2 + X[7]*X[8] + X[8]^2,
    -2*X[1]*X[7] + X[3]*X[7] - X[4]^2 - X[4]*X[8] + X[5]^2 + X[6]^2 + 
        X[6]*X[7] - X[7]*X[8],
    -X[1]*X[7] + X[2]*X[5] + 2*X[3]*X[7] + X[4]^2 - X[4]*X[5] - X[4]*X[6] 
        - X[6]*X[8] + 2*X[7]^2,
    -X[1]*X[5] + 2*X[1]*X[6] + 2*X[1]*X[7] + X[2]*X[6] + X[2]*X[7] - 
        X[3]*X[5] - X[3]*X[6] - X[3]*X[7] + X[4]^2 + X[4]*X[5] - X[4]*X[6]
        + X[4]*X[8] + X[6]^2 + X[6]*X[7] + X[6]*X[8] - X[7]^2 + X[7]*X[8],
    -X[1]*X[4] + 2*X[1]*X[7] - X[2]*X[7] + 2*X[3]*X[4] - X[3]*X[7] + 
        X[4]^2 + X[4]*X[7] - X[4]*X[8] - X[5]^2 + 2*X[6]^2 - 2*X[7]^2 + 
        X[7]*X[8],
    2*X[1]*X[4] - X[1]*X[5] - X[2]*X[6] - X[3]*X[4] - X[3]*X[5] + 
        X[4]*X[6] - X[4]*X[7] + X[4]*X[8] + X[5]*X[8] - 2*X[6]^2 - 
        2*X[6]*X[7] - 3*X[7]^2,
    -X[1]*X[5] - 2*X[1]*X[6] + X[2]*X[6] + 2*X[3]*X[5] + X[3]*X[6] - 
        X[4]*X[5] + X[4]*X[6] - X[5]^2 + 2*X[5]*X[7] - X[5]*X[8] + 
        2*X[6]^2 - 2*X[6]*X[7] - X[6]*X[8],
    -X[1]*X[2] - 2*X[1]*X[5] + X[1]*X[6] + X[1]*X[7] - X[2]*X[3] + 
        X[2]*X[8] - 2*X[3]*X[5] + X[3]*X[6] + X[3]*X[7] - 3*X[4]*X[8] - 
        X[5]*X[8] - X[6]*X[8] - X[7]*X[8],
    -2*X[1]*X[2] - X[1]*X[6] - 2*X[1]*X[7] + X[2]*X[3] - X[2]*X[8] - 
        3*X[3]*X[5] + 2*X[3]*X[6] + X[3]*X[7] + 3*X[3]*X[8] - 2*X[4]^2 + 
        X[4]*X[6] - X[5]^2 - X[5]*X[6] - X[5]*X[7] - X[6]^2 - X[6]*X[8],
    2*X[1]^2 + X[1]*X[3] - X[1]*X[6] - X[1]*X[8] + 3*X[2]*X[4] - X[3]^2 - 
        X[3]*X[6] + 2*X[3]*X[8] - X[4]*X[5] + 3*X[4]*X[6] + 3*X[4]*X[7] + 
        2*X[4]*X[8] - X[5]^2 + 3*X[5]*X[7] + X[5]*X[8] - 2*X[6]*X[8] - 
        X[8]^2,
    4*X[1]^2 - 4*X[1]*X[3] - X[1]*X[5] + X[1]*X[6] - 2*X[1]*X[7] + 
        4*X[1]*X[8] + 2*X[2]*X[4] + X[2]*X[5] + X[2]*X[6] + X[2]*X[7] - 
        X[2]*X[8] + X[3]^2 - X[3]*X[5] - 2*X[3]*X[6] + X[3]*X[7] - 
        2*X[3]*X[8] + X[4]*X[5] + 4*X[4]*X[6] + X[4]*X[7] + 2*X[4]*X[8] - 
        X[5]^2 + 4*X[5]*X[6] + X[5]*X[7] + X[5]*X[8] + X[6]^2 + X[6]*X[7] 
        + X[6]*X[8] - X[7]^2 + X[7]*X[8] + X[8]^2,
    2*X[1]^2 + X[1]*X[3] + 2*X[1]*X[4] - 2*X[1]*X[5] + 5*X[1]*X[8] - 
        X[2]^2 + X[2]*X[4] + X[2]*X[5] - X[2]*X[6] - 4*X[2]*X[7] - 
        X[2]*X[8] - X[3]^2 - X[3]*X[4] + X[3]*X[5] - X[3]*X[8] - X[4]*X[5]
        + 4*X[4]*X[7] + 2*X[4]*X[8] - X[5]*X[7] - X[5]*X[8] + 2*X[6]^2 + 
        X[6]*X[7] + 2*X[7]^2 + 2*X[8]^2,
    X[1]^2 - X[1]*X[3] - X[1]*X[4] + X[1]*X[5] + 2*X[1]*X[6] - X[1]*X[7] -
        X[2]*X[4] + 2*X[2]*X[5] - 2*X[2]*X[6] - 2*X[3]^2 + 2*X[3]*X[4] - 
        2*X[3]*X[5] - X[3]*X[6] - X[3]*X[7] + 3*X[3]*X[8] + 2*X[4]^2 - 
        X[4]*X[5] + X[4]*X[6] - 4*X[4]*X[7] - X[4]*X[8] + X[5]*X[6] + 
        3*X[5]*X[7] + X[5]*X[8] + X[6]^2 + 2*X[6]*X[7] + X[6]*X[8] - 
        X[7]^2 - X[8]^2
];

assert ideal<PP|quads> eq ideal<PP|printed_quads>;

// Making changes of co-ordinates corresponding to the 
// trivialisation and the fudge factor 1/yT
T1 := Matrix(Q,9,9,[mapRtoV(<0,r[2]*yT>): r in basR]);
T2 := Matrix(Q,9,9,[Eltseq(mat):mat in MM]);
mat := Transpose(T2^(-1)*T1);

assert (1/4)*mat eq Matrix(Rationals(),9,9,[ 
  0, 0, 0, 0, 0, 0, 0, 0, 0, 
  -70, -54, 144, 12, 86, 423, 6, -48, -16, 
  360, 615, 1128, -444, -510, 663, -372, 207, 150, 
  160, -180, -621, 222, -257, 648, 174, 291, 97, 
  -75, -48, 195, -195, 21, -285, -69, -81, 54, 
  3, -24, -267, 24, -192, -21, 87, 81, 189, 
  -69, -81, -108, -9, 57, -135, 117, 36, 12, 
  -105, -81, -270, 18, 129, 27, 9, -72, -24, 
  252, -72, -801, 72, -333, 423, 261, 243, 81 ]);

SS<z11,z12,z13,z21,z22,z23,z31,z32,z33> := PolynomialRing(Q,9);
subst := [&+[mat[i,j]*SS.j: j in [1..9]]: i in [2..9]];
quads := [Evaluate(q,subst) : q in quads];
quads := TidyUpPolynomials(quads,MD(SS,2));
assert #quads eq 18;

// Substituting z_ij = x_i y_j
RR<x1,x2,x3,y1,y2,y3> := PolynomialRing(Q,6);
zz := [[x1,x2,x3][i]*[y1,y2,y3][j]: i,j in [1..3]];
quads := [Evaluate(f,zz): f in quads];

// Extracting a cubic
R<x,y,z> := PolynomialRing(Q,3);
monsx := [Evaluate(mon,[x1,x2,x3]): mon in MD(R,3)];
monsy := [Evaluate(mon,[y1,y2,y3]): mon in MD(R,2)];
mons := [mx*my : mx in monsx, my in monsy];
forms := [x*q : x in [x1,x2,x3], q in quads];
mat := Matrix(54,50,[MC(f,mons[i]): i in [11..60], f in forms]);
km := KernelMatrix(mat);
cubic := &+[km[1,i]*forms[i]: i in [1..54]];
cubic := (-1/1548)*Evaluate(cubic,[x,y,z,1,0,0]);
model := GenusOneModel(cubic);
assert IsIsomorphic(Jacobian(model),E);
assert Discriminant(model) eq Discriminant(E0); // already minimal
assert forall{p : p in BadPrimes(E) | IsLocallySoluble(model,p)};
assert cubic eq 3*x^3 - 13*x^2*y + 4*x^2*z + 2*x*y^2 + x*y*z 
                  - y^3 - 5*y^2*z - y*z^2 + z^3;
assert Reduce(cubic) eq x^3 + 2*x^2*z + 4*x*y^2 + x*y*z - 4*x*z^2 
                    + 3*y^3 - 9*y^2*z + 3*y*z^2 - 6*z^3;

FF := [
  3*x^3 - 13*x^2*y + 4*x^2*z + 2*x*y^2 + x*y*z 
	- y^3 - 5*y^2*z - y*z^2 + z^3,
  x^3 + 6*x^2*y + 4*x^2*z + 4*x*y^2 + 5*x*y*z + 2*x*z^2 + y^3 
	- 3*y^2*z + 7*y*z^2 + 6*z^3,
  x^3 - 2*x^2*y - x^2*z - 7*x*y*z + 8*x*z^2 + 4*y^3 
        - 5*y^2*z + 6*y*z^2 + z^3,
  x^3 - 2*x^2*z + 4*x*y^2 + 3*x*y*z - 5*x*z^2 - y^3 + 6*y^2*z 
        + 2*y*z^2 + 7*z^3 
];

assert cubic eq FF[1];
print "F1 =",cubic;

// print "Computing the other elements in III[3]";
// // SetVerbose("ThreeDescent",2);
// // SetClassGroupBounds(500);
// time three_desc := ThreeDescent(E0);
// print "Checking answer";
// cubics0 := [GenusOneModel(x): x in FF];
// cubics1 := [GenusOneModel(x): x in three_desc];
// assert #cubics1 eq 4;
// time for m0 in cubics0 do 
//   assert exists{m1 : m1 in cubics1 | IsEquivalent(m0,m1:E := E0)};
// end for;

//////////////////////////////////////////////////////////////

printf "\nSection 8.2. Adding ternary cubics\n";

Q := Rationals();
R<x,y,z> := PolynomialRing(Q,3); 

E := EllipticCurve("126a3");
cubics := [
    x^3 + x*y^2 - x*y*z + x*z^2 + y^3 + 3*y^2*z - 6*y*z^2 + z^3,
    2*x^2*y + 2*x^2*z + 3*x*y^2 - x*y*z + 2*x*z^2 + 2*y*z^2 + 2*z^3
];
models := [GenusOneModel(c): c in cubics];
assert cInvariants(models[1]) eq cInvariants(E);
assert cInvariants(models[2]) eq cInvariants(E);
assert Discriminant(E) eq -2^6*3^6*7^3;
assert [IsLocallySoluble(models[1],p): p in [2,3,7]] eq [false,true,false];

K<zeta3> := CyclotomicField(3);
RK := PolynomialRing(K,3);
sigma := hom<K->K|zeta3^2>;
P := E![1,13];
Q := E(K)![-6,21*zeta3 + 13];
assert 3*P eq E!0;
assert 3*Q eq E!0;
assert E(K)![sigma(x): x in Eltseq(P)] eq P;
assert E(K)![sigma(x): x in Eltseq(Q)] eq -Q;

TT := ThreeTorsionPoints(E);
assert TT eq <P,-P,-Q,-Q+P,-Q-P>;

MM := ThreeTorsionMatrices(E,models[1]);
MP := (-1/336)*MM[1];
MQ := (1/252*(zeta3 - 1))*Matrix(K,3,3,[sigma(MM[3][i,j]):i,j in [1..3]]);

assert MP eq Matrix(3,3,[ 1, 0, 0, 0, -1, 1, 0, -1, 0 ]);
assert MQ eq Matrix(3,3,[ 0, 1 - zeta3, 1 + 2*zeta3, 
  1, 1, -1 - zeta3, 1 + zeta3, -zeta3, -1 ]);
assert cubics[1]^MP eq Determinant(MP)*cubics[1];
assert (RK!cubics[1])^MQ eq Determinant(MQ)*(RK!cubics[1]);
assert Determinant(MP) eq 1;
assert Determinant(MQ) eq -zeta3^2*(1-zeta3)^3;

NN := ThreeTorsionMatrices(E,models[2]);
NP := (1/48)*NN[1];
NQ := 1/72*(zeta3 + 2)*Matrix(K,3,3,[sigma(NN[3][i,j]):i,j in [1..3]]);

assert NP eq Matrix(3,3,[ 0, -3, -1, 2, 2, -2, -2, 0, -2 ]);
assert NQ eq Matrix(3,3,[ 1, 4 + 2*zeta3, -4*zeta3, 
  4 + 2*zeta3, 1 - zeta3, 2 + 4*zeta3, -2*zeta3, 1 + 2*zeta3, -2 + zeta3]);
assert cubics[2]^NP eq Determinant(NP)*cubics[2];
assert (RK!cubics[2])^NQ eq Determinant(NQ)*(RK!cubics[2]);
assert Determinant(NP) eq -28;
assert Determinant(NQ) eq (4-zeta3)^3*(1-2*zeta3^2)/(1-2*zeta3);

///////////////////////////////////////////////////////////////

// Second example

Q := Rationals();
R<x,y,z> := PolynomialRing(Q,3); 

E := EllipticCurve("1722f1");
assert ThreeTorsionType(E) eq "Generic";
print "E =",CremonaReference(E);

// A 3-torsion point on E
P<X> := PolynomialRing(Q);
f := X^8 + 234*X^4 + 1256*X^2 - 4563;
L<u> := NumberField(f);
OL := RingOfIntegers(L);
T := E(L)![1/192*(u^6 + 9*u^4 + 315*u^2 + 1979), 
   1/44928*(-643*u^7 - 117*u^6 - 1755*u^5 - 1053*u^4 
	    - 166257*u^3 - 36855*u^2 - 888689*u - 254007) ];
assert 3*T eq E!0;
E`ThreeTorsionPoints := <T>;

cubics := [
    x^3 - 2*x^2*z + 2*x*y^2 + x*y*z + 3*x*z^2 + y^3 + 3*y^2*z - y*z^2 + 2*z^3,
    3*x^2*y + x^2*z - x*y^2 + 3*x*y*z - 2*x*z^2 + y^3 + 6*y*z^2 + z^3
];
printf "F1 = %o\nF2 = %o\n",cubics[1],cubics[2];
models := [GenusOneModel(c): c in cubics];
assert cInvariants(models[1]) eq cInvariants(E);
assert cInvariants(models[2]) eq cInvariants(E);
assert Discriminant(E) eq -2^8*3^3*7^3*41; 
assert [IsLocallySoluble(models[1],p): p in [2,3,7,41]] 
  eq [true,false,false,true];

alphas := [ThreeSelmerElement(E,model)[1]: model in models];
alpha1 := 1/13312*(-11*u^7 - 65*u^6 - 39*u^5 - 117*u^4 
           - 2561*u^3 - 16419*u^2 - 20173*u - 126503);
alpha2 := 1/6656*(-253*u^7 + 364*u^6 - 793*u^5 + 1092*u^4 
            - 58695*u^3 + 81172*u^2 - 457635*u + 616252);
assert IsPower(alphas[1]/alpha1,3) or IsPower(alphas[1]*alpha1,3);
assert IsPower(alphas[2]/alpha2,3) or IsPower(alphas[2]*alpha2,3);

alpha := alpha1*alpha2;

// Setting up fields L and M and the maps between them
P<Y> := PolynomialRing(L);
g := Y^6 + u^2*Y^4 + (u^4 + 234)*Y^2 + u^6 + 234*u^2 + 1256;
M<v> := NumberField(g);
sigma := hom<L->L|-u>;
tau := hom<M->M|-v>;
u10 := 1/52396032*(-157*u^6 - 1521*u^4 - 30615*u^2 + 177957)*v^5 + 
    1/52396032*(157*u^7 + 1521*u^5 + 30615*u^3 + 493787*u)*v^4 + 
    1/26198016*(-1521*u^6 + 6123*u^4 + 39277*u^2 - 716391)*v^3 + 
    1/671744*(39*u^7 - 157*u^5 + 7605*u^3 + 18369*u)*v^2 + 
    1/52396032*(-18369*u^6 + 493787*u^4 - 3581955*u^2 + 47018985)*v + 
    1/1343488*(471*u^7 + 4563*u^5 + 91845*u^3 + 2153105*u);
u01 := tau(u10);
assert {r[1]: r in Roots(f,M)} eq {u,-u,v,-v,u10,-u10,u01,-u01};
iota10 := hom<L->M|u10>;
iota01 := hom<L->M|u01>;
T10 := E(M)![iota10(x): x in Eltseq(T)];
T01 := E(M)![iota01(x): x in Eltseq(T)];
assert T01 + T10 eq T;

// The etale algebras R and (R tensor R)
R := car<Q,L>;
RR := car<Q,L,L,L,L,M>;

// Converting alpha to rho 
s := 1/128*(-5*u^6 - 17*u^4 - 1183*u^2 - 8779);
t := 1/26198016*(-2513*u^7 - 23205*u^5 - 490035*u^3 + 358176*u^2 
  - 4099543*u + 3939936)*v^4 + 1/13099008*(-245*u^7 + 98007*u^5 
  + 179088*u^4 + 1743105*u^3 - 1612387*u - 33131280)*v^2 
  + 1/26198016*(-799141*u^7 + 358176*u^6 - 2714985*u^5 
  - 196664559*u^3 + 83813184*u^2 - 1533143171*u + 596363040);

assert s^3 eq alpha*sigma(alpha);
assert t^3 eq sigma(alpha)*iota10(alpha)*iota01(alpha);
rho := RR!<1, 1, 1, sigma(alpha)/s,s,t/s>;

// The Weil pairing as an element of (R tensor R)
zeta3 := 1/52396032*(39*u^7 - 157*u^5 + 7605*u^3 + 18369*u)*v^5 + 
    1/26198016*(-157*u^7 - 1521*u^5 - 30615*u^3 - 157915*u)*v^3 + 
     1/52396032*(4563*u^7 - 18369*u^5 + 1561529*u^3 + 2149173*u)*v - 1/2;
assert zeta3^3 eq 1;
eps := RR!<1, 1, 1, 1, 1, zeta3 >;

// Computing a good basis for R (depending on alpha)
alphaconj := Conjugates(alpha);
aa := ideal<OL|alpha>;
ff := Factorization(aa);
FracOL := FieldOfFractions(OL);
ccinv := &*[PowerIdeal(FracOL)| f[1]^(-Floor(f[2]/3)): f in ff];
basL := [L!x : x in Basis(ccinv)];
for ctr in [1..1] do  // perhaps iterate a few times?
  MM := Matrix(8,8,[Conjugates(x): x in basL]);
  MMconj := Matrix(8,8,[ComplexConjugate(x): x in Eltseq(MM)]);
  D := DiagonalMatrix([Abs(x)^(2/3): x in alphaconj]);
  A := MM*D*Transpose(MMconj);
  A := Matrix(8,8,[Re(A[i,j]):i,j in [1..8]]);
  A := A + Transpose(A);
  _,cob := LLLGram(A);
  basL := [&+[cob[i,j]*basL[j]: j in [1..8]]: i in [1..8]];
end for;
basR := [R!<1,0>] cat [R|<0,x>: x in basL];
assert basR eq [ R | 
    <1, 0>,
    <0, 1/39936*(15*u^7 + 13*u^6 - 65*u^5 + 13*u^4 + 4069*u^3 + 
        3055*u^2 - 2675*u - 10569)>,
    <0, 1/19968*(-7*u^7 + 26*u^6 + 65*u^5 - 130*u^4 - 1885*u^3 + 
        6422*u^2 + 5859*u - 8814)>,
    <0, 1/19968*(-4*u^7 - 13*u^6 + 26*u^5 + 65*u^4 - 1144*u^3 - 
        5083*u^2 - 318*u + 8775)>,
    <0, 1/39936*(-9*u^7 - 65*u^6 + 91*u^5 + 403*u^4 - 2795*u^3 - 
        17459*u^2 + 15385*u + 52065)>,
    <0, 1/19968*(4*u^7 + 13*u^6 - 26*u^5 - 65*u^4 + 1144*u^3 + 2587*u^2
        + 318*u + 3705)>,
    <0, 1/39936*(3*u^7 + 13*u^6 - 65*u^5 - 143*u^4 + 169*u^3 + 2119*u^2
        - 9323*u - 1989)>,
    <0, 1/19968*(7*u^7 - 13*u^6 - 13*u^5 + 143*u^4 + 1781*u^3 - 
        4615*u^2 - 3311*u + 14469)>,
    <0, 1/13312*(3*u^7 + 13*u^6 - 13*u^5 - 195*u^4 + 481*u^3 + 3471*u^2
        + 1129*u - 7449)>
  ];

// Structure constants for the obstruction algebra (checked below)
strconsts := [
    [ 1, 0, 0, 0, 0, 0, 0, 0, 0 ],
    [ 0, 1, 0, 0, 0, 0, 0, 0, 0 ],
    [ 0, 0, 1, 0, 0, 0, 0, 0, 0 ],
    [ 0, 0, 0, 1, 0, 0, 0, 0, 0 ],
    [ 0, 0, 0, 0, 1, 0, 0, 0, 0 ],
    [ 0, 0, 0, 0, 0, 1, 0, 0, 0 ],
    [ 0, 0, 0, 0, 0, 0, 1, 0, 0 ],
    [ 0, 0, 0, 0, 0, 0, 0, 1, 0 ],
    [ 0, 0, 0, 0, 0, 0, 0, 0, 1 ],
    [ 0, 1, 0, 0, 0, 0, 0, 0, 0 ],
    [ 56, -2, -2, 0, 2, 3, -2, 2, 1 ],
    [ 1, 6, 4, 7, 2, -8, 7, -4, 1 ],
    [ -36, 4, 6, 2, 0, -3, -2, -2, 1 ],
    [ -18, 0, -3, -4, -1, 7, 0, 2, 0 ],
    [ 19, 12, 8, 7, 9, -10, 4, -9, -2 ],
    [ -9, 0, 6, 0, 3, -3, -4, -3, -1 ],
    [ 28, -6, -4, -1, -5, 1, -5, 7, 2 ],
    [ -39, 12, 6, 3, -1, -2, 0, -1, 4 ],
    [ 0, 0, 1, 0, 0, 0, 0, 0, 0 ],
    [ 1, -6, -2, 1, -5, -1, -5, -2, 0 ],
    [ -22, -5, -5, -2, -5, 8, -3, 9, 1 ],
    [ -57, -8, -5, 0, -2, 5, 5, 4, -3 ],
    [ -18, 6, 1, 1, 7, -5, 4, 1, -5 ],
    [ -37, -8, -3, 2, -2, -1, 0, 4, 2 ],
    [ 9, -2, 6, -10, 0, -2, 2, 3, 2 ],
    [ -103, 2, -1, -4, -2, 2, 0, 3, -4 ],
    [ 15, -1, 0, -10, 15, 4, -4, -2, 0 ],
    [ 0, 0, 0, 1, 0, 0, 0, 0, 0 ],
    [ -36, 2, -3, 6, -4, 0, -1, -1, -2 ],
    [ -57, -3, -2, -4, 2, 2, -5, 6, 3 ],
    [ -60, -8, -1, -4, 4, -2, 3, 1, 0 ],
    [ 54, -5, -1, -1, -1, 0, 0, 2, -5 ],
    [ -57, -6, 3, -5, -1, 4, 3, 8, 0 ],
    [ 9, 0, -1, -4, -1, 0, 5, 4, 2 ],
    [ -30, -11, -17, 4, -9, 2, -10, -2, 1 ],
    [ -21, 4, -9, 1, 15, 8, 1, 2, -2 ],
    [ 0, 0, 0, 0, 1, 0, 0, 0, 0 ],
    [ -18, 4, 0, -3, -2, -2, -2, 0, 0 ],
    [ -18, -4, -5, 0, -1, 7, -9, 0, 4 ],
    [ 54, -2, 5, -7, 2, 0, 3, -1, 4 ],
    [ 90, 1, 5, 5, 0, -4, 0, -6, -1 ],
    [ -18, -6, -3, -6, 0, 7, -3, 2, 6 ],
    [ -27, -5, -7, 3, -3, 1, 6, -2, 1 ],
    [ 27, 1, 1, -3, -6, 0, -4, -3, 3 ],
    [ 42, -8, 6, 4, -6, 5, 4, 5, -2 ],
    [ 0, 0, 0, 0, 0, 1, 0, 0, 0 ],
    [ 19, -10, -10, 0, -7, 4, -3, 6, -1 ],
    [ -37, 11, 3, 0, 1, 1, 1, 3, 3 ],
    [ -57, 10, 0, 5, 4, -2, 10, -3, 0 ],
    [ -18, 19, 6, 4, 8, -5, 7, -6, -3 ],
    [ -58, 6, 6, 4, -2, -3, 6, 2, 2 ],
    [ -87, -6, 10, -6, 0, -4, -1, -3, -2 ],
    [ -46, 8, -6, 6, -2, -3, 2, -2, 0 ],
    [ -42, 15, 12, 12, 17, -20, 3, -13, 0 ],
    [ 0, 0, 0, 0, 0, 0, 1, 0, 0 ],
    [ -9, -4, -6, -1, -6, 1, 4, 1, 1 ],
    [ 9, 11, 3, 6, 3, -9, 0, -4, -3 ],
    [ 9, 9, 5, 5, 5, -9, 2, -5, -4 ],
    [ -27, -2, -1, -3, 0, 7, 0, 4, 1 ],
    [ -87, 4, -2, -2, -6, 1, -3, 2, 2 ],
    [ -18, 0, 4, -4, 0, -3, 1, -1, 0 ],
    [ 87, -4, 2, 6, 10, -6, 11, -3, -2 ],
    [ -9, 13, 3, 4, 5, -15, 13, -8, -3 ],
    [ 0, 0, 0, 0, 0, 0, 0, 1, 0 ],
    [ 28, 10, 2, 6, 1, 0, -1, -9, 3 ],
    [ -103, -8, -7, 4, 0, 1, -4, 0, 0 ],
    [ -30, -7, 4, -7, -1, -4, 3, 5, 1 ],
    [ 27, -10, -2, -2, -7, 0, 0, 7, -6 ],
    [ -46, 0, -6, -8, 12, 6, 0, -1, 0 ],
    [ 87, 10, 2, -10, 4, 1, 1, 4, 6 ],
    [ 14, -10, -10, -2, -12, 0, -8, 5, 2 ],
    [ 48, 8, 0, -4, 8, 5, 2, 9, -1 ],
    [ 0, 0, 0, 0, 0, 0, 0, 0, 1 ],
    [ -39, 6, 6, -6, 5, -2, 9, 2, 1 ],
    [ 15, -4, 3, -1, -15, -5, -7, 1, 3 ],
    [ -21, -14, 0, 1, -12, -10, 4, -7, -2 ],
    [ 42, -11, -15, -8, -6, 2, -8, -1, -2 ],
    [ -42, -9, -3, -3, -16, 4, -3, -1, 3 ],
    [ -9, 1, -15, 13, -10, 6, -8, 4, 0 ],
    [ 48, -16, -3, -10, -10, -7, 5, 9, -4 ],
    [ 246, 18, 3, 9, 3, -3, 6, -12, -9 ]
];

function AlgebraToVectorSpace(basR)
  V := VectorSpace(Q,9);
  bb := [&cat[Eltseq(x): x in b]: b in basR];
  psi := hom<V->V|Matrix(9,9,bb)>;
  mapRtoV := func<r|V!(&cat[Eltseq(x): x in r]) @@ psi>;
  return V,mapRtoV;
end function;

function TraceMap(rr)
  x1 := rr[1] + Trace(rr[5]);
  x2 := rr[2] + rr[3] + rr[4] + Trace(rr[6]);
  return R!<x1,x2>;
end function;

function TensorMap(r1,r2)
  a1,b1 := Explode(r1);
  a2,b2 := Explode(r2);
  bb1 := sigma(b1);
  bb2 := sigma(b2);
  c := iota10(b1)*iota01(b2);
  return elt<RR|a1*a2,b1*a2,a1*b2,bb1*bb2,b1*bb2,c>;
end function;

// Computing the structure constants
BB := [TensorMap(basR[i],basR[j]):i,j in [1..9]];
V,mapRtoV := AlgebraToVectorSpace(basR);
cc := [ mapRtoV(TraceMap(<eps[i]*rho[i]*bb[i] : i in [1..6]>)): bb in BB];
// Creating the obstruction algebra
A := AssociativeAlgebra< V | cc : Check := true >;
// Checking structure contstants are as above.
assert strconsts eq [Eltseq(A.i*A.j): i,j in [1..9]];

assert Identity(A) eq A.1;
assert A.2*A.2 eq 56*A.1 - 2*A.2 - 2*A.3 + 2*A.5 + 3*A.6 
                                        - 2*A.7 + 2*A.8 + A.9; 
assert A.2*A.3 eq A.1 + 6*A.2 + 4*A.3 + 7*A.4 + 2*A.5 - 8*A.6 
                                        + 7*A.7 - 4*A.8 + A.9; 
assert A.2*A.4 eq -36*A.1 + 4*A.2 + 6*A.3 + 2*A.4 
                                        - 3*A.6 - 2*A.7 - 2*A.8 + A.9; 
assert A.2*A.5 eq -18*A.1 - 3*A.3 - 4*A.4 - A.5 + 7*A.6 + 2*A.8; 
assert A.2*A.6 eq 19*A.1 + 12*A.2 + 8*A.3 + 7*A.4 + 9*A.5 - 10*A.6 
                                        + 4*A.7 - 9*A.8 - 2*A.9; 
assert A.2*A.7 eq -9*A.1 + 6*A.3 + 3*A.5 - 3*A.6 - 4*A.7 - 3*A.8 - A.9; 
assert A.2*A.8 eq 28*A.1 - 6*A.2 - 4*A.3 - A.4 - 5*A.5 + A.6  
                                        - 5*A.7 + 7*A.8 + 2*A.9; 
assert A.2*A.9 eq -39*A.1 + 12*A.2 + 6*A.3 + 3*A.4 - A.5 
                                        - 2*A.6 - A.8 + 4*A.9; 

// Printing out some structure constants
// S<a1,a2,a3,a4,a5,a6,a7,a8,a9> := PolynomialRing(Q,9);
// for i in [2..2] do
// for j in [2..9] do
//   printf "a%o*a%o = %o \\\\ \n",i,j,&+[Eltseq(A.i*A.j)[k]*S.k: k in [1..9]];
// end for; end for;
// print [MinimalPolynomial(x): x in Basis(A)];

pp := Decomposition(OL,7);
assert exists(p){p[1] : p in pp | Valuation(alpha,p[1]) mod 3 eq 1};
assert exists(q){p[1] : p in pp | Valuation(alpha,p[1]) mod 3 eq 2};
assert IsPower(alpha*(p*q^2)^(-1),3);
assert Norm(p) eq 7^3;
assert Norm(q) eq 7^3;

assert Abs(Discriminant(OL)) eq 2^4*3^7*41^4;
assert Abs(MyDiscriminant(A)) eq 2^4*3^16*7^6*41^4;
// print [MinimalPolynomial(b): b in Basis(A)];

// A basis for a maximal order (checked below)
mat0 := Matrix(9,9,[
    [ 1/3, 0, 0, 0, 0, 2/3, 2/3, 2, 17/3 ],
    [ 0, 1/246, 0, 5/123, 11/82, 130/123, 37/246, 713/246, 10201/246 ],
    [ 0, 0, 1/246, 16/123, 0, 12/41, 139/246, 359/246, 5843/123 ],
    [ 0, 0, 0, 1/3, 0, 2/3, 0, 4/3, 119/3 ],
    [ 0, 0, 0, 0, 1/3, 2/3, 0, 7/3, 22/3 ],
    [ 0, 0, 0, 0, 0, 1, 0, 2, 34 ],
    [ 0, 0, 0, 0, 0, 0, 1, 0, 32 ],
    [ 0, 0, 0, 0, 0, 0, 0, 1, 16 ],
    [ 0, 0, 0, 0, 0, 0, 0, 0, 1 ]
]);
bb0 := [&+[mat0[i,j]*A.j: j in [1..9]]: i in [1..9]];

// Computing maximal order (and repeating until the answer agrees 
// with the one I got earlier)
while true do 
  bb := Basis(MaximalOrder(A));  
  if bb eq bb0 then break; end if;
  delete A`MaximalOrder; // stop answer being cached
end while;

assert Abs(MyDiscriminant(A:bas:=bb)) eq 3^6*7^6;
print "Maximal order has discriminant 3^6*7^6";

// Rewriting as a norm equation
mat := LLL(Matrix(Q,9,9,[Eltseq(b): b in bb]));

bb1 := [A!mat[i]: i in [1..9]];
assert Abs(MyDiscriminant(A:bas:=bb)) eq 3^6*7^6;
assert [MinimalPolynomial(b): b in bb1] eq [ 
   X^3 + 2*X + 1, X^3 - 7*X - 2, X^3 + 3*X + 12, X^3 - X + 2, 
   X^3 - 7*X - 26, X^3 + 2*X - 1, X^3 + 5, X^3 - 7*X + 2, 
   X^3 - X^2 - X + 2];

uu := (1/246)*(A![ -82, -21, -32, -4, -37, 26, 23, 25, 13 ]);
vv := (1/246)*(A![ 0, 256, 225, -80, 166, -224, -335, -81, 104 ]);

poly := X^3 + X^2 - 4*X + 1;
assert MinimalPolynomial(uu) eq poly;
assert Discriminant(poly) eq 13^2;

L<u> := NumberField(poly);
sigma := hom<L->L|u^2 + u - 3>;
assert MinimalPolynomial(u) eq MinimalPolynomial(sigma(u));
assert vv*uu*vv^(-1) eq uu^2 + uu - 3;
assert vv^3 eq 3^2*5*7^2*(A!1);
assert #Decomposition(L,3) eq 1;
assert #Decomposition(L,7) eq 1;
assert NormEquation(L,3^2*5*7^2) eq false;
