
/////////////////////////////////////////////////////////////////
//                                                             //
//           Explicit 5-descent on elliptic curves             //
//                       T.A. Fisher                           //
//                                                             //
/////////////////////////////////////////////////////////////////

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

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

printf "\nSection 6. An example of full 5-descent\n\n";

Q := Rationals();
P<X> := PolynomialRing(Q);

function MyRowReduce(mat)
  L := PureLattice(Lattice(mat));
  return LLL(Matrix(Q,[Eltseq(x): x in Basis(L)]));
end function;

function MyKernelMatrix(mat)
  return MyRowReduce(KernelMatrix(mat));
end function;

function ActionOnSubspace(mat1,mat2)
  K := BaseRing(mat2);
  n1 := Nrows(mat1);
  n2 := Ncols(mat1);
  assert Nrows(mat2) eq n2;
  assert Ncols(mat2) eq n2;
  U := VectorSpace(K,n1);
  V := VectorSpace(K,n2);
  mat1 := ChangeRing(mat1,K);
  map := hom<U->V|mat1>;
  return Matrix(K,n1,n1,[Eltseq((mat1[i]*mat2) @@ map): i in [1..n1]]);
end function;

function AlgebraHomomorphism(mats1,mats2,dimhom)
  m := #mats1;
  n1 := Nrows(mats1[1]);
  n2 := Nrows(mats2[1]);
  AA := [Matrix(Q,n1,n2,[<i,j,1>]): j in [1..n2],i in [1..n1]];
  bigmat := ZeroMatrix(Q,n1*n2,0);
  for k in [1..m] do
    newmat := Matrix(Q,n1*n2,n1*n2,[Eltseq(mats1[k]*A-A*mats2[k]): A in AA]);
    bigmat := HorizontalJoin(bigmat,newmat);
    if k gt 1 and Rank(bigmat) eq n1*n2-dimhom then break; end if;
  end for;
  km := MyKernelMatrix(bigmat);
  pi := Matrix(Q,n1,n2,Eltseq(km[1]));
  assert forall{i : i in [1..m]| mats1[i]*pi eq pi*mats2[i]};
  return pi;
end function;

function Trivialisation(A,x)
  nn := Dimension(A);
  assert exists(n){n : n in [3,5] | n^2 eq nn};
  _,mm := MatrixAlgebra(A);
  MM := [mm(A.i): i in [1..nn]];
  rmat := Matrix(Q,nn,nn,[Eltseq(x*A.i): i in [1..nn]]);
  matA := MyRowReduce(rmat);
  matB := MyKernelMatrix(rmat);
  rank := Integers()!(Nrows(matA)/n);
  if rank eq 1 then 
    M := [ActionOnSubspace(matA,m): m in MM];
  elif rank eq n-1 then 
    M := [ActionOnSubspace(matB,m): m in MM];
  end if;
  if n eq 5 and rank in {2,3} then 
    if rank eq 2 then 
      M10 := [ActionOnSubspace(matA,m): m in MM];
      M15 := [ActionOnSubspace(matB,m): m in MM];
    else
      M10 := [ActionOnSubspace(matB,m): m in MM];
      M15 := [ActionOnSubspace(matA,m): m in MM];
    end if;
    pi := AlgebraHomomorphism(M15,M10,6);
    if Rank(pi) eq 10 then 
      mat := MyKernelMatrix(pi);
      M := [ActionOnSubspace(mat,m): m in M15];
    else 
      mat := MyRowReduce(pi);
      M := [ActionOnSubspace(mat,m): m in M10];
    end if;
  end if;
  return M;
end function;

function Alt2Matrix(mat)
  L := BaseRing(mat);
  n := Nrows(mat);
  assert Ncols(mat) eq n;
  nn := Binomial(n,2);
  altmats := [Matrix(L,n,n,[<i,j,1>,<j,i,-1>]): i,j in [1..n]| i lt j];
  return Matrix(L,nn,nn,[[B[i,j]: i,j in [1..n]| i lt j] 
     where B is Transpose(mat)*A*mat : A in altmats]);
end function;

function Sym2Matrix(mat)
  L := BaseRing(mat);
  n := Nrows(mat);
  assert Ncols(mat) eq n;
  nn := Binomial(n+1,2);
  symmats := [Matrix(L,n,n,[<i,j,1>,<j,i,1>]): i,j in [1..n]| i le j];
  return Matrix(L,nn,nn,[[B[i,j]: i,j in [1..n]| i le j] 
     where B is Transpose(mat)*A*mat : A in symmats]);
end function;

// If the action of E[5] on C is given by MW then the action of E[5] 
// on a genus one model defining C is given by [MV,MW] for some MV.
// This rountine converts the matrix MW to the transformation [MV,MW].
function MatrixToTransformation(MW)
  n := Nrows(MW);
  assert n in {3,5};
  alpha := Determinant(MW);
  if n eq 3 then 
    tr := <1/alpha,MW>;
  else
    M10 := Alt2Matrix(MW); 
    M15 := Sym2Matrix(MW); 
    mats1 := [Matrix(Q,10,10,[M10[i,j][k]: i,j in [1..10]]): k in [1..24]];
    mats2 := [Matrix(Q,15,15,[M15[i,j][k]: i,j in [1..15]]): k in [1..24]];
    pi := AlgebraHomomorphism(mats2,mats1,6);
    if Rank(pi) eq 10 then 
      mat := MyKernelMatrix(pi);
      mats := [ActionOnSubspace(mat,m): m in mats2];
    else 
      mat := MyRowReduce(pi);
      mats := [ActionOnSubspace(mat,m): m in mats1];
    end if;
    L := BaseRing(MW);
    MV := &+[ChangeRing(mats[i],L)*Basis(L)[i]: i in [1..24]];
    assert Determinant(MV) eq alpha^2;
    tr := <MV,(1/alpha)*MW>; 
  end if;
  return tr;
end function;

// Writes the field generator in terms of xT and yT.
function ReversePolynomial(T)
  R<x,y> := PolynomialRing(Q,2);
  L := Ring(Parent(T));
  d := Degree(L);
  mat := Matrix(1,d,Eltseq(L!1));
  n := 1;
  bb := [R|1];
  while Rank(mat) lt d do
    n +:= 1;
    if n mod 2 eq 0 then 
      b := x^m where m is Integers()!(n/2);
    else
      b := x^m*y where m is Integers()!((n-3)/2);
    end if;
    bb := bb cat [b];
    row := Eltseq(L!Evaluate(b,[T[1],T[2]]));
    mat := VerticalJoin(mat,Matrix(1,d,row));
  end while;
  vec := Vector(d,Eltseq(L.1));
  soln := Solution(mat,vec);
  f := &+[soln[i]*bb[i] : i in [1..#bb]];
  assert Evaluate(f,[T[1],T[2]]) eq L.1;
  return f;
end function;  

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

E := EllipticCurve("6727a1");
print "E =",CremonaReference(E);
assert BadPrimes(E) eq [7,31];
assert TamagawaNumber(E,7) eq 1;
assert TamagawaNumber(E,31) eq 2;

// Setting up the field L
poly := X^12 + 4*X^11 - 6*X^10 - 20*X^9 + 15*X^8 - 303*X^7 
           + 323*X^6 + 303*X^5 + 15*X^4 + 20*X^3 - 6*X^2 - 4*X + 1;
Ldash<u> := NumberField(poly);
vv := 1/3*(2*u^11 + 9*u^10 - 8*u^9 - 46*u^8 + 10*u^7 - 591*u^6 
           + 343*u^5 + 928*u^4 + 331*u^3 + 60*u^2 + 8*u - 9);
PP<Y> := PolynomialRing(Ldash); 
LL<v> := NumberField(Y^2 - vv);
L := AbsoluteField(LL);
OL := MaximalOrder(L);
FracOL := FieldOfFractions(OL);
assert Discriminant(OL) eq 5^23*7^16*31^18;

// A 5-torsion point on E
xT := 1/15*(-372*u^11 - 1364*u^10 + 2635*u^9 + 6355*u^8 - 7440*u^7
  + 115971*u^6 - 159588*u^5 - 43090*u^4 - 8060*u^3 - 4495*u^2 + 2697*u - 3786);
yT := 1/10*(1922*u^11 + 6727*u^10 - 14415*u^9 - 28830*u^8 + 43245*u^7 
  - 611196*u^6 + 926404*u^5 - 33635*u^4 + 48050*u^3 + 9610*u^2 
  - 16337*u - 1922)*v + 1/30*(372*u^11 + 1364*u^10 - 2635*u^9 - 6355*u^8 
  + 7440*u^7 - 115971*u^6 + 159588*u^5 + 43090*u^4 + 8060*u^3 + 4495*u^2 
  - 2697*u + 3786);
T := E(L)![xT,yT];
assert 5*T eq E(L)!0;
f := ReversePolynomial(T);
uu := Evaluate(f,[Eltseq(2*T)[i] : i in [1..2]]);
sigma2 := hom<L->L|uu>;
assert E(L)![sigma2(x): x in Eltseq(T)] eq 2*T;
assert v eq -31*(2*yT + xT)/(xT^2 + 480*xT + 87391);
assert u eq v/sigma2(v);

// An algebraic representative for a 5-Selmer element
alpha := 1/294*(4600*u^11 + 8325*u^10 - 72155*u^9 - 50035*u^8 
  + 289975*u^7 - 1450795*u^6 + 4510595*u^5 - 592350*u^4 - 3962957*u^3 
  - 1755928*u^2 - 811953*u - 191035)*v + 1/294*(158985*u^11 
  + 661975*u^10 - 836070*u^9 - 3280275*u^8 + 1784950*u^7 - 48064875*u^6 
  + 43645605*u^5 + 52498690*u^4 + 14516335*u^3 + 7628705*u^2 + 310520*u 
  - 311257);
beta := 1/294*(452*u^11 + 1935*u^10 - 2186*u^9 - 9743*u^8 + 4070*u^7 
  - 135379*u^6 + 108106*u^5 + 172665*u^4 + 54912*u^3 + 14840*u^2 - 4879*u 
  - 12762)*v + 1/294*(-1983*u^11 - 9082*u^10 + 7240*u^9 + 46137*u^8 
  - 7149*u^7 + 585937*u^6 - 289205*u^5 - 957562*u^4 - 338134*u^3 
  - 139997*u^2 - 62943*u + 7646);
ff := Factorization(ideal<OL|alpha>);
cc := ff[1][1];
assert ff eq [<cc,5>];
assert Norm(cc) eq 76045490831;
assert alpha*sigma2(alpha)^2 eq beta^5;

prec := 400;
SetKantPrecision(L,prec); 
SetKantPrecision(OL,prec);
SetKantPrecision(FracOL,prec);

C_approx<i> := ComplexField(10);
PP := E(C_approx)![ 1996.324513, -87675.66269 ];
QQ := E(C_approx)![ -643.5460756, 321.7730378 - 13079.33254*i ];

// We record gamma to low precision. (This is enough to specify 
// the choice of 5th root.)
C_approx<i> := ComplexField(3);
gamma_approx := Matrix(C_approx,5,5,[ 
  1.0000, -3.9638 + 0.90392*i, 1.3886 - 4.0571*i, 
           1.3886 + 4.0571*i, -3.9638 - 0.90392*i, 
  -0.91911,5.8658 - 2.1851*i,  2.3876 + 1.9576*i, 
           2.3876 - 1.9576*i,  5.8658 + 2.1851*i, 
  -2.2011, 4.4066 + 3.0044*i, -3.5591 - 4.1895*i, 
          -3.5591 + 4.1895*i,  4.4066 - 3.0044*i, 
  -2.1265,-7.1327 - 4.3290*i, -0.28979+ 3.7545*i, 
          -0.28979- 3.7545*i, -7.1327 + 4.3290*i, 
   4.4363, 0.14277- 0.12419*i,-0.96307- 0.47726*i, 
          -0.96307+ 0.47726*i, 0.14277+ 0.12419*i ]);

print "Computing the real and complex embeddings";
xTconj := Conjugates(T[1]);
yTconj := Conjugates(T[2]);
// Matching up the group law on E(C)[5] with the embeddings L -> C
function emb(r,s)
  if [r,s] eq [0,0] then return 0; end if;
  TT := (r mod 5)*PP + (s mod 5)*QQ;
  scores := [Abs(TT[1]-xTconj[j])+Abs(TT[2]-yTconj[j]): j in [1..24]];
  eps,k := Minimum(scores);
  assert eps lt 10^(-2);
  return k;
end function; 

print "Computing gamma to high precision"; 
alphaconj := Conjugates(L!alpha);
gammaconj := [x^(1/5) : x  in alphaconj];
C<i> := Universe(gammaconj);
zeta := Exp(2*Pi(C)*i/5);
pi := Pi(C);
for r,s in [0..4] do
  if [r,s] ne [0,0] then 
    j := emb(r,s);
    gg := [zeta^k*gammaconj[j]: k in [0..4]];
    g_approx := gamma_approx[r+1,s+1];
    mm := [Abs(1-(g/g_approx)): g in gg];
    eps,k := Minimum(mm);
    assert eps lt 10^(-2);
    gammaconj[j] := gg[k];
  end if;
end for;  

// Checking   (i) gamma*sigma(gamma)^2 = beta, and 
//            (ii) gamma is Gal(C/R)-equivariant
betaconj := Conjugates(L!beta);
for r,s in [0..4] do
  if [r,s] ne [0,0] then 
    z := gammaconj[emb(r,s)]*gammaconj[emb(2*r,2*s)]^2/betaconj[emb(r,s)];
    assert Abs(z-1) lt 10^(-200);
    z := ComplexConjugate(gammaconj[emb(r,s)])/gammaconj[emb(r,-s)];
    assert Abs(z-1) lt 10^(-200);
  end if;
end for;

print "Computing a good basis for L (depending on alpha)";

// Making sign switch. (As explained in the paper this is to make the
// basis elements less messy to print out.)
alpha := 1/alpha;
beta := 1/beta;
alphaconj := [1/x: x in alphaconj];
gammaconj := [1/x: x in gammaconj];

// Computing a good basis for L
Lbasis := [L!x : x in Basis(cc)];
M := Matrix(24,24,[Conjugates(x): x in Lbasis]);
Mconj := Matrix(24,24,[ComplexConjugate(x): x in Eltseq(M)]);
D := DiagonalMatrix([Abs(x)^(2/5): x in alphaconj]);
A := M*D*Transpose(Mconj);
A := Matrix(24,24,[Re(A[i,j]): i,j in [1..24]]);
A := A + Transpose(A);
_,cob := LLLGram(A);
Lbasis := [&+[cob[i,j]*Lbasis[j]: j in [1..24]]: i in [1..24]];

assert [LL!Lbasis[i]: i in [1..2]] eq [ 
  1/147*(906*u^11 + 3697*u^10 - 5099*u^9 - 18382*u^8 + 11847*u^7 
  - 274284*u^6 + 271264*u^5 + 284304*u^4 + 51522*u^3 + 31261*u^2 
  - 4247*u - 3174),
  1/294*(640*u^11 + 2621*u^10 - 3565*u^9 - 13051*u^8 + 8154*u^7 
  - 193589*u^6 + 188894*u^5 + 204155*u^4 + 40745*u^3 + 21338*u^2 
  - 5548*u - 2903)*v 
  + 1/294*(-221*u^11 - 943*u^10 + 1135*u^9 + 4972*u^8 - 2330*u^7 
  + 65086*u^6 - 53197*u^5 - 99488*u^4 - 12061*u^3 + 13094*u^2 
  + 5473*u + 4980) ];

// Computing dual basis for L with respect to the trace form
mat := Matrix(Q,24,24,[Trace(b1*b2):b1 in Lbasis,b2 in Basis(L)])^(-1);
Ldualbasis := [L!Eltseq(mat[i]): i in [1..24]];

print "Computing the structure constants";
mat1 := Matrix(C,5,5,[<(i mod 5)+1,i,1>: i in [1..5]]);
mat2 := DiagonalMatrix(C,5,[zeta^i: i in [0..4]]);
hmats := [ZeroMatrix(C,5,5): i in [1..24]];
for r,s in [0..4] do
  if [r,s] ne [0,0] then 
    hmats[emb(r,s)] := zeta^(3*r*s)*mat1^r*mat2^s;
  end if;
end for;
basis := [Conjugates(x): x in Lbasis];  
idC := IdentityMatrix(C,5);
mats := [gammaconj[j]*hmats[j]: j in [1..24]];
mats := [idC] cat [&+[b[j]*mats[j]: j in [1..24]]: b in basis];
mats := [Matrix(5,5,[Re(mat[i,j]): i,j in [1..5]]): mat in mats];
bigmat := Matrix(25,25,[Eltseq(mat): mat in mats])^(-1);
strconsts := [];
for m1 in mats do
  approxsc := [Vector(Eltseq(m1*m2))*bigmat: m2 in mats];
  strconsts cat:= [[Round(Re(x)): x in Eltseq(y)]: y in approxsc];
end for;
V := VectorSpace(Q,25);
A := AssociativeAlgebra< V | strconsts : Check := true >;

// Taking a look at the structure constants
ss := [Abs(x): x in &cat(strconsts)];
assert forall{x : x in ss | IsIntegral(x)};
assert Maximum(ss) eq 448;
assert (&+ss)/#ss eq 22.650048; 

function MyDiscriminant(A:bas := Basis(A))
  _, mm := MatrixAlgebra(A);
  trd := func<x|(1/5)*Trace(mm(x))>;
  trmat := Matrix(25,25,[trd(bas[i]*bas[j]): i,j in [1..25]]);
  return Determinant(trmat);
end function;

disc := Integers()!MyDiscriminant(A);
assert Abs(disc) eq 5^48*7^16*31^18;
assert Abs(Discriminant(OL)) eq 5^23*7^16*31^18;

assert [MinimalPolynomial(A.i): i in [1..3]] eq 
  [ X - 1, X^5 + 435*X^3 + 7315*X^2 + 835*X + 32172, 
    X^5 - 390*X^3 - 4885*X^2 + 17560*X + 1407822 ];

print "Computing a maximal order in CSA of dimension 25";
bb := Basis(MaximalOrder(A));
cob := LLL(Matrix(Q,25,25,[Eltseq(b): b in bb]));
bas := [A!cob[i] : i in [1..25]];
assert Abs(MyDiscriminant(A:bas := bas)) eq 1;

// Random choices are made in MaximalOrder. One time we ran it we obtained
// the following basis. (This is the basis whose minimal polynomials are 
// recorded in the paper.)

cob := (1/1085)*Matrix(Q,25,25,[ 
0, 3, 6, -9, -4, 3, 4, 5, 3, 0, 3, -1, -7, 4, 4, 2, 3, 3, 10, 6, 4, -6, 1, -1,
-1, 0, -1, 4, -4, -2, 4, -5, -9, -5, 1, -11, -4, -11, 7, 4, -2, 4, 0, 11, 2, 4, 
1, -3, -1, -5, 0, 4, -5, 6, 4, -1, 1, -2, -4, 5, 9, 3, -3, 4, 9, 7, -1, 9, 5, 0,
1, -16, 0, -5, 8, 0, -1, -15, 2, -3, 3, -5, 3, -6, 4, -4, -9, 4, 9, 8, 16, -5, 
1, 17, -2, 2, 10, 5, -5, 0, 0, 4, 7, 15, -5, 4, -5, -4, -7, -10, -6, -1, 6, 6, 
6, 3, 4, -6, -1, 7, -5, -9, 2, -2, 15, 0, -8, 8, 1, -7, 10, 3, 6, -3, -1, 3, 
-14, 19, -6, -4, 7, 1, 0, 4, -2, 2, -1, -5, -3, 2, 0, -7, -13, 6, 1, 1, 10, 8, 
5, 11, 4, 2, 6, 6, -2, -5, 1, 4, 10, -21, 10, 7, 13, -5, 7, 0, 4, 3, 7, 16, -4, 
-15, 7, -9, 9, 12, -2, 2, 11, 7, 9, 5, 4, 9, 7, 5, -1, -4, 14, -11, 0, 4, -8, 
-2, -5, -24, -11, 23, 1, 22, -3, -1, 16, -5, 1, -13, -1, -2, 9, -6, 6, 5, 2, 6, 
9, 0, -6, 1, -8, 2, 4, 13, -2, 17, 11, -4, -2, 14, -3, -12, -10, -7, 6, -9, 5, 
14, 11, -10, -5, 21, 0, -21, 9, 12, 5, 1, 4, -3, 1, 14, 10, 8, -9, -10, -2, -8, 
-8, -7, 8, -12, 2, 2, -14, 6, 0, 0, 10, -2, -10, -10, -10, 1, 6, 10, 0, 5, -13, 
16, -23, 13, 1, 16, 3, 11, -2, -16, 4, -10, 2, -10, 0, -7, 5, 4, 1, -4, 1, -13, 
18, -1, -6, 2, 9, -17, 5, -6, 2, -14, 0, 13, 0, -4, 9, 8, -12, 0, 3, -4, 8, 6, 
-9, -23, 3, -4, 14, -9, -4, 9, -4, 22, 1, 11, -3, -5, 7, -13, -1, 20, -7, 6, 0, 
2, 7, -7, 14, -8, 2, -8, -12, -3, 5, 7, -13, 13, -28, -17, -8, -1, 7, 1, -3, 4, 
-3, 9, 5, 0, 12, 0, -7, 16, 4, 4, 5, -2, -4, -1, 7, 17, 5, -6, 12, 7, -5, 12, 
-5, 7, -15, -3, -10, 3, 0, 5, -9, -3, 21, 3, -22, 9, 12, -1, 8, 5, -1, -6, 13, 
-12, 10, 6, 16, 7, -5, -10, -7, -4, 3, 0, -5, 6, 7, -12, -13, 14, 0, -2, -18, 
-2, 1, -20, 1, 9, 9, -1, -29, -8, 16, 7, -2, -7, -18, -12, 0, -16, -13, -22, 
-10, 5, -7, 13, 1, -3, 2, -17, -12, 4, 9, -15, 2, -32, -10, -16, -3, -12, -8, 
-7, 0, 0, 9, 14, 12, 6, 6, -1, 12, 12, -17, -28, -5, -6, -1, 4, -1, -1, 4, -3, 
-6, -17, -2, 6, -12, -18, 0, -4, -13, 3, 2, -14, -13, -1, 4, 2, 1, -22, -6, -14,
-8, -1, 15, 4, -6, 11, 22, 3, -20, 17, 11, 0, 23, -8, 9, -14, 17, 0, -6, 5, -2, 
9, -9, -4, -7, 0, -14, -14, 6, 7, -9, -15, 1, 1, 2, -20, 0, -21, -10, 7, -7, 10,
7, 15, 6, 2, -20, 16, -2, -9, -25, 18, 11, 3, 3, -3, -13, -13, 12, -13, -6, 0, 
1, -20, 12, 15, 17, 10, 20, -21, -3, -12, -22, 8, 0, 26, -20, 2, 8, 10, 0, 4, 
-10, 17, 25, 4, -217, -4, -5, -9, 5, 0, 4, -8, 6, -3, -5, 3, 2, -8, -13, -6, 
 -11, -6, 5, 2, -4, 5, -6, -6, 7 ]);

bas := [A!cob[i] : i in [1..25]];
assert Abs(MyDiscriminant(A:bas := bas)) eq 1;

assert [MinimalPolynomial(b): b in bas] eq [ 
  X^2, X^2, X^2, X^4, X^2, X^3, X^2, X^3 - X, X^5 - 2*X^3 + X, X^4 - X^2, 
  X^4 - X^2, X^5 + X^3, X^2, X^4 - X^2, X^3, X^4 - 2*X^2, X^4 - X^2, 
  X^5 - X^3 + X^2 + X, X^5 - X^3 - 4*X^2 + 4*X, X^4 - 2*X^2 - X, 
  X^4 + X^3 - X^2 - X, X^5 - X^3, X^5 - 2*X^3, X^5 - 5*X^2 + X, X^3 + X^2 ];

// for b in bas do
//   mp := P!MinimalPolynomial(b);
//   if not IsIrreducible(mp) then
//     f := Factorization(mp)[1][1];
//     a := Evaluate(f,b);
//     break;
//   end if;
// end for;
// printf "Computing trivialisation\n";
// mats := Trivialisation(A,a);
// print "Checking trivialisation";
// assert forall{[i,j] : i,j in [1..25] |
//   mats[i]*mats[j] eq &+[(A.i*A.j)[k]*mats[k]: k in [1..25]]};

// The trivialisation used in the paper was as follows.

trivdata := [
    [ 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 
    0, 0, 0, 0, 1 ],
    [ 13, -5, -20, -20, -15, -40, -22, 40, 20, 10, -20, -35, 3, 
    -15, -15, -15, 0, 30, 13, 0, 15, 15, -10, -5, -7 ],
    [ -12, 5, 0, 10, 5, -30, -42, 35, 5, 0, -50, -35, 38, 20, 5, 
    -110, -50, 65, 8, 5, 45, 45, -30, 0, 8 ],
    [ 0, 0, 0, -10, -10, -20, -15, 30, 25, 10, -15, -10, 25, 10, 0,
    -70, -45, 45, 5, -10, 5, -20, 10, -5, -15 ],
    [ -59, -30, 45, 5, -20, -10, -14, -5, -15, 25, -90, -55, 51, 
    -10, -10, 45, 30, -15, 11, 20, -115, -50, 80, 10, 11 ],
    [ 18, -10, -30, -15, -15, -55, -22, 70, 10, 10, -20, -30, 23, 
    -15, -20, 170, 90, -125, -22, -10, -50, -50, 10, -10, 3 ],
    [ -11, -10, 30, -10, -20, -10, 9, -25, 15, 20, -15, 0, 4, 0, 
    -10, 90, 65, 25, -1, -40, -135, -90, 95, 0, -1 ],
    [ -30, -10, 20, 5, 0, 85, 30, -25, -5, -20, 65, 30, -10, 5, 
    -25, -90, -45, 30, 5, 60, -55, -40, 40, 10, 5 ],
    [ -106, -55, 50, 5, 35, 90, 44, -45, -15, -30, -90, -45, 44, 
    -5, 30, -35, 0, 5, 24, 15, 15, 5, -10, -5, -6 ],
    [ -5, 5, 5, 5, -5, 50, 25, -30, 5, 5, 25, 20, -5, 10, -5, 70, 
    40, -45, -20, 20, -10, -10, -10, 5, 5 ],
    [ 53, 25, -25, -5, 0, -95, -52, 50, 10, 5, -35, -20, 23, 5, 5, 
    120, 40, -100, -22, -35, 85, 55, -50, -10, -2 ],
    [ 128, 80, -80, 5, -10, -110, -92, 75, -10, -5, 75, 20, -27, 0,
    -20, -10, 25, -5, 8, 15, 90, 50, -75, -5, -17 ],
    [ 26, 5, 5, 0, -25, 10, 26, -50, -5, 20, 50, 25, -34, -5, -20, 
    -15, -5, 25, 1, 15, 10, 20, -20, -5, -19 ],
    [ -89, -60, 65, 0, 25, 135, 96, -115, -5, -20, 0, 5, -9, 0, 20,
    40, 5, -15, -14, -15, -30, -10, -15, -10, 16 ],
    [ -6, 20, 0, 10, 15, 10, -21, -15, -20, -10, 20, 15, -26, 0, 5,
    15, 30, 35, 19, -20, -55, -5, 25, 0, 34 ],
    [ -16, -15, 25, 10, -25, -50, -21, 0, -15, 35, -90, -55, 49, 0,
    5, 30, 30, -10, 9, 0, 30, 10, 0, 10, -21 ],
    [ 49, 35, -5, 0, -25, -90, -51, 20, -5, 45, -45, -15, 29, -5, 
    15, 155, 90, -100, -26, -35, 10, 15, -5, 10, -1 ],
    [ 31, 10, -30, -10, -10, -60, -29, 50, 10, 15, -20, -15, 6, -5,
    0, -135, -70, 100, 6, -20, 10, -5, -20, -15, -14 ],
    [ 28, 10, 5, 0, -30, -40, -12, 0, -5, 25, 5, 10, -2, 0, -15, 
    -35, -20, 45, 18, 20, 25, 5, 5, 0, -32 ],
    [ 69, 35, -30, -10, -30, -65, -26, 15, 10, 25, 25, 20, -16, -5,
    -10, 55, 15, -25, -6, -5, -20, -30, 15, -5, -21 ],
    [ 81, 45, -60, 0, -10, -70, -44, 50, 5, 15, 35, 10, -34, 0, 5, 
    -100, -80, 35, 6, -10, 75, 55, -25, 10, -9 ],
    [ 79, 50, -35, -5, -30, -165, -101, 70, 0, 40, -70, -45, 24, 
    -5, 5, 60, 30, 10, 14, -40, 35, 15, -10, -10, -16 ],
    [ 11, 5, -35, -5, 5, -60, -44, 75, 10, -5, -30, -30, 16, 0, -5,
    -15, 0, -25, -9, -5, 15, 25, -35, 0, 26 ],
    [ -40, -25, 20, 10, -5, 75, 50, -30, -20, -5, 25, 20, -10, -5, 
    -20, 40, 25, -5, -20, 10, -70, -45, 40, 10, 20 ],
    [ 7, 15, -10, 10, 5, -25, -23, 15, -25, -5, -40, -15, 12, -5, 
    5, 45, 25, -20, 7, 15, 0, -5, 25, 10, -3 ]
];

mats := [Matrix(Q,5,5,x) : x in trivdata];

print "Checking trivialisation";
assert forall{[i,j] : i,j in [1..25] |
  mats[i]*mats[j] eq &+[(A.i*A.j)[k]*mats[k]: k in [1..25]]};

// Action of E[5] on the curve C is now given by
MW := &+[Ldualbasis[i]*mats[i+1]: i in [1..24]];
assert IsPower(Determinant(MW)/alpha,5);

print "Computing the action of 5-torsion on the space of equations";
trans := MatrixToTransformation(MW);

print "Solving for the Heisenberg invariant models";
I := IdentityMatrix(L,50);
basis := [GenusOneModel(5,Eltseq(I[i])): i in [1..50]];
ans := [];
bigmat := ZeroMatrix(Q,50,0);
_,tr := IsTransformation(5,<trans[1],trans[2]>);
_,tr1 := IsTransformation(5,<Transpose(trans[1]),Transpose(trans[2])>);

for k in [1..50] do
  vv := Vector(ModelToSequence(tr1*basis[k])) - I[k];
  newmat := Matrix(Q,50,24,[Eltseq(x): x in Eltseq(vv)]);
  bigmat := HorizontalJoin(bigmat,newmat);
  rk := Rank(bigmat);
  if k gt 1 and rk eq 48 then break; end if;
end for;
mat := MyKernelMatrix(bigmat);
assert Nrows(mat) eq 2;
pencil := [GenusOneModel(5,Eltseq(mat[i])): i in [1..2]];

print "Extracting the required model from the pencil";
model := pencil[1];
c4,c6,Delta := Invariants(model);
DD,cc4,cc6 := HessePolynomials(5,1,[c4,c6]);
ff := cc4^3 - jInvariant(E)*Delta*DD^5;
f := Evaluate(ff,[X,1]);
rts := Roots(f);
if Degree(f) lt TotalDegree(ff) then 
  assert #rts eq 0;
else
  assert #rts eq 1;
  model := rts[1][1]*model + Hessian(model);
end if;
dd := RationalGCD(Eltseq(model));
model0 := GenusOneModel(5,[x/dd: x in Eltseq(model)]);

// Checking model is fixed by trans
modelL := ChangeRing(model0,L);
assert tr*modelL eq modelL;

// printf "\nMinimising and Reducing\n";
// model1 := Reduce(Minimise(model0));

R<x1,x2,x3,x4,x5> := PolynomialRing(model);
model1 := GenusOneModel(Matrix(5,5,[ 
  0, -x1 + x2 + x3, x1 + 3*x2 + x4, -2*x2 + x3 + x5, 2*x2 - 2*x3 + x5, 
  x1 - x2 - x3, 0, -x1 - x2 - x3 + x5, x2 - x4 + x5, -x2 + x3 + x4, 
  -x1 - 3*x2 - x4, x1 + x2 + x3 - x5, 0, -x3 + x5, -x1 + x3 - x5,
  2*x2 - x3 - x5, -x2 + x4 - x5, x3 - x5, 0, x4, 
  -2*x2 + 2*x3 - x5, x2 - x3 - x4, x1 - x3 + x5, -x4, 0 ]));

subst := [ -x1 + 2*x2 + x3 - 2*x4 + x5, 
  -x1 + x2 + x3 - x4 + x5, x1 - x2 + x4, x5, x4 ];
quads0 := [Evaluate(f,subst): f in Equations(model0)];
quads1 := Equations(model1);
assert ideal<R|quads0> eq ideal<R|quads1>;

printf "%o\n\n",model1;

assert cInvariants(model1) eq cInvariants(E);
P2 := ProjectiveSpace(Rationals(),2);
C := Image(map<Curve(model1)->P2|[x1,x2,x3]>);
assert forall{p : p in BadPrimes(E) | IsLocallySoluble(C,p:Smooth)};

/*

Magma V2.17-3     Thu Mar 31 2011 16:56:19 on coot     [Seed = 1100523644]
Type ? for help.  Type <Ctrl>-D to quit.
> load "fivedesc.m";
Loading "fivedesc.m"

Section 6. An example of full 5-descent

E = 6727a1
Computing the real and complex embeddings
Computing gamma to high precision
Computing a good basis for L (depending on alpha)
Computing the structure constants
Computing a maximal order in CSA of dimension 25
Checking trivialisation
Computing the action of 5-torsion on the space of equations
Solving for the Heisenberg invariant models
Extracting the required model from the pencil
Genus one model of degree 5 defined over Rational Field given by the matrix 
[0   -x1 + x2 + x3   x1 + 3*x2 + x4   -2*x2 + x3 + x5   2*x2 - 2*x3 + x5]
[x1 - x2 - x3   0   -x1 - x2 - x3 + x5   x2 - x4 + x5   -x2 + x3 + x4]
[-x1 - 3*x2 - x4   x1 + x2 + x3 - x5   0   -x3 + x5   -x1 + x3 - x5]
[2*x2 - x3 - x5   -x2 + x4 - x5   x3 - x5   0   x4]
[-2*x2 + 2*x3 - x5   x2 - x3 - x4   x1 - x3 + x5   -x4   0]

> 
*/
