/**********************************************************
 * Some bits and pieces necessary for obtaining cubics    *
 * representing elements in the 3-Selmer group of an      *
 * elliptic curve.                                        *
 *                                                        *
 * Michael Stoll                                          *
 * started 2002-04-28                                     *
 **********************************************************/

// TO DO:
//
// * Rewrite MaximalOrder to really give a maximal order in all cases

declare verbose ThreeDescent, 3;
declare verbose CSAMaximalOrder, 3;

intrinsic Conjugates(x::FldRatElt) -> SeqEnum
{}
  return [RealField() | x];
end intrinsic;

// A function that splits an affine algebra into fields and
// finds an integral basis if over the rationals
function SplitAlgebra(A, pol, gen, flag)
  vprintf ThreeDescent, 1: "  SplitAlgebra -- splitting algebra of dim = %o\n",
                           Dimension(A);
  assert Degree(pol) eq Dimension(A);
  assert Evaluate(pol, gen) eq 0;
  n := Degree(pol);
  V, ma := VectorSpace(A);
  K := BaseRing(A);
  if K cmpeq Rationals() then
    // do some basic reduction
    lcm := LCM([Denominator(c) : c in Coefficients(pol)]);
    pol := lcm^(Degree(pol))*Evaluate(pol, Parent(pol).1/lcm);
    gen := lcm*gen;
    gcd := GCD([Integers() | Coefficient(pol, i) : i in [0..Degree(pol)-1]]);
    primes := PrimeDivisors(gcd);
    for p in primes do
      v := Min([Floor(Valuation(c, p)/i) : i in [1..Degree(pol)] | c ne 0
                 where c := Coefficient(pol, Degree(pol)-i) ]);
      pol := p^(-v*Degree(pol))*Evaluate(pol, Parent(pol).1*p^v);
      gen := p^-v*gen;
    end for;
  end if;
  if IsIrreducible(pol) then // IsField sometimes leads to ``Monomial overflow''
    vprintf ThreeDescent, 3: "    Characteristic polynomial is irreducible.\n";
    K2 := ext<K | pol>;
    if flag and K cmpeq Rationals() then
      vprintf ThreeDescent, 2:
              "    Computing integral basis of degree %o number field...\n",
              Degree(pol);
      KO := OptimisedRepresentation(K2);
      vprintf ThreeDescent, 2:
              "    Using better field description, given by\n      %o\n",
              DefiningPolynomial(KO);
      flag, iso := IsIsomorphic(K2, KO);
      assert flag;
      basK := [b @@ iso : b in IntegralBasis(KO)];
      basA := [ma(Evaluate(Parent(pol)!Eltseq(K2!b), gen)) : b in basK];
    else
      if flag then
        hKA := hom<V -> V | [ma(gen^j) : j in [0..n-1]]>;
        basK := [K2!Eltseq(V.j @@ hKA) : j in [1..n]];
      end if;
      basA := Basis(V);
    end if;
    mb := hom<V -> V | basA>;
    m := ma*Inverse(mb);
    if flag then
      // define msplit: A --> K2 and msplitinv: K2 --> A
      msplit := func<a | <&+[m(a)[i]*basK[i] : i in [1..n]]>>;
      msplitinv := func<tup | Evaluate(Parent(pol)!Eltseq(tup[1]), gen)>;
      vprintf ThreeDescent, 1: "  SplitAlgebra -- done.\n";
      return V, m, msplit, msplitinv;
    else
      vprintf ThreeDescent, 1: "  SplitAlgebra -- done.\n";
      return V, m, _, _;
    end if;
  else
    // split A into fields
    vprintf ThreeDescent, 2: "    Factoring characteristic polynomial...\n";
    fact := Factorisation(pol);
    assert forall{a : a in fact | a[2] eq 1} and #fact ge 2;
    // set up basis
    genbas := [ma(gen^i) : i in [0..n-1]];
    genh := hom<V -> V | genbas>;
    fact := [a[1] : a in fact];
    bas2 := [V | &cat[[Coefficient(m, j) : j in [0..Degree(f)-1]]
                        where m := Parent(f).1^i mod f
                        : f in fact] : i in [0..n-1]];
    genh2 := hom<V -> V | bas2>;
    hh := Inverse(genh)*genh2;
    parts := [&+[Integers() | Degree(f) : f in fact[1..i]] : i in [0..#fact]];
    zero := <>;
    for f in fact do
      Kf := ext<K | f>;
      Append(~zero, Kf!0);
    end for;
    // define msplit : A -> K21 x K22 x ... and its inverse msplitinv
    msplit := function(a)
                b := Eltseq(hh(ma(a)));
                res := <>;
                for i := 1 to #zero do
                  Append(~res, Parent(zero[i])!(b[parts[i]+1..parts[i+1]]));
                end for;
                return res;
              end function;
    msplitinv := function(tup)
                   s := &cat[Eltseq(Parent(zero[i])!tup[i]) : i in [1..#zero]];
                   return (V!s) @@ hh @@ ma;
                 end function;
    bases := <>;
    for i := 1 to #zero do
      if flag and K cmpeq Rationals() then
        vprintf ThreeDescent, 2:
                "    Computing integral basis of degree %o number field...\n",
                Degree(fact[i]);
        Append(~bases, IntegralBasis(Parent(zero[i])));
      else
        Append(~bases, [Kf.1^j : j in [0..Degree(fact[i])-1]]
                         where Kf := Parent(zero[i]));
      end if;
    end for;
    bas := [V | ];
    for i := 1 to #bases do
      over := function(b)
                zo := zero; zo[i] := b; return msplitinv(zo);
              end function;
      bas cat:= [ma(over(b)) : b in bases[i]];
    end for;
    assert sub<V | bas> eq V;
    mb := hom<V -> V | bas>;
    m := ma*Inverse(mb);      
    vprintf ThreeDescent, 1: "  SplitAlgebra -- done.\n";
    return V, m, msplit, msplitinv;
  end if;
end function;

// A workaround for broken WeilPairing function
function wp3(pt1, pt2)
  E := Curve(pt1);
  K := Parent(pt1[1]);
  assert 3*pt1 eq E!0 and 3*pt2 eq E!0;
  // trivial cases
  if pt1 eq E(K)!0 or pt2 eq E(K)!0 then return K!1; end if;
  if pt1[1] eq pt2[1] then return K!1; end if;
  // nontrivial case
  a1,a2,a3,a4,a6 := Explode(aInvariants(E));
  // set up tangent lines  y = li*x + mi,  i = 1,2
  x1 := pt1[1]; y1 := pt1[2];
  x2 := pt2[1]; y2 := pt2[2];
  pt3 := pt2-pt1;
  x3 := pt3[1]; y3 := pt3[2];
  l1 := (3*x1^2 + 2*a2*x1 + a4 - a1*y1)/(2*y1 + a1*x1 + a3);
  m1 := y1 - l1*x1;
  l2 := (3*x2^2 + 2*a2*x2 + a4 - a1*y2)/(2*y2 + a1*x2 + a3);
  m2 := y2 - l2*x2;
  y1i := -a1*x1 - a3 - y1;
  e := (y3-l1*x3-m1)*(y1-l2*x1-m2)/(y1i-l1*x1-m1)/(y1i-l2*x1-m2);
  assert e^3 eq 1;
  return e;
end function;

// first step: Get the algebras, plus the tensor, Delta and trace maps:
// returns  
// (1) A
// (2) A tensor A (as vector spaces)
// (3) a tuple consisting of:
//     (a) the tensoring map as a map<car<V,V> -> VV>
//     (b) Delta : A -> A tensor A as a linear map
//     (c) Trace : A tensor A -> A as a linear map
//     (d) mult : A x A -> A as a map as above
//     (e) mult : AA x AA -> AA as a map as above
// (4) <xT/yT,1/yT, e>, the coordinates of a generic 3-torsion point
//                      and the Weil pairing element
// (5) <D, r> as functions on 2 variables (point coordinates)
intrinsic SetupAlgebras(E::CrvEll) 
  -> ModTupFld, ModTupFld, Tup, Tup, Tup
{Construct the various algebras and maps associated with a 3-descent on E.}
  // Check that E has the simple form  y^2 = x^3 + ax + b
  vprintf ThreeDescent, 1: "SetupAlgebras for %o\n", E;
  as := aInvariants(E);
  a1, a2, a3, a4, a6 := Explode(as);
  K := FieldOfFractions(BaseRing(E));
  PK := PolynomialRing(K);
  PK2 := PolynomialRing(K, 2);
  f := PK2.2^2 + a1*PK2.1*PK2.2 + a3*PK2.2
        - PK2.1^3 - a2*PK2.1^2 - a4*PK2.1 -a6;
  phi3 := 1/3*PK!DivisionPolynomial(E, 3);
  // Set up A1 = algebra of E[3] \ 0
  vprintf ThreeDescent, 2: "  Defining degree 8 algebra...\n";
  A1<xT,yT> := AffineAlgebra<K, xT,yT | Evaluate(f, [xT, yT]),
                                        Evaluate(phi3, xT)>;
  vprintf ThreeDescent, 3: "  Looking for minimal polynomial of degree 8...\n";
  gen := yT;
  polA1 := MinimalPolynomial(gen);
  while Degree(polA1) ne 8 do
    gen := Random(-2,2)*xT + Random(-2,2)*yT;
    polA1 := MinimalPolynomial(gen);
  end while;
  vprintf ThreeDescent, 3: "    Found %o\n", polA1;
  polA1 *:= 1/LeadingCoefficient(polA1);
  V1, m1, msplit1, msplitinv1 := SplitAlgebra(A1, polA1, gen, true);
  A1IsField := IsIrreducible(polA1);
  V := KSpace(K, 9);
  m := func<pair | pair[1]*V.1 + &+[v1[j]*V.(j+1) : j in [1..8]]
                   where v1 := m1(pair[2])>;
  minv := func<v | <v[1], (&+[v[j+1]*V1.j : j in [1..8]]) @@ m1>>;
  conj := func<v | [v[1]] cat &cat[Conjugates(w[i]) : i in [1..#w]]
                   where w := msplit1(minv(v)[2])>;
  // m maps pairs <k,a>, k in K, a in A1 into V
  // Set up A2 = algebra of pairs of independent 3-torsion points
  vprintf ThreeDescent, 2: "  Defining degree 48 algebra...\n";
  A2<xT1,yT1,xT2,yT2>
    := AffineAlgebra<K, xT1,yT1,xT2,yT2 |
                       Evaluate(f, [xT1, yT1]), Evaluate(f, [xT2, yT2]),
                       Evaluate(phi3, xT1),
                       Evaluate(phi3d, [xT1,xT2])>
     where phi3d := ExactQuotient(Evaluate(phi3,P2.1)-Evaluate(phi3,P2.2),
                                  P2.1-P2.2)
     where P2 := PolynomialRing(K, 2);
  vprintf ThreeDescent, 3: "  Looking for minimal polynomial of degree 48...\n";
  gen := yT1 + 2*yT2;
  polA2 := MinimalPolynomial(gen);
  while Degree(polA2) ne 48 do
    gen := Random(-2,2)*xT1 + Random(-2,2)*xT2
            + Random(-2,2)*yT1 + Random(-2,2)*yT2;
    polA2 := MinimalPolynomial(gen);
  end while;
  // vprintf ThreeDescent, 3: "    Found %o\n", polA2;
  polA2 *:= 1/LeadingCoefficient(polA2);
  V2, m2, msplit2, msplitinv2 := SplitAlgebra(A2, polA2, gen, false);
  A2IsField := IsIrreducible(polA2);
  // a helpful function
  function invert(z)
    pol := MinimalPolynomial(z);
    zinv := -1/Coefficient(pol, 0)
              * &+[Coefficient(pol,i)*z^(i-1) : i in [1..Degree(pol)]];
    assert z*zinv eq 1;
    return zinv;
  end function;
  // represent AA as K x A1 x A1 x A1 x A1 x A2, corresponding to orbits of
  // (O,O), (O,T), (T,O), (T,T), (T,-T), (T1,T2), resp.
  // Basis elements 1, 2..9, 10..17, 18..25, 26..33, 34..81
  VV := KSpace(K, 81);
  mm := func<tup | tup[1]*VV.1
                    + &+[v2[j]*VV.(j+1) : j in [1..8]]
                    + &+[v3[j]*VV.(j+9) : j in [1..8]]
                    + &+[v4[j]*VV.(j+17) : j in [1..8]]
                    + &+[v5[j]*VV.(j+25) : j in [1..8]]
                    + &+[v6[j]*VV.(j+33) : j in [1..48]]
                   where v2 := m1(tup[2])
                   where v3 := m1(tup[3])
                   where v4 := m1(tup[4])
                   where v5 := m1(tup[5])
                   where v6 := m2(tup[6])>;
  mminv := func<vv | <vv[1],
                      (&+[vv[j+1]*V1.j : j in [1..8]]) @@ m1,
                      (&+[vv[j+9]*V1.j : j in [1..8]]) @@ m1,
                      (&+[vv[j+17]*V1.j : j in [1..8]]) @@ m1,
                      (&+[vv[j+25]*V1.j : j in [1..8]]) @@ m1,
                      (&+[vv[j+33]*V2.j : j in [1..48]]) @@ m2>>;
  vprintf ThreeDescent, 2: "  Setting up basic elements...\n";
  if A1IsField then
    T := E(A1)![xT,yT];
    lambdaT := -Evaluate(Derivative(f, 1), [xT, yT])
                 / Evaluate(Derivative(f, 2), [xT, yT]);
    // determine xT/yT=zT, 1/yT=wT as elements of V
    wTV := m(<0, 2/(2*yT + a1*xT + a3)>);
    zTV := m(<1, xT/(2*yT + a1*xT + a3)>);
  else
    yTinv := invert(2*yT + a1*xT + a3);
    lambdaT := -Evaluate(Derivative(f, 1), [xT, yT])*yTinv;
    wTV := m(<0, 2*yTinv>);
    zTV := m(<1, xT*yTinv>);
  end if;
  op := hom<A2 -> A2 | xT2, yT2, xT1, yT1>;
  if A2IsField then
    T1 := E(A2)![xT1,yT1];
    T2 := E(A2)![xT2,yT2];
    lambdaT1T2 := (yT2-yT1)/(xT2-xT1);
    Tsum := T1+T2;
    xTsum := Tsum[1]/Tsum[3];
    yTsum := Tsum[2]/Tsum[3];
  else
    dx := xT1 - xT2;
    dxinv := invert(dx);
    lambdaT1T2 := (yT1 - yT2)*dxinv; 
    xTsum := lambdaT1T2^2 + a1*lambdaT1T2 + a2 - xT1 - xT2;
    yTsum := -(yT1 + lambdaT1T2*(xTsum - xT1));
  end if;
  // now for the maps...
  neg := hom<A1 -> A1 | xT, -yT - a1*xT - a3>;
  inj1 := hom<A1 -> A2 | xT1, yT1>;
  inj2 := hom<A1 -> A2 | xT2, yT2>;
  tensor := func<pair1, pair2 | <pair1[1]*pair2[1],
                                 pair1[1]*pair2[2],
                                 pair1[2]*pair2[1],
                                 pair1[2]*pair2[2],
                                 pair1[2]*neg(pair2[2]),
                                 inj1(pair1[2])*inj2(pair2[2])>>;
  tensorV := map<car<V,V> -> VV |
                 pair :-> mm(tensor(minv(pair[1]),minv(pair[2])))>;
  del := hom<A1 -> A2 | xTsum, yTsum>;
  delta := func<pair | <pair[1],
                        pair[2],
                        pair[2],
                        neg(pair[2]),
                        pair[1],
                        del(pair[2])>>;
  deltaV := hom<V -> VV | [mm(delta(minv(b))) : b in Basis(V)]>;
  Mat1, mmat1 := MatrixAlgebra(A1);
  tr1 := func<a | Trace(mmat1(a))>;
  // the trace A2 -> A1 is more complicated to obtain --
  // first find a Delta(A1)-basis of A2
  vprintf ThreeDescent, 2: "  Building the trace map...\n";
  B1 := [V2 | ];
  B2 := [V2 | ];
  S := sub<V2 | >;
  for j := 1 to 6 do
    vprintf ThreeDescent, 3: "    j = %o: ", j;
    repeat // try random linear combinations
      vprintf ThreeDescent, 3: "*";
      if A1IsField then
        v := V2.Random(1,48);
      else
        v := &+[Random(-1,1)*V2.i : i in [1..48]];
      end if;
      Badd := [m2(del(b @@ m1)*(v @@ m2)) : b in Basis(V1)];
      Sadd := sub<V2 | S, Badd>;
    until Dimension(Sadd) eq 8*j;
    vprintf ThreeDescent, 3: "\n";
    Append(~B1, v);
    B2 cat:= Badd;
    S := Sadd;
  end for;
  h := hom<V2 -> V2 | B2>;
  B1A := [b @@ m2 : b in B1];
  tr2 := func<a | (&+[&+[s[i+1,8*i+j+1] : i in [0..5]]*V1.(j+1) : j in [0..7]])
                     @@ m1
                  where s := [m2(a*B1A[i]) @@ h : i in [1..6]]>;
  trace := func<tup | <tup[1]+tr1(tup[5]),
                       tup[2]+tup[3]+neg(tup[4])+tr2(tup[6])>>;
  vprintf ThreeDescent, 3: "    Building the trace map homomorphism...\n";
  traceV := hom<VV -> V | [m(trace(mminv(b))) : b in Basis(VV)]>;
  // multiplication maps
  vprintf ThreeDescent, 2: "  Building multiplication maps...\n";
  mult1V := map<car<V,V> -> V |
                pair :-> m(<p1[1]*p2[1], p1[2]*p2[2]>)
                where p1 := minv(pair[1]) where p2 := minv(pair[2])>;
  mult2V := map<car<VV,VV> -> VV |
                pair :-> mm(<p1[1]*p2[1], p1[2]*p2[2], p1[3]*p2[3],
                             p1[4]*p2[4], p1[5]*p2[5], p1[6]*p2[6]>)
                where p1 := mminv(pair[1]) where p2 := mminv(pair[2])>;
  // determine Weil Pariring
  vprintf ThreeDescent, 2: "  Computing Weil pairing element...\n";
  if A2IsField then
//    e0 := WeilPairing(T1,T2,3);
    e0 := wp3(T1,T2);
    e := mm(<1, A1!1, A1!1, A1!1, A1!1, e0>);
  else
    // Weil pairing is more difficult ...
    xT1i := msplit2(xT1);
    yT1i := msplit2(yT1);
    xT2i := msplit2(xT2);
    yT2i := msplit2(yT2);
    seqs := <>; 
    for i := 1 to #xT1i do
      Kf := Parent(xT1i[i]);
      KfT1 := E(Kf)![xT1i[i], yT1i[i]];
      KfT2 := E(Kf)![xT2i[i], yT2i[i]];
//      ef := WeilPairing(KfT1, KfT2, 3);
      ef := wp3(KfT1, KfT2);
      Append(~seqs, ef);
    end for;
    e0 := msplitinv2(seqs);
    e := mm(<1, A1!1, A1!1, A1!1, A1!1, e0>);
  end if;
  // now the maps D and r
  DV := func<x,y | m(<1, y-yT - lambdaT*(x-xT)>)>;
  if A1IsField then
    rV4 := func<x,y | (y-yT)/(x-xT) - lambdaT>;
  else
    rV4 := func<x,y | (y-yT)*invert(x-xT) - lambdaT>;
  end if;
  if A2IsField then
    rV6 := func<x,y | (y+yTsum+a1*xTsum+a3)/(x-xTsum) - lambdaT1T2>;
  else
    rV6 := func<x,y | (y+yTsum+a1*xTsum+a3)*invert(x-xTsum) - lambdaT1T2>;
  end if;
  rV := func<x,y | mm(<1, A1!1, A1!1, rV4(x,y), x-xT, rV6(x,y)>)>;
  negV := hom<V -> V | [m(<v[1], neg(v[2])>) where v := minv(b)
                         : b in Basis(V)]>;
  return V, VV,
         <tensorV, deltaV, traceV, mult1V, mult2V, m, minv, mm, mminv, conj, neg>, 
         <zTV, wTV, e, lambdaT, lambdaT1T2, xTsum, yTsum, xT, yT, 2*yT+a1*xT+a3>,
         <DV, rV>;
end intrinsic;

// Given an element r_a = <1,1,1, ...> of AA such that there is a = <1,a1>
// in A with  a tensor a = r_a^3 Delta(a), construct the associative
// algebra structure A_a on A.
intrinsic GetCSAlgebra(A::ModTupFld, maps::Tup, coords::Tup, r::ModTupFldElt)
  -> AlgAss, AlgAssElt //, Map
{Construct A_a, given e, r_a.}
  // set up the structure constants
  vprintf ThreeDescent, 1: "GetCSAlgebra ...\n";
  tensor := maps[1];
  mult2 := maps[5];
  trace := maps[3];
  conj := maps[10];
  er := mult2(coords[3], r);
  vprintf ThreeDescent, 2: "  Setting up first algebra ...\n";
  mult := func<x,y | trace(mult2(er, tensor(x, y)))>;
  scs := [[Eltseq(mult(ei,ej)) : ej in Basis(A)] : ei in Basis(A)];
  AlgA := AssociativeAlgebra<A | scs>;
  one := AlgA.1;
  /*
  newsize := &+[Abs(c) : c in &cat &cat scs];
  // do some reduction
  yTconjabs := [Modulus(c) : c in conj(maps[6](<0, coords[9]>))[2..9]];
  alpha := maps[6](<1,rr[4]*rr[5]>) where rr := maps[9](r);
  alphaconjabs := [Modulus(c)^(1/3) : c in conj(alpha)[2..9]];
  ays := [yTconjabs[i]*alphaconjabs[i] : i in [1..8]];
  h0 := hom<A -> A | Basis(A)>;
  repeat
    size := newsize;
    vprintf ThreeDescent, 2: "  Reducing structure constants ...\n";
    basconj := [conj(AlgA.i @@ h0)[2..9] : i in [2..9]];
    Bre := Matrix(RealField(),
                  [[ays[i]*Real(basconj[j,i]) : i in [1..8]] : j in [1..8]]);
    Bim := Matrix(RealField(),
                  [[ays[i]*Imaginary(basconj[j,i]) : i in [1..8]]
                    : j in [1..8]]);
    B := Transpose(Bre)*Bre + Transpose(Bim)*Bim;
    _, T := LLLGram(B);
    AlgA, one, h, newsize
     := ChangeBasisCSAlgebra(AlgA, AlgA.1, 
                             [AlgA.1] cat [&+[T[j,i]*AlgA.(i+1) : i in [1..8]]
                                           : j in [1..8]],
                             true);
    h0 := h*h0;
  until newsize ge size;
  */
  vprintf ThreeDescent, 2: "GetCSAlgebra -- done.\n";
  return AlgA, one; //, h0;
end intrinsic;

intrinsic ChangeBasisCSAlgebra(A::AlgAss, one::AlgAssElt, bas::SeqEnum, 
                               flag::BoolElt)
  -> AlgAss, AlgAssElt, Map, .
{Change structure constants to be with respect to another basis.}
  require Universe(bas) eq A: "Basis elements are not in A.";
  V := Module(A);
  basV := ChangeUniverse(bas, V);
  n := Dimension(A);
  require #bas eq n: "Basis sequence has wrong length.";
  require sub<V | basV> eq V: "Basis elements are not independent.";
  h := hom<V -> V | basV>;
  multtab := [[Eltseq((bas[i]*bas[j]) @@ h) : j in [1..n]] : i in [1..n]];
  newA := AssociativeAlgebra<V | multtab>;
  newone := newA!(one @@ h);
  if flag then
    return newA, newone, h, &+[Abs(c) : c in &cat &cat multtab];
  else
    return newA, newone, h, _;
  end if;
end intrinsic;

intrinsic MaximalOrder(A::AlgAss) -> SeqEnum
{Computes a maximal Z-order in the central simple Q-algebra A.}
  return MaximalOrder(A, []);
end intrinsic;

intrinsic MaximalOrder(A::AlgAss, primelist::[RngIntElt]) -> SeqEnum
{Computes an order that is maximal at the primes in the primelist.}
  require BaseRing(A) cmpeq Rationals(): "Algebra must be over the rationals.";
  vprintf CSAMaximalOrder, 1: "MaximalOrder(primelist = %o)\n", primelist;
  V := Module(A);
  n := Dimension(V);
  flag, nn := IsSquare(n);
  require flag: "Dimension of A must be a square.";
  bas := [A!V.i : i in [1..n]];
  require forall{b : b in  bas | bas[1]*b eq b}:
          "First basis element must be unit.";
  multtab := [[bas[i]*bas[j] : j in [1..n]] : i in [1..n]];
  dps := &join{&join{Seqset(PrimeDivisors(Denominator(c))) : c in Eltseq(e)}
                : e in &cat multtab};
  lcm := LCM(&cat[[Denominator(c) : c in Eltseq(e)] : e in &cat multtab]);
  bas := [bas[1]] cat [lcm*b : b in bas[2..n]];
  Mat, mmat := MatrixAlgebra(A);
  tr0 := map<A -> Rationals() | x :-> 1/nn*Trace(mmat(x))>;
  Q1 := KSpace(Rationals(), 1);
  tr1 := hom<V -> Q1 | [Q1 | [tr0(A!V.i)] : i in [1..n]]>;
  tr := map<A -> Rationals() | x :-> tr1(V!x)[1]>;
  trmat := Matrix([[tr(bas[i]*bas[j]) : j in [1..n]] : i in [1..n]]);
  det := Determinant(trmat);
  if IsEmpty(primelist) then
    primelist := [p : p in PrimeDivisors(Integers()!det)];
  end if;
  primes := Setseq(Seqset(primelist) join dps);
  vprintf CSAMaximalOrder, 2: "  New primelist = %o\n", primes;
  if IsVerbose("CSAMaximalOrder", 2) then
    printf "    Valuations of discriminant: ";
    for p in primes do
      printf "<%o, %o> ", p, Valuation(det, p);
    end for;
    printf "\n";
  end if;
  for p in primes do
    vprintf CSAMaximalOrder, 2: "  Maximizing at prime p = %o\n", p;
    F := GF(p);
    VF := KSpace(F, n);
    val := Valuation(det, p) div 2;
    // some useful functions
    lift := func<a | 2*l gt p select l-p else l where l := Integers()!a>;
    function annihilator(idl, bas, h)
      liftA := func<v | A!&+[lift(v[i])*bas[i] : i in [1..n]]>;
      if Dimension(idl) eq 0 then return idl; end if;
      idlbas := [V | liftA(b) : b in Basis(idl)];
      L := Lattice(Matrix([V | p*b : b in bas] cat idlbas));
      BL := Basis(L);
      hBL := hom<V -> V | BL>;
      multtabL := [Matrix([VF!(((A!V!BL[i])*(A!V!BL[j])) @@ hBL)
                            : j in [1..n]])
                    : i in [1..n]];
      // find annihilator
      ann := &meet[Kernel(m) : m in multtabL];
      annB := [VF | &+[b[i]*BL[i] : i in [1..n]] @@ h : b in Basis(ann)];
      return sub<VF | annB>;
    end function;
    ind := function(b)
             i := 1; while b[i] eq 0 do i +:= 1; end while;
             assert b[i] eq 1;
             return i;
           end function;
    function newbasis(ann, bas)
      liftA := func<v | A!&+[lift(v[i])*bas[i] : i in [1..n]]>;
      newB := [<ind(b), V!liftA(b)> : b in Basis(ann)];
      if newB[1,1] eq 1 then
        i := 2; while Basis(ann)[1,i] eq 0 do i +:= 1; end while;
        b1 := 1/Basis(ann)[1,i]*Basis(ann)[1];
        newB := Append([newB[i] : i in [2..#newB]], <i, V!liftA(b1)>);
      end if;
      if IsVerbose("CSAMaximalOrder", 2) then
        printf "      Adding new basis elements\n";
        for b in newB do printf "        %o/%o\n", b[2], p; end for;
      end if;
      bas1 := bas;
      for b in newB do
        bas1[b[1]] := 1/p*b[2];
      end for;
      return bas1;
    end function;
    // now enter the loop
    while val ne 0 do
      // bas is a Z-basis of the order currently considered
      // liftA lifts an element of this order mod p, expressed on bas,
      //  to an element of A, expressed on its given basis
      liftA := func<v | A!&+[lift(v[i])*bas[i] : i in [1..n]]>;
      h := hom<V -> V | [V | b : b in bas]>;
      multtab := [[(V!(bas[i]*bas[j])) @@ h : j in [1..n]] : i in [1..n]];
      assert LCM(&cat[[Denominator(c) : c in Eltseq(e)] : e in &cat multtab])
               eq 1;
      // find radical
      vprintf CSAMaximalOrder, 2: "    Computing radical ...\n";
      rad := VF;
      for m := 0 to Floor(Log(p, nn)) do
        basrad := Basis(rad);
        basloc := [liftA(b) : b in basrad];
        trmat1 := Matrix(F, [[tr((b1*b2)^(p^m))/p^m : b2 in basloc]
                              : b1 in basloc]);
        rad := sub<VF | [&+[b[i]*basrad[i] : i in [1..#basrad]]
                          : b in Basis(Kernel(trmat1))]>;
      end for;
      dim := Dimension(rad);
      vprintf CSAMaximalOrder, 2: "    Radical has dimension %o\n", dim;
      if dim gt 0 and IsVerbose("CSAMaximalOrder", 3) then
        printf "      Basis of radical:\n";
        for b in Basis(rad) do
          printf "        %o\n", b;
        end for;
      end if;
      // try if radical leads to a larger order
      ann := annihilator(rad, bas, h);
      dimann := Dimension(ann);
      vprintf CSAMaximalOrder, 2: "    Annihilator has dimension %o\n", dimann;
      if dimann gt 0 and IsVerbose("CSAMaximalOrder", 3) then
        printf "      Basis of annihilator:\n";
        for b in Basis(ann) do
          printf "        %o\n", b;
        end for;
      end if;
      if dimann gt 0 then
        bas := newbasis(ann, bas);
        val -:= dimann;
      else // dimann eq 0
        // look for minimal ideals above radical and try them in turn
        vprintf CSAMaximalOrder, 2:
                "    Computing minimal ideals above radical ...\n";
        // set up quotient algebra
        R, red := quo<VF | rad>;
        Rbas := [liftA(b @@ red) : b in Basis(R)];
        scs := [[red(VF!((V!(b1*b2)) @@ h)) : b2 in Rbas] : b1 in Rbas];
        AR := AssociativeAlgebra<R | scs>;
        // find center
        CR := Center(AR);
        vprintf CSAMaximalOrder, 3: "      Center has dimension %o.\n", 
                                    Dimension(CR);
        VCR := Module(CR);
        // find subspace fixed under x |-> x^p
        soc := Kernel(hom<VCR -> VCR | [(CR!b)^p - (CR!b) : b in Basis(VCR)]>);
        dimsoc := Dimension(soc);
        vprintf CSAMaximalOrder, 2:
                "    There are %o minimal ideals.\n", dimsoc;
        // split this
        if Dimension(soc) eq 1 then
          es := [AR!CR!(Basis(soc)[1])];
        else
          work := [soc];
          es := [];
          while not IsEmpty(work) do
            sp := work[#work]; Prune(~work);
            repeat
              el := CR!Random(sp);
              pol := MinimalPolynomial(el);
            until Degree(pol) gt 1;
            elpwrs := [el^i : i in [0..Degree(pol)-1]];
            roots := [Parent(pol).1 - r[1] : r in Roots(pol)];
            spl := [sub<VCR | [VCR | e*CR!c : c in Basis(sp)]
                              where e := &+[Coefficient(q,i)*elpwrs[i+1]
                                             : i in [0..Degree(q)]]
                              where q := ExactQuotient(pol, r)>
                     : r in roots];
            work cat:= [s : s in spl | Dimension(s) gt 1];
            es cat:= [AR!CR!(Basis(s)[1]) : s in spl | Dimension(s) eq 1];
          end while;
        end if;
        // es is list of generators
        idlsR := [sub<R | [R | r1*e*r2 : r1, r2 in Basis(AR)]> : e in es];
        // lift back to A
        idls := [I @@ red : I in idlsR];
        if IsVerbose("CSAMaximalOrder", 3) then
          Icount := 0;
          for I in idls do
            Icount +:= 1;
            printf "      Ideal #%o (dimension %o) has basis\n", 
                   Icount, Dimension(I);
            for b in Basis(I) do
              printf "        %o\n", b;
            end for;
          end for;
        end if;
        flag := false;
        Icount := 0;
        for I in idls do
          Icount +:= 1;
          ann := annihilator(I, bas, h);
          dimann := Dimension(ann);
          vprintf CSAMaximalOrder, 2:
                  "    Annihilator of ideal #%o has dimension %o\n",
                  Icount, dimann;
          if dimann gt 0 and IsVerbose("CSAMaximalOrder", 3) then
            printf "      Basis of annihilator:\n";
            for b in Basis(ann) do
              printf "        %o\n", b;
            end for;
          end if;
          if dimann gt 0 then
            bas := newbasis(ann, bas);
            val -:= dimann;
            flag := true;
            break;
          end if;
        end for;
        if not flag then
          vprintf CSAMaximalOrder, 1: "  Order is maximal at p = %o.\n", p;
          break;
        end if;
      end if; // dimann gt 0
    end while;
  end for;
  return bas;
end intrinsic;

function em(i,j,k,N) // elementary matrix
  return IdentityMatrix(Integers(),N)
          + Matrix(N, [(ii eq i and jj eq j) select k else 0
                          : ii,jj in [1..N]]);
end function;

// set up a list of (monoid) generators
gens := func<N | gs1 cat gs2
                 where gs1 := [em(i,j,1,N) : i,j in [1..N] | i ne j]
                 where gs2 := [em(i,j,-1,N) : i,j in [1..N] | i ne j]>;

function matsize(mat)
  return &+[c^2 : c in Eltseq(mat)];
end function;

function siz(matseq)
  return &+[matsize(mat) : mat in matseq];
end function;

function act(matseq, T)
  return [&+[iT[i,j]*seq0[j] : j in [1..#seq0]] : i in [1..#seq0]]
         where iT := T1^-1
         where seq0 := [Transpose(T1)*mat*T1 : mat in matseq]
         where T1 := ChangeRing(T, CoefficientRing(Universe(matseq)));
end function;

function actelm(matseq, i, j, k)
  seq0 := [AddColumn(AddRow(m, k, i, j), k, i, j) : m in matseq];
  seq0[i] -:= k*seq0[j]; 
  return seq0;
end function;

function act1(matseq, T)
  return [Transpose(T1)*mat*T1 : mat in matseq]
         where T1 := ChangeRing(T, CoefficientRing(Universe(matseq)));
end function;

function act1elm(matseq, i, j, k)
  seq0 := [AddColumn(AddRow(m, k, i, j), k, i, j) : m in matseq];
  return seq0;
end function;

function red0(matseq, gens)
  pols := [<seq, siz(seq), g> where seq := act(matseq, g) : g in gens];
  return Sort(pols, func<a, b | a[2]-b[2]>);
end function;

function red1(matseq)
  // LLL-reduce set of matrices
  n := NumberOfColumns(matseq[1]);
  mat := LLL(Matrix([Eltseq(m) : m in matseq]));
  return [Matrix(n, Eltseq(mat[i])) : i in [1..NumberOfRows(mat)]];
end function;

function red2(matseq)
  id := IdentityMatrix(Rationals(), NumberOfRows(matseq[1]));
  pols := [<seq, siz(seq), id>
            where seq := [ii eq i select matseq[i] + k*matseq[j]
                                  else matseq[ii]
                           : ii in [1..#matseq]]
             : k in [-1,1], j, i in [1..#matseq] | i ne j];
  return pols;
end function;

function red3(matseq, old)
  n := NumberOfRows(matseq[1]);
  T := Matrix(Rationals(),
              [[i eq j select 0
                       else &+[matseq[k,i,l]*matseq[k,j,l]
                                + matseq[k,l,i]*matseq[k,l,j]
                                - matseq[i,k,l]*matseq[j,k,l] : k,l in [1..n]]
                : j in [1..n]] : i in [1..n]]);
  T := -1/&+[matseq[k,i,j]^2 : i,j,k in [1..n]]*T;
  Tr := [Round(c) : c in Eltseq(T)];
  s := &+[Abs(c) : c in Tr]; // vprintf ThreeDescent, 3: "%o ", s;
  while s lt 10 do
    T *:= 2;
    Tr := [Round(c) : c in Eltseq(T)];
    s := &+[Abs(c) : c in Tr]; // vprintf ThreeDescent, 3: "%o ", s;
  end while;
  while s gt 100 do
    T *:= 1/2;
    Tr := [Round(c) : c in Eltseq(T)];
    s := &+[Abs(c) : c in Tr]; // vprintf ThreeDescent, 3: "%o ", s;
  end while;
  Trm := Matrix(n, Tr);
  gens := [<i, j, Sign(Trm[i,j]), Abs(Trm[i,j])> : j, i in [1..n]
                                                 | Trm[i,j] ne 0];
  gens1 := &cat[[<i, j, s> : s in [-1, 1] | s ne Sign(Trm[i,j])]
                  : j, i in [1..n] | i ne j];
  new1 := [<seq, siz(seq), em(g[1], g[2], g[3], n)>
             : g in gens1 | seq notin old
                            where seq := act1elm(matseq, g[1], g[2], g[3])];
  Append(~new1, <seq, siz(seq), IdentityMatrix(Integers(), n)>
                where seq := red1(matseq));
  Sort(~new1, func<a, b | a[2]-b[2]>);
  new1 := new1[1..Min(10, #new1)];
  tr := IdentityMatrix(Integers(), n);
  agenda := [<matseq, siz(matseq), gens, tr, true>];
  while exists(i){i : i in [1..#agenda] | agenda[i,5]} do
    entry := agenda[i];
    agenda[i,5] := false;
    gl := entry[3];
    if not IsEmpty(gl) then
      new := [<seq, siz(seq), gl1, tr*em(g[1], g[2], g[3], n), true>
                 : i in [1..#gl] | seq notin old
               where seq := act1elm(matseq, g[1], g[2], g[3])
               where gl1 := g[4] eq 1 select Remove(gl, i) 
                                      else f()
                                           where f := function()
                                                        gl1 := gl;
                                                        gl1[i,4] -:= 1;
                                                        return gl1;
                                                      end function
               where g := gl[i]];
      Sort(~new, func<a, b | a[2]-b[2]>);
      k := 1;
      for newe in new do
        while k le #agenda and agenda[k,2] lt newe[2] and k le 10 do
          k +:= 1;
        end while;
        if k ge 10 then break; end if;
        Insert(~agenda, k, newe);
      end for;
      if #agenda gt 10 then agenda := agenda[1..10]; end if;
    end if;
  end while;
  // vprintf ThreeDescent, 3: "min = %o\n", agenda[1,2];
  new2 := [<a[1], a[2], a[4]> : a in agenda[1..Min(10, #agenda)]];
  new := new1 cat new2;
  new cat:= [<seq, siz(seq), a[3]> where seq := red1(a[1]) : a in new];
  _, pos := Min([a[2] : a in new]);
  return new[pos];
end function;

function red4(matseq)
  n := NumberOfColumns(matseq[1]);
  Us := [Upper(m) : m in matseq];
  new := [<seq, siz(seq), Transpose(T)>
          where seq := act1(matseq, Transpose(T))
          where _, T := LLLGram(U + Transpose(U))
           : U in Us];
  new cat:= [<seq, siz(seq), a[3]>
             where seq := [Matrix(n, Eltseq(mat[i])) : i in [1..n]]
             where mat := LLL(Matrix([Eltseq(m) : m in a[1]]))
              : a in new];
  _, pos := Min([a[2] : a in new]);
  return new[pos];
end function;

function rand(matseq, n, k)
  new := [<seq, siz(seq), g>
           where seq := act1(matseq, g)
           where g := RandomSLnZ(n, 5, k)
           : j in [1..20]];
  _, pos := Min([a[2] : a in new]);
  return new[pos];
end function;

intrinsic Trivialize(Aa::AlgAss, bas::SeqEnum) -> Map
{Given a split central simple algebra over Q and a basis of a maximal order, 
 find an explicit isomorphism of this maximal order with the matrix algebra.}
  vprintf ThreeDescent, 1: "Trivialize ...\n";
  n := Dimension(Aa);
  require n eq 9: "Dimension of Aa has to be 9.\n";
  _, mm := MatrixAlgebra(Aa);
  trmat := Matrix([[1/3*Trace(mm(bas[i]*bas[j])) : j in [1..n]] : i in [1..n]]);
  require Abs(Determinant(trmat)) eq 1:
          "Aa must be split, and bas must be a basis of a maximal order.\n";
  V := Module(Aa);
  h := hom<V -> V | bas>;
  sss := [[Eltseq((bas[i]*bas[j]) @@ h) : j in [1..n]] : i in [1..n]];
  mmats := [Matrix([[sss[i,j,k] : j in [1..n]] : i in [1..n]]) : k in [1..n]];
  size := siz(mmats);
  vprintf ThreeDescent, 2: "    Size = %o\n", size;
  // first reduce trace form
  vprintf ThreeDescent, 2: "  Reduce trace form ...\n";
  U := Upper(trmat);
  _, T := LLLGram(U + Transpose(U));
  bas := [&+[T[i,j]*bas[j] : j in [1..n]] : i in [1..n]];
  T := Transpose(T);
  mmats := act(mmats, T);
  size := siz(mmats);
  vprintf ThreeDescent, 2: "    Size = %o\n", size;
  // reduce sum-of-squares of forms for x^2 = 0
  PP := PolynomialRing(Rationals(), n);
  vprintf ThreeDescent, 2: "  Reduce forms for x^2 = 0 ...\n";
  // put in new method, acting directly on the matrices
  repeat
    olds := size;
    U := &+[Upper(m + Transpose(m)) : m in mmats];
    U := U + Transpose(U);
    while not IsPositiveDefinite(U) do U *:= 10; end while;
    _, T := LLLGram(U);
    bas := [&+[T[i,j]*bas[j] : j in [1..n]] : i in [1..n]];
    T := Transpose(T);
    mmats := act(mmats, T);
    size := siz(mmats);
    vprintf ThreeDescent, 2: "    Size = %o\n", size;
  until size ge olds;
  vprintf ThreeDescent, 2: "  Done with x^2 = 0 reduction.\n", size;
  //
  // forms := [&+[PP.i*mmats[k,i,j]*PP.j : i,j in [1..n]] : k in [1..n]];  
  // forms, tr := RedForms(forms);
  // mats := [Matrix(Integers(),
  //                [[(i eq j select 2 else 1)*MonomialCoefficient(f, PP.i*PP.j)
  //                   : j in [1..n]] : i in [1..n]])
  //          : f in forms];
  // bas := [&+[bas[j]*tr[j,i] : j in [1..n]] : i in [1..n]];
  mats := [m + Transpose(m) : m in mmats];
  size := siz(mats);
  vprintf ThreeDescent, 2: "    Size (symmetrized) = %o\n", size;
  // now do rest with graph search
  vprintf ThreeDescent, 2: "  Do graph search ...\n";
  m := 10;
  _, sn := IsSquare(n);
  id := IdentityMatrix(Integers(), n);
  // the following function finds an element t /= 0 with t^2 = 0
  function reduce(mats, bas, size)
    current := mats;     // current sequence of quadratic forms (as matrices)
    best := current;     // best so far
    trans := id;         // transformation from best to current
    target := size;      // size(best), to be improved upon
    origsize := size;    // size of original sequence
    steps := 0;          // counts steps
    randcount := 0;      // counts unsuccessful consecutive random steps
    randk := 3;          // affects the size of the random transformations
    while size gt 2*(n+1)*sn do
      vprintf ThreeDescent, 3: "    current = %o, best = %o\n", size, target;
      //
      h := hom<V -> V | bas>;
      sss := [[Eltseq((V!(b1*b2)) @@ h) : b2 in bas] : b1 in bas];
      mats := [Matrix(Integers(), [[sss[i,j,k]+sss[j,i,k] : j in [1..n]]
                                    : i in [1..n]])
                : k in [1..n]];
      assert Lattice(Matrix([Eltseq(m) : m in mats]))
              eq Lattice(Matrix([Eltseq(m) : m in best]));
      //
      new := red3(current, {}); // get best descendant
      steps +:= 1;
      if new[2] lt target then
        // improvement upon best so far
        T := trans*new[3];
        bas := [&+[bas[i]*T[i,j] : i in [1..n]] : j in [1..n]];
        best := new[1];
        target := new[2];
        trans := id;
        current := best;
        size := target;
        randcount := 0;
      elif new[2] lt size then
        // improvement upon previous
        current := new[1];
        size := new[2];
        trans := trans*new[3];
      else
        // try something else
        if exists(t){b : b in bas | b^2 eq 0} then
          vprintf ThreeDescent, 2: "  Solution found.\n";
          return t;
        end if;
        forms := [&+[PP.i*best[k,i,j]*PP.j : i,j in [1..n]] : k in [1..n]];
        if target le 2000 then
          // try find a point
          vprintf ThreeDescent, 2: "  Try to find solution ...\n";
          I := ideal<PP | forms>;
          v := [];
          for i := 1 to 20 do
            I1 := ideal<PP | I, [PP.i - Random(-5,5) : i in [6..9]]>;
            if Dimension(I1) eq 0 then
              v := [w : w in Variety(I1)
                      | exists{j : j in [1..#w] | w[j] ne 0}];
              if not IsEmpty(v) then break; end if;
            end if;
          end for; // i
          if not IsEmpty(v) then
            v := [v[1, i] : i in [1..#v[1]]];
            d := LCM([Denominator(c) : c in v]);
            v := [Integers() | d*c : c in v];
            g := GCD(v);
            v := [ExactQuotient(c, g) : c in v];
            t := &+[v[i]*bas[i] : i in [1..#v]];
            assert t^2 eq 0;
            vprintf ThreeDescent, 2: "  Solution found.\n";
            return t;
          end if; // not IsEmpty(v)
        end if; // target le 2000
        // try reduction of individual matrices
        vprintf ThreeDescent, 2: "  Try individual reduction ...\n";
        new := red4(best);
        // try forms reduction
        if new[2] ge size then
          vprintf ThreeDescent, 2: "  Try quartic forms reduction ...\n";
          forms1, tr := RedForms(forms);
          mats := [Matrix(Integers(),
                    [[(i eq j select 1 else 1/2)
                       *MonomialCoefficient(f, PP.i*PP.j)
                       : j in [1..n]] : i in [1..n]])
                    : f in forms1];
          new := <mats, siz(mats), tr>;
        end if; // new[2] ge size
        if new[2] lt target then
          // improvement upon best so far
          T := new[3];
          bas := [&+[bas[i]*T[i,j] : i in [1..n]] : j in [1..n]];
          best := new[1];
          target := new[2];
          trans := id;
          current := best;
          size := target;
          randcount := 0;
        elif new[2] lt size then
          // improvement upon previous
          current := new[1];
          size := new[2];
          trans := trans*new[3];
        else
          // no improvement
          vprintf ThreeDescent, 2: "  Randomize ...\n";
          randcount +:= 1;
          new := rand(best, n, randcount*randk);
          current := new[1];
          size := new[2];
          trans := new[3];
        end if; // sm lt target
      end if; // new[2] lt target
    end while;
    vprintf ThreeDescent, 2: "  %o reduction steps, average gain: %o\n",
                             steps, 
                             steps eq 0 select 1
                                        else (origsize/(2*(n+1)*sn))^(1/steps);
    assert exists(t){b : b in bas | b^2 eq 0};
    return t;
  end function;
  if not exists(t){b : b in bas | b^2 eq 0} then
    t := reduce(mats, bas, size);
  end if;
  vprintf ThreeDescent, 2: "  Found element t with t^2 = 0 in maximal order.\n";
  iso := ConstructIsomorphism(Aa, Aa.1, t, bas);
  vprintf ThreeDescent, 1: "Trivialize -- done.\n";
  return iso;
end intrinsic;

intrinsic Trivialize1(Aa::AlgAss, bas::SeqEnum) -> Map
{Given a split central simple algebra over Q and a basis of a maximal order, 
 find an explicit isomorphism of this maximal order with the matrix algebra.}
  vprintf ThreeDescent, 1: "Trivialize1 ...\n";
  n := Dimension(Aa);
  require n eq 9: "Dimension of Aa has to be 9.\n";
  _, mm := MatrixAlgebra(Aa);
  trmat := Matrix([[1/3*Trace(mm(bas[i]*bas[j])) : j in [1..n]] : i in [1..n]]);
  require Abs(Determinant(trmat)) eq 1:
          "Aa must be split, and bas must be a basis of a maximal order.\n";
  V := Module(Aa);
  h := hom<V -> V | bas>;
  sss := [[Eltseq((bas[i]*bas[j]) @@ h) : j in [1..n]] : i in [1..n]];
  mmats := [Matrix([[sss[i,j,k] : j in [1..n]] : i in [1..n]]) : k in [1..n]];
  size := siz(mmats);
  vprintf ThreeDescent, 2: "    Size = %o\n", size;
  // first reduce trace form
  vprintf ThreeDescent, 2: "  Reduce trace form ...\n";
  U := Upper(trmat);
  _, T := LLLGram(U + Transpose(U));
  bas := [&+[T[i,j]*bas[j] : j in [1..n]] : i in [1..n]];
  // reduce Tr(x^3)
  PP := PolynomialRing(Rationals(), n);
  F := &+[1/3*Trace(mm(bas[i]*bas[j]*bas[k]))*PP.i*PP.j*PP.k : i,j,k in [1..n]];
  size := &+[c^2 : c in Coefficients(F)];
  vprintf ThreeDescent, 2: "  Size of cubic trace form = %o\n", size;
  vprintf ThreeDescent, 2: "  Reduce cubic trace form ...\n";
  F, T := RedForm(F);
  bas := [&+[bas[j]*T[j,i] : j in [1..n]] : i in [1..n]];
  return Trivialize(Aa, bas);
end intrinsic;

// Given a central simple algebra of dimension 3^2 and an element t
// in it satisfying t^2 = 0, construct an isomorphism with the
// matrix algebra.
intrinsic ConstructIsomorphism(Aa::AlgAss, one::AlgAssElt, t::AlgAssElt, 
                               bas::SeqEnum)
  -> Map, SeqEnum
{Find isomorphism of Aa with Mat(3,K), given t with t^2 = 0.}
  require Aa eq Parent(t): "t must live in Aa.";
  require Dimension(Aa) eq 9: "Aa must be of dimension 9.";
  require t*t eq Aa!0 and t ne Aa!0: "t must satisfy t^2 = 0 and be nonzero.";
  require Universe(bas) eq Aa and #bas eq 9
           and Dimension(sub<Module(Aa) | ChangeUniverse(bas, Module(Aa))>)
                eq 9:
          "bas must be a basis of Aa.";
  V := Module(Aa);
  h := hom<V -> V | bas>;
  th := t @@ h;
  require forall{c : c in Eltseq(th) | Denominator(c) eq 1}:
          "t must be in order generated by bas.";
  th := 1/GCD(ChangeUniverse(Eltseq(th), Integers()))*th;
  for i := 1 to 9 do if t[i] ne 0 then ti := i; break; end if; end for;
  quott := func<u | u[ti]/t[ti]>;
  // first step: set up linear form l such that t*x*t = l(x)*t
  vprintf ThreeDescent, 1: "ConstructIsomorphism: t = %o\n", th;
  lin := [quott(t*b*t) : b in bas];
  g, gs := XGCD(ChangeUniverse(lin, Integers()));
  assert g eq 1;
  e1 := t*Aa!h(V!gs);
  assert e1^2 eq e1;
  e1h := e1 @@ h;
  vprintf ThreeDescent, 2: "  Found idempotent e1 = %o\n", e1h;
  // second step: set up linear form L such that e1*x*t = L(x)*t 
  lin1 := [quott(e1*b*t) : b in bas];
  sol := Solution(Transpose(Matrix(Integers(), [lin, lin1])), Vector([1, 0]));
  e2 := (Aa!h(V!sol))*t;
  e2h := e2 @@ h;
  vprintf ThreeDescent, 2: "  Found idempotent e2 = %o\n", e2h;
  assert e2^2 eq e2 and e1*e2 eq Aa!0 and e2*e1 eq Aa!0;
  e3 := one - e1 - e2;
  e3h := e3 @@ h;
  assert e3*e3 eq e3 and e3*e1 eq Aa!0 and e1*e3 eq Aa!0
          and e3*e2 eq Aa!0 and e2*e3 eq Aa!0;
  vprintf ThreeDescent, 2: "  Found idempotent e3 = %o\n", e3h;
  // now take e1 to [[100][000][000]], e2 to [[000][010][000]],
  // e3 to [[000][000][001]]
  // find other basis elements as intersections of left and right ideals
  lme1 := Matrix(Integers(), [Eltseq((e1*b - b) @@ h) : b in bas]);
  rme1 := Matrix(Integers(), [Eltseq((b*e1 - b) @@ h) : b in bas]);
  lme2 := Matrix(Integers(), [Eltseq((e2*b - b) @@ h) : b in bas]);
  rme2 := Matrix(Integers(), [Eltseq((b*e2 - b) @@ h) : b in bas]);
  lme3 := Matrix(Integers(), [Eltseq((e3*b - b) @@ h) : b in bas]);
  rme3 := Matrix(Integers(), [Eltseq((b*e3 - b) @@ h) : b in bas]);
  B12 := Basis(Kernel(lme1) meet Kernel(rme2));
  assert #B12 eq 1;
  e12 := Aa!h(B12[1]);
  B23 := Basis(Kernel(lme2) meet Kernel(rme3));
  assert #B23 eq 1;
  e23 := Aa!h(B23[1]);
  B31 := Basis(Kernel(lme3) meet Kernel(rme1));
  assert #B31 eq 1;
  e31 := Aa!h(B31[1]);
  for i := 1 to 9 do if e1[i] ne 0 then e1i := i; break; end if; end for;
  e1231 := e12*e23*e31;
  s := e1[e1i]/e1231[e1i];
  assert s in {-1, 1};
  e31 := s*e31;
  e21 := e23*e31;
  e32 := e31*e12;
  e13 := e12*e23;
  assert e12*e23*e31 eq e1 and e13*e32*e21 eq e1;
  basis := [e1,e12,e13, e21,e2,e23, e31,e32,e3];
  Mat := MatrixAlgebra(BaseRing(Aa), 3);
  h := Inverse(hom<V -> V | ChangeUniverse(basis, V)>);
  iso:= map<Aa -> Mat | a :-> &+[v[i]*Basis(Mat)[i] : i in [1..9]]
                                where v := h(a)>;
  mat := [V | Eltseq(iso(Aa.i)) : i in [1..9]];
  hh := Inverse(hom<V -> V | mat>);
  bas1 := [Aa | hh(V.i) : i in [1..9]];
  assert [iso(b) : b in bas1] eq Basis(Mat);
  return iso, bas1;
end intrinsic;

// Now try to approach the curves


intrinsic CoveringCurve1(A::ModTupFld, maps::Tup, coords::Tup, 
                         r::ModTupFldElt) 
  -> SeqEnum
{Gives model of curve as degree 9 curve in P^8.}
  vprintf ThreeDescent, 1: "CoveringCurve1\n";
  P9<[z]> := PolynomialRing(BaseRing(A), 9);
  lT := coords[4];
  lT1T2 := coords[5];
  A1 := Parent(lT);
  yT := A1.2;
  A2 := Parent(lT1T2);
  AA := Codomain(maps[1]);
  tensor := maps[1];
  delta := maps[2];
  del := hom<A1 -> A2 | coords[6], coords[7]>;
  neg := maps[11]; // hom<A1 -> A1 | A1.1, -A1.2>;
  multAA := maps[5];
  m := maps[6]; minv := maps[7];
  mm := maps[8]; mminv := maps[9];
  mons2 := [z[i]*z[j] : j in [i..9], i in [1..9]];
  vprintf ThreeDescent, 2: "  Setting up basis...\n";
  base := [A.1] cat [m(<0,coords[10]*minv(A.(i+1))[2]>) : i in [1..8]];
  // compute r*(z tensor z)
  vprintf ThreeDescent, 2: "  Setting up r*(z tensor z)...\n";
  rzzc := [(i eq j select 1/2 else 1)
            * multAA(r, tensor(base[i], base[j]) + tensor(base[j], base[i]))
            : j in [i..9], i in [1..9]];
  rP1 := mm(<1,1,1,-lT,-A1.1,-lT1T2>);
  rP2 := mm(<0,0,0,1,1,1>);
  vprintf ThreeDescent, 2: "  Setting up equations...\n";
  rhsc1 := [multAA(rP1, delta(base[j])) : j in [1..9]]
             cat [AA | 0 : j in [10..#mons2]];
  rhsc2 := [multAA(rP2, delta(base[j])) : j in [1..9]]
             cat [AA | 0 : j in [10..#mons2]];
  vprintf ThreeDescent, 3: "    eqns5a...\n";
  eqns5a := [[rzzc[j][i] - rhsc1[j][i] : j in [1..#mons2]]
              : i in [26..33]];
  rzzspl := [mminv(c) : c in rzzc];
  assert forall{0 : s in rzzspl | neg(s[5]) eq s[5]};
  assert rhsc2[1][26] ne 0;
  assert forall{0 : j in [2..#mons2] | mminv(rhsc2[j])[5] eq 0};
  vprintf ThreeDescent, 3: "    eqns5...\n";
  eqns5 := [[rhsc2[1][26]*eqns5a[i,j]-rhsc2[1][25+i]*eqns5a[1,j]
              : j in [1..#mons2]]
             : i in [2..#eqns5a]];
  vprintf ThreeDescent, 3: "    eqns6a...\n";
  eqns6a := [mminv(rzzc[j]-rhsc1[j])[6] : j in [1..#mons2]];
  vprintf ThreeDescent, 3: "    eqns4a...\n";
  eqns4a := [del(neg(mminv(Eltseq(rzzc[j]-rhsc1[j]))[4])) : j in [1..#mons2]];
  assert forall{j : j in [1..#mons2] |
                    del(neg(mminv(Eltseq(rhsc2[j]))[4]))
                     eq mminv(Eltseq(rhsc2[j]))[6]};
  vprintf ThreeDescent, 3: "    eqns6...\n";
  etemp := [mm(<0,0,0,0,0,eqns6a[j]-eqns4a[j]>) : j in [1..#mons2]];
  eqns6 := [[etemp[j, i] : j in [1..#mons2]] : i in [34..81] ];
  eqnsmat := Matrix(eqns5 cat eqns6);
  // eliminate unnecessary equations and eliminate z[1]
  vprintf ThreeDescent, 2: "  Eliminating...\n";
  eqnsmat1 := Submatrix(eqnsmat, 1, 1, NumberOfRows(eqnsmat), 9);
  neweqns := Basis(Image(KernelMatrix(eqnsmat1)*eqnsmat));
  assert #neweqns eq 18;
  return [z[1]] cat [&+[b[i]*mons2[i] : i in [1..#mons2]] : b in neweqns];
end intrinsic;

intrinsic CoveringCurve2(eqns1::SeqEnum, iso::Map) -> RngMPolElt
{Convert into cubic form.}
  vprintf ThreeDescent, 1: "CoveringCurve2\n";
  P9<[u]> := PolynomialRing(BaseRing(Universe(eqns1)), 9);
  vprintf ThreeDescent, 2: "  Transferring equations to Mat(3)...\n";
  mat := Matrix([Eltseq(iso(Domain(iso).i)) : i in [1..9]]);
  mati := mat^-1;
  subst := [&+[u[j]*mati[j,i] : j in [1..9]] : i in [1..9]];
  eqns2 := [Evaluate(e, subst) : e in eqns1];
  // do some reduction first.
  mons2 := MonomialsOfDegree(P9, 2);
  eqns := eqns2[2..19];
  M := Matrix([[MonomialCoefficient(e, m) : m in mons2] : e in eqns]);
  vprintf ThreeDescent, 2: "  Computing equation lattice...\n";
  L := PureLattice(Lattice(M));
  vprintf ThreeDescent, 2: "  LLL reduction...\n";
  L := LLL(L);
  eqns := [&+[b[i]*mons2[i] : i in [1..#mons2]] : b in Basis(L)];
  P3<x,y,z> := PolynomialRing(BaseRing(P9), 3);
  vprintf ThreeDescent, 2: "  Computing elimination ideal...\n";
  // Iel := EliminationIdeal(ideal<P9 | eqns, 3*eqns2[1]>, {u[1],u[4],u[7]});
  // B := Basis(Iel);
  // do it with linear algebra instead
  eqns3 := &cat[[u[i]*e : i in [1..9]] : e in eqns]
            cat [u[i]*u[j]*(u[1]+u[5]+u[9]) : j in [i..9], i in [1..9]];
    // cubic forms in ideal
  // want to get rid of anything that contains something besides u[1],u[4],u[7]
  mons1 := [u[1]^i*u[4]^j*u[7]^(3-i-j) : j in [0..3-i], i in [0..3]];
  mons := [m : m in MonomialsOfDegree(P9, 3) | m notin mons1];
  M := Matrix([[MonomialCoefficient(e, m) : m in mons] : e in eqns3]);
  Ker := Kernel(M);
  B := BasisMatrix(Ker); // Basis(Iel);
  vprintf ThreeDescent, 3: "    Kernel dim: %o\n", Dimension(Ker);
  B := Basis(Image(B*Matrix([[MonomialCoefficient(e, m) : m in mons1]
                              : e in eqns3])));
  vprintf ThreeDescent, 3: "    Image dim: %o\n", #B;
  assert #B eq 1;
  cubic := Evaluate(&+[B[1,i]*mons1[i] : i in [1..#mons1]],
                    [x,0,0,y,0,0,z,0,0]);
  // assert #B eq 1;
  // cubic := Evaluate(B[1], [x,0,0,y,0,0,z,0,0]);
  // assert forall{m : m in Monomials(cubic) | TotalDegree(m) eq 3};
  if BaseRing(Parent(cubic)) cmpeq Rationals() then
    vprintf ThreeDescent, 2: "  Making coefficients integral...\n";
    cubic *:= LCM([Denominator(c) : c in Coefficients(cubic)]);
  end if;
  return cubic;
end intrinsic;