// 3coverings.m
//
// Michael Stoll, started 2002-05-15
//
// Some glueing code in order to produce a list of ELS 3-coverings for
// a given elliptic curve.

declare verbose ThreeCoverings, 3;

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;

function cbrt(z)
  pol := MinimalPolynomial(z);
  K := NumberField(pol);
  assert Evaluate(pol, K.1) eq 0;
  /* -- this too often is wrong and then takes a long time
  flag, rt := IsPower(K.1, 3);
  if flag then
    cbrtz := Evaluate(Parent(pol)!Eltseq(rt), z);
    assert cbrtz^3 eq z;
    return cbrtz;
  end if;
  printf "IsPower(_, 3) returned wrong result!\n";
  // next attempt
  */
  flag, rt := HasRoot(PolynomialRing(K)![-K.1, 0, 0, 1]);
  if flag then
    cbrtz := Evaluate(Parent(pol)!Eltseq(rt), z);
    assert cbrtz^3 eq z;
    return cbrtz;
  end if;
  printf "HasRoot returned wrong result!\n";
  fact := Factorisation(Evaluate(pol, Parent(pol).1^3));
  assert exists(f){f : f in fact | Degree(f[1]) eq Degree(pol)};
  flag, rt := HasRoot(f[1], K);
  if flag then
    cbrtz := Evaluate(Parent(pol)!Eltseq(rt), z);
    assert cbrtz^3 eq z;
    return cbrtz;
  end if;
  error "HasRoot returned wrong result again!\n";
end function;

// Transform elements as returned by ThreeSelmerGroup into r's as needed
// by other programs
function AlphaToR(A, maps, isoKA, alpha)
  // A: vector space representing algebra A
  // maps: tuple of maps, as given by SetupAlgebras
  // isoKA: isomorphism of Parent(alpha) with A1
  // alpha: Selmer group element as given by ThreeSelmerGroup
  a1 := isoKA(alpha);
  A1 := Parent(a1);
  neg := maps[11]; // hom<A1 -> A1 | A1.1, -A1.2>;
  na1 := a1*neg(a1);
  cna1 := cbrt(na1);
  a := maps[6](<1, a1>); // in A
  d6ic := maps[9](maps[2](maps[6](<1, invert(cna1)>)))[6];
  d6 := maps[9](maps[2](maps[6](<1, na1>)))[6];
  ai := maps[6](<1, invert(a1)>);
  r3 := maps[5](maps[1](a, a), maps[2](ai)); // (a tensor a)/Delta(a)
  rr3 := maps[9](r3); // rr3 = <1, 1, 1, rr3[4], rr3[5], rr3[6]>
  assert rr3[1] eq 1 and rr3[2] eq 1 and rr3[3] eq 1;
  assert rr3[4]*na1 eq a1^3 and rr3[5] eq na1;
  rr := <1, 1, 1, a1*invert(cna1), cna1, d6ic*cbrt(d6*rr3[6])>;
  return maps[8](rr);
end function;

function GetIsoKA(K, A)
  assert Degree(K) eq Dimension(A);
  pol := MinimalPolynomial(A.2);
  n := Degree(K);
  assert Degree(pol) eq n;
  rt := Roots(pol, K)[1,1];
  V, m := VectorSpace(A);
  h := hom<V -> V | [V | Eltseq(rt^i) : i in [0..n-1]]>;
  genA := Evaluate(Parent(pol)!Eltseq(V.2 @@ h), A.2);
  return hom<K -> A | genA>;
end function;

intrinsic ThreeCoverings(E::CrvEll : Bound := 0) -> SeqEnum
{Computes a list of all everywhere locally soluble nontrivial 3-coverings 
 for the elliptic curve E over Q, as reduced and minimal plane cubics.}
  // compute Selmer group fisrt
  vprintf ThreeCoverings, 1: "ThreeCoverings for E =\n%o\n", E;
  vprintf ThreeCoverings, 2: "  Computing 3-Selmer group ...\n";
  rb, r, V, mAS3 := ThreeSelmerGroup(E : Bound := Bound);
  if r ne 0 then
    K := Parent(mAS3(V!0));
    // set up algebras
    vprintf ThreeCoverings, 2: "  Setting up algebras and stuff ...\n";
    E1 := E; // WeierstrassModel(E);
    A, AA, maps, coords := SetupAlgebras(E1);
    A1 := Parent((maps[7](A!0))[2]);
    vprintf ThreeCoverings, 2: "  Getting isomorphism K --> A1 ...\n";
    isoKA := GetIsoKA(K, A1);
    elts := {v : v in V | v ne 0};
    elements := [V | ];
    while not IsEmpty(elts) do
      e := Rep(elts);
      elts diff:= {e, -e};
      Append(~elements, e);
    end while;
    vprintf ThreeCoverings, 2:
            "  There are %o different covering spaces.\n", #elements;
    vprintf ThreeCoverings, 2: "  Computing r_alpha's ...\n";
    rlist := [AlphaToR(A, maps, isoKA, mAS3(alpha)) : alpha in elements];
    vprintf ThreeCoverings, 2:
            "  Determining bad primes for maximal order ...\n";
    badprimes := {2, 3} join Seqset(BadPrimes(E));
    for v in Basis(V) do
      nv := Norm(mAS3(v));
      flag, nrt := IsPower(nv, 3);
      assert flag;
      badprimes join:= Seqset(PrimeDivisors(Numerator(nrt))
                               cat PrimeDivisors(Denominator(nrt)));
    end for;
    badprimes := Setseq(badprimes);
    vprintf ThreeCoverings, 2: "  Bad primes: %o\n", badprimes;
    P3<x,y,z> := PolynomialRing(Rationals(), 3);
    result := [P3 | ];
    for i := 1 to #rlist do
      r := rlist[i];
      vprintf ThreeCoverings, 2:
              "    Computing structure constants for A_alpha ...\n";
      Aa, Aa1 := GetCSAlgebra(A, maps, coords, r);
      vprintf ThreeCoverings, 2: "    Computing maximal order ...\n";
      bas := MaximalOrder(Aa, badprimes);
      vprintf ThreeCoverings, 2: "    Trivializing A_alpha ...\n";
      iso := Trivialize(Aa, bas);
      vprintf ThreeCoverings, 2: "    Computing equations in P^7 ...\n";
      eqns1 := CoveringCurve1(A, maps, coords, r);
      vprintf ThreeCoverings, 2: "    Computing plane cubic ...\n";
      C := CoveringCurve2(eqns1, iso);
      vprintf ThreeCoverings, 2: "    Minimising cubic ...\n";
      C := Minimise(C);
      vprintf ThreeCoverings, 2: "    Reducing cubic ...\n";
      C := Reduce(C);
      vprintf ThreeCoverings, 2: "    Cubic is\n    %o\n", C;
      Append(~result, C);
    end for;
    return [Sign(Coefficients(f)[1])*f : f in result];
  else
    return [];
  end if;
end intrinsic;