/* A WuthrichType example is a semisimple elliptic curve over the rationals and a cyclic extension K/Q
   of odd prime degree l such that E(K)_l is a permutation module of the form \Zl^{n_0} \oplus \ZlG.
   We also require all the assumptions used in the paper etncIV (essentially checked in ConditionB).

   The intrinsic DoWuthrichTypeExample produces all the data and writes it on the fiel filename that is
   needed to check the ETNC at l over the maximal order (using the results of etncIV).

   These computations were not included in the paper since by Theorem 2.8 the ETNC at l over the
   maximal order is equivalent to BSD at l for all the subfields. In the paper we wanted to concentrate on
   examples which go beyond BSD, however, these computations are  still of some interest. At least they provide
   evidence for BSD.

*/

intrinsic DoWuthrichTypeExample(filename, data) -> BoolElt
{}
    local i;

    E := data[1]; p := data[2]; l := data[3]; n := data[4]; prec := data[5]; Bound := data[6]; CheckSaturation := data[7];
    fprintf filename, "****************************************************************************************\n";

    SetDefaultRealFieldPrecision(prec);

    C<zeta> := CyclotomicField(p);
    w := PrimitiveRoot(p);
    w := w^(l^n) mod p;
    delta1 := &+[zeta^(w^i) : i in [0..((p-1)/l^n)-1]];
    K := sub<C | delta1>;

    /* Compute Artin representations etc. */
    print "Compute the Artin representations";
    A := ArtinRepresentations(K);
    G := K`artinrepdata`G;
    h := CompatibleIso(K, G);
    OK := MaximalOrder(K);
    QG := InitGroupAlgebra(G);
    delta := [h(g)(delta1) : g in G | Abs( Conjugates(h(g)(K!delta1))[1] - Conjugates(delta1)[1]  ) le 10^-10 ][1];

    print "Condition B is ", ConditionB(E, OK, l);
    fprintf filename, "Condition B is %o\n", ConditionB(E, OK, l);

    EK := BaseChange(E, K);
    MW, iota := MordellWeilGroup(E);

    fprintf filename, "E = %o\n", E;
    fprintf filename, "K defined by %o\n", MinimalPolynomial(K.1);
    fprintf filename, "****************************************************************************************\n";

    print "E(Q) = ", MW, " ---- ", [iota(g) : g in Generators(MW)];
    fprintf filename, "E(Q) = %o\n ------- %o\n", MW, [iota(g) : g in Generators(MW)];
    Qpts := [ iota(g) : g in Generators(MW) | Order(g) eq 0 ];
    pts := [ EK ! [K!x : x in ElementToSequence(iota(g))] : g in Generators(MW) | Order(g) eq 0 ];
    reg := MyRegulator(pts);
    r := #pts;
    print "The rank of E(Q) is ", r;

    print "Compute resolvents, periods and L-values";
    /* Compute the resolvents */
    Res := Resolvents(QG, h, delta);
    fprintf filename, "Resolvents =\n%o\n", Res;

    /* Compute the real and purely imaginary period of E/Q */
    Ovals := CompOvals(QG, E, OK, h);
    fprintf filename, "Periods =\n%o\n", Ovals;


    /* Compute the leading terms of the L-series and the order of vanishing */
    SetVerbose("LSeries", 1);
    Lvec := InitLSeries(QG, E, K, prec);
    Lvals,orders := Evaluate(Lvec);
    E`Lvals := Lvals; E`orders := orders; E`OK := OK; E`h := h;
    print "The analytic rank conjecturally is ", orders;
    fprintf filename, "The analytic rank conjecturally is %o\n", orders;
    fprintf filename, "Leading terms =\n%o\n", Lvals;

    print "Search for a point";
    Kpts := MyPoints(E, EK, 10);   /* Bound=10 is very ad hoc; should be changed; Bound should be set depending on the example */
    Kpts := [P : P in Kpts | not P[1] in Rationals()];
    if #Kpts eq 0 then
        fprintf filename, "No point found\n";
        print "No point found";
        return false;
    end if;
    /* We simple take the first point in Kpts; probably one should use a point of minimal height to have better performance */
    R := Kpts[1];
    Q1 := E ! Trace(R, h);
    log := NaiveLogModTorsion(Qpts,Q1, Bound);  /* The bound has to vary according to the example */
    print "log = ", log;
    fprintf filename, "R = %o\n", R;
    fprintf filename, "Tr(R) = %o\nlog = %o\n", Q1, log;

    /* We search an index i0 such that log[i0] is minimal and ne 0; 
       we should also require, if possible, that log[i0] is prime to l; this is not yet done, so if 
       log[i0] is not prime to l then this must be visible in the output and disturb the validity of the ETNC
    */
    d, i0 := Max([Abs(l) : l in log]);
    for i:=1 to r do
        if Abs(log[i]) lt d and log[i] ne 0 then
            i0 := i; d := Abs(log[i]);
        end if;
    end for;
    Qs := [EK ! Q1] cat [EK ! Qpts[i] : i in [1..r] | i ne i0];

    reg2 := MyRegulator(Qs);
    print "d = ", d, " reg2 / reg = ", reg2/reg;  /* Simple check: reg2 / reg must equal d^2 */

    sigma := G.1;
    gens := [Qs[i] : i in [2..r]] cat [GAction(g, R, h) : g in G];
    if CheckSaturation then
        fprintf filename, "V is %o-saturated : %o\n", l, IslSaturated(gens, l);
    else
        fprintf filename, "l-Saturation of V not checked\n";
    end if;
    fprintf filename, "V = %o\n", gens;
    Kreg := MyRegulator(gens);

    RR := [ [Kreg] ] cat [ [1 : chi in QG`X[j]] : j in [2..#QG`X] ];   /* Ad hoc, only needed to compute the order of Sha */
    /* Note that (conjecturally by BSD) ts is the order of Sha times the square of the index of gens in E(K) */
    ts := AnalyticTateShafarevicGroup(E, QG, Lvals, RR, OK);
    fprintf filename, "#Analytic order of Sha times [E(K) : <gens>]^2 = %o\n", ts;


    /* We compute the regulators as in etncIV for the non trivial characters */
    Rvals := lCyclicRegulator(QG, h, R , prec);
    /* Now consider the trivial character */
    
    RegMat3 := Matrix(RealField(), r-1,1, [HeightPairing(Qs[1], Qs[k], K) :  k in [2..r]]);
    RegMat2 := Matrix(RealField(), r-1,1, [HeightPairing(Qs[k], R, K) / Degree(K) :  k in [2..r]]);
    RegMat1 := Matrix(RealField(), r-1,r-1, [HeightPairing(Qs[i], Qs[k], K) / Degree(K) : i in [2..r], k in [2..r]]);

    RegMat := ZeroMatrix(RealField(), r, r);
    for i:=1 to r-1 do
        for j:=1 to r-1 do
            RegMat[i,j] := RegMat1[i,j];
        end for;
    end for;
    for j:=1 to r-1 do RegMat[j,r] := RegMat2[j,1]; end for;
    for j:=1 to r-1 do RegMat[r, j] := RegMat3[j, 1]; end for;
    RegMat[r,r] := HeightPairing(R, Trace(R, h),K);

    Rvals[1,1] := Determinant(RegMat);

    C<I> := ComplexField();
    pi := Pi(C);
    zeta := Exp( 2*pi*I / l);

    /* Simple check: the regulator must equal the product of the equivariant regulators */
    print "Difference of regulators: ", Kreg - Abs( Rvals[1,1] * Rvals[2,1]^(l-1) * (zeta - 1)^((r-1)*(l-1)) );
    fprintf filename, "Difference of regulators: %o\n", Kreg - Abs( Rvals[1,1] * Rvals[2,1]^(l-1) * (zeta - 1)^((r-1)*(l-1)) );

    Z := ComputeZ(QG, E, Lvals, Res, Ovals, Rvals);
    print "Z = ", Z;
    fprintf filename, "Z = %o\n", Z;
    RatZ := MakeRational(Z, QG);
    print "RatZ = ", RatZ;
    fprintf filename, "RatZ = %o\n", RatZ;

    eul := DualLocalEulerFactors(QG, E, E`OK, E`h, p);
    u := [* [ RatZ[i,j] * eul[i,j] : j in [1..#RatZ[i]] ] : i in [1..#RatZ] *];

    /* The u here is NOT the element u in the proof of Corollary 2.11 because it is not adapted by the elements
       delta_psi. The element calL is therefore just the sum over all non-trivial components of RatZ */

    t0 := 1;
    a := ShaAnnihilator(QG, t0, u);
    print "Annihilator = ", a;
    g0 := G.1;
    fprintf filename, "Ann = ";
    for i:=0 to #G-1 do
        print Coefficient(a, g0^i), "*S^", i, " + ";
        if i eq #G-1 then
             fprintf filename, "%o*S^%o\n", Coefficient(a, g0^i), i;
        else
            fprintf filename, "%o*S^%o + ", Coefficient(a, g0^i), i;
        end if;
    end for;

    /* We check (i) of Theorem 2.12. This is not included in the paper because for degree l extensions the result is
       implied from the validity of the ETNC over the maximal order which, in turn, is implied from BSD for all the subfields.
    */
    h := r-1;
    sigma := QG ! G.1;
    ZG := ModuleInit( IdentityMatrix(Rationals(), #G) );
    I := ModuleScalarProd(QG, ZG, (sigma - 1)^h);
    J := ModuleScalarProd(QG, ZG, (sigma - 1)^(h+1));
    a_in_I := IsModuleElt(QG, a, I);
    a_in_J := IsModuleElt(QG, a, J);
    fprintf filename, "Ann in I_{G,p}^%o ---------> %o\n", h, a_in_I;
    fprintf filename, "Ann in I_{G,p}^%o ---------> %o\n", h+1, a_in_J;
    printf "Ann in I_{G,p}^%o ---------> %o\n", h, a_in_I;
    printf "Ann in I_{G,p}^%o ---------> %o\n", h+1, a_in_J;

    /* We check (iii) of Theorem 2.8. Not included in the paper for the same reason as above. */
    OC := MaximalOrder( Parent(RatZ[2,1]) );
    P := Factorization(l*OC)[1,1];
    print "val(RatZ[1,1]) = ", Valuation(Rationals()!RatZ[1,1], l), "   expected ", 0;
    print "val(RatZ[2,1]) = ", Valuation(RatZ[2,1], P), "   expected ", r-1;
    fprintf filename, "val(RatZ[1,1]) = %o  expected %o\n", Valuation(Rationals()!RatZ[1,1], l), 0;
    fprintf filename, "val(RatZ[2,1]) = %o  expected %o\n", Valuation(RatZ[2,1], P), r-1;

    fprintf filename, "********************** Example done ****************************************************\n\n\n";
    fprintf filename, "****************************************************************************************\n";


    return r-1 eq Valuation(RatZ[2,1], P);

end intrinsic;


/* A lSquared example is a semisimple elliptic curve over the rationals and a cyclic extension K/Q
   of odd prime power degree l^2 such that E(K)_l is a permutation module of the form \Zl^{n_0} \oplus \Zl[G/H] where
   H is the subgroup of order l. In all examples so far we had l = 3.
   We also require all the assumptions used in the paper etncIV (essentially checked in ConditionB).

   The intrinsic DolSquaredExample produces all the data and writes it on the file filename that is
   needed to check the ETNC at l over the maximal order (using the results of etncIV).

   We did not include examples with n_0 = 0 in the paper because here calL is essentially the Mazur/Tate element
   and BSD implies the assertions of part (i) of Cor. 2.11 and Th. 2.12.
*/

intrinsic DolSquaredExample(filename, data) -> BoolElt
{}
    local i;

    E := data[1]; p := data[2]; l := data[3]; n := data[4]; prec := data[5]; Bound := data[6]; CheckSaturation := data[7];
    fprintf filename, "****************************************************************************************\n";

    SetDefaultRealFieldPrecision(prec);

    C<zeta> := CyclotomicField(p);
    w := PrimitiveRoot(p);
    w := w^(l^n) mod p;
    delta1 := &+[zeta^(w^i) : i in [0..((p-1)/l^n)-1]];
    K := sub<C | delta1>;

    /* Compute Artin representations etc. */
    print "Compute the Artin representations";
    A := ArtinRepresentations(K);
    G := K`artinrepdata`G;
    h := CompatibleIso(K, G);
    OK := MaximalOrder(K);
    QG := InitGroupAlgebra(G);
    delta := [h(g)(delta1) : g in G | Abs( Conjugates(h(g)(K!delta1))[1] - Conjugates(delta1)[1]  ) le 10^-10 ][1];

    fprintf filename, "E = %o\n", E;
    fprintf filename, "K defined by %o\n", MinimalPolynomial(K.1);
    fprintf filename, "****************************************************************************************\n";

    print "Condition B is ", ConditionB(E, OK, l);
    fprintf filename, "Condition B is %o\n", ConditionB(E, OK, l);

    EK := BaseChange(E, K);
    MW, iota := MordellWeilGroup(E);

    print "E(Q) = ", MW, " ---- ", [iota(g) : g in Generators(MW)];
    fprintf filename, "E(Q) = %o\n ------- %o\n", MW, [iota(g) : g in Generators(MW)];

    Qpts := [ iota(g) : g in Generators(MW) | Order(g) eq 0 ];
    pts := [ EK ! [K!x : x in ElementToSequence(iota(g))] : g in Generators(MW) | Order(g) eq 0 ];
    reg := MyRegulator(pts);
    r := #pts;
    print "The rank of E(Q) is ", r;

    print "Compute resolvents, periods and L-values";
    /* Compute the resolvents */
    Res := Resolvents(QG, h, delta);
    fprintf filename, "Resolvents =\n%o\n", Res;

    /* Compute the real and purely imaginary period of E/Q */
    Ovals := CompOvals(QG, E, OK, h);
    fprintf filename, "Periods =\n%o\n", Ovals;


    /* Compute the leading terms of the L-series and the order of vanishing */
    SetVerbose("LSeries", 1);
    Lvec := InitLSeries(QG, E, K, prec);
    Lvals,orders := Evaluate(Lvec);
    E`Lvals := Lvals; E`orders := orders; E`OK := OK; E`h := h;
    print "The analytic rank conjecturally is ", orders;
    fprintf filename, "The analytic rank conjecturally is %o\n", orders;
    fprintf filename, "Leading terms =\n%o\n", Lvals;

    print "Search for a point";
    /* We search the point R in the subfield F with [F:Q] = l */
    S := Subfields(K,l);
    F := S[1,1];
    EF := BaseChange(E, F);
    Kpts := MyPoints(E, EF, 50);
    Kpts := [EK ! P : P in Kpts | not P[1] in Rationals()];
    if #Kpts eq 0 then
        print "No point found";
        fprintf filename, "No point found\n";
        return false;
    end if;
    /* Among the points found in E(F) we look for one with minimal height. */
    i0 := 1; hmin := HeightPairing(Kpts[1], Kpts[1]);
    for i:=2 to #Kpts do
        height := HeightPairing(Kpts[i], Kpts[i]);
	if height lt hmin then
	    hmin := height;
	    i0 := i;
	end if;
    end for;
    R := Kpts[i0];
    sigma := G.1;
    Q1 := E ! &+[ GAction(sigma^i, R, h) : i in [0..l-1] ]; 
    log := NaiveLogModTorsion(Qpts,Q1, Bound);
    print "log = ", log;
    fprintf filename, "R = %o\n", R;
    fprintf filename, "Tr(R) = %o\nlog = %o\n", Q1, log;

    d, i0 := Max([Abs(l) : l in log]);
    for i:=1 to r do
        if Abs(log[i]) lt d and log[i] ne 0 then
            i0 := i; d := Abs(log[i]);
        end if;
    end for;
    Qs := [EK ! Q1] cat [EK ! Qpts[i] : i in [1..r] | i ne i0];

    reg2 := MyRegulator(Qs);
    print "d = ", d, " reg2 / reg = ", reg2/reg; /* A simple check: reg2/reg must be d^2 */
    fprintf filename, "d = %o  reg2/reg = %o (we should have d^2 = reg2/reg)\n", d, reg2/reg; 

    H := [S`subgroup : S in Subgroups(G) | S`order eq l][1];
    GmodH, pi := quo<G | H>;
    GmodH := [(pi^-1)(g) : g in GmodH];
    gens := [Qs[i] : i in [2..r]] cat [GAction(g, R, h) : g in GmodH];
    if CheckSaturation then
        bool := IslSaturated(gens, l);
        print "Is saturated for l = ", l, "..............", bool;
        fprintf filename, "V is %o-saturated : %o\n", l, bool;
    else
        fprintf filename, "l-Saturation of V not checked\n";
    end if;
    fprintf filename, "V = %o\n", gens;

    Kreg := MyRegulator(gens);

    RR := [ [Kreg] ] cat [ [1 : chi in QG`X[j]] : j in [2..#QG`X] ];
    ts := AnalyticTateShafarevicGroup(E, QG, Lvals, RR, OK);
    print "Analytic order of Sha times [E(K) : <gens>]^2 = ", ts;
    fprintf filename, "Analytic order of Sha times [E(K) : <gens>]^2 = %o\n", ts;

    sig := [orders[1,1]-1, orders[2,1], orders[3,1]];   /* sig = [m_0, m_1, m_2] in the notation of etncIV. 
                                                           Attention: Works only for m_2 = 0 (which is the only case where we use it so far
							*/
    print "sig = ", sig;
    fprintf filename, "(m_0, m_1, m_2) = %o\n", sig;
    if sig[1] eq 0 and sig[2] eq 1 and sig[3] eq 0 then
        P := [R];
    end if;
    if sig[1] eq 1 and sig[2] eq 1 and sig[3] eq 0 then
        P := [gens[1], R];
    end if;
    Rvals := lSquaredCyclicRegulator(QG, h, P , sig, prec);

    Z := ComputeZ(QG, E, Lvals, Res, Ovals, Rvals);
    print "Z = ", Z;
    fprintf filename, "Z = %o\n", Z;
    RatZ := MakeRational(Z, QG);
    print "RatZ = ", RatZ;
    fprintf filename, "RatZ = %o\n", RatZ;

    eul := DualLocalEulerFactors(QG, E, E`OK, E`h, p);
    u := [* [ RatZ[i,j] * eul[i,j] : j in [1..#RatZ[i]] ] : i in [1..#RatZ] *];
    /* The u here is NOT the element u in the proof of Corollary 2.11 because it is not adapted by the elements
       delta_psi. The element calL in this therefore just the sum over all non-trivial components of RatZ */

    t0 := 2;
    a := ShaAnnihilator(QG, t0, u);
    print "Annihilator = ", a;
    g0 := G.1;
    fprintf filename, "Ann = ";
    for i:=0 to #G-1 do
        print Coefficient(a, g0^i), "*S^", i, " + ";
        if i eq #G-1 then
             fprintf filename, "%o*S^%o\n", Coefficient(a, g0^i), i;
        else
            fprintf filename, "%o*S^%o + ", Coefficient(a, g0^i), i;
        end if;
    end for;

    /* We now check part (i) of Cor. 2.11 or Th. 2.12 */
    h := &+sig;
    sigma := QG ! G.1;
    ZG := ModuleInit( IdentityMatrix(Rationals(), #G) );
    I := ModuleScalarProd(QG, ZG, (sigma - 1)^h);
    J := ModuleScalarProd(QG, ZG, (sigma - 1)^(h+1));
    a_in_I := IsModuleElt(QG, a, I);
    a_in_J := IsModuleElt(QG, a, J);
    fprintf filename, "Ann in I_{G,p}^%o ---------> %o\n", h, a_in_I;
    fprintf filename, "Ann in I_{G,p}^%o ---------> %o\n", h+1, a_in_J;
    printf "Ann in I_{G,p}^%o ---------> %o\n", h, a_in_I;
    printf "Ann in I_{G,p}^%o ---------> %o\n", h+1, a_in_J;


    /* We also check part (iii) of Theorem 2.8, however, this is not included in the paper since it is implied by BSD for the subfields. */
    OC2 := MaximalOrder( Parent(RatZ[3,1]) );
    P2 := Factorization(l*OC2)[1,1];

    OC1 := MaximalOrder( Parent(RatZ[2,1]) );
    P1 := Factorization(l*OC1)[1,1];

    if sig[1] eq 0 and sig[2] eq 1 and sig[3] eq 0 then
        print "val(RatZ[1,1]) = ", Valuation(Rationals()!RatZ[1,1], l), "  expected  ", 0;
        print "val(RatZ[2,1]) = ", Valuation(RatZ[2,1], P1), "  expected  ", 0;
        print "val(RatZ[3,1]) = ", Valuation(RatZ[3,1], P2), "  expected  ", l;
        fprintf filename, "val(RatZ[1,1]) = %o expected %o\n", Valuation(Rationals()!RatZ[1,1], l), 0;
        fprintf filename, "val(RatZ[2,1]) = %o expected %o\n", Valuation(RatZ[2,1], P1), 0;
        fprintf filename, "val(RatZ[3,1]) = %o expected %o\n", Valuation(RatZ[3,1], P2), l;
        return l eq Valuation(RatZ[3,1], P2) and 0 eq Valuation(RatZ[2,1], P1) and 0 eq Valuation(Rationals()!RatZ[1,1], l);
    end if;
    if sig[1] eq 1 and sig[2] eq 1 and sig[3] eq 0 then
        print "val(RatZ[1,1]) = ", Valuation(Rationals()!RatZ[1,1], l), "  expected  ", 0;
        print "val(RatZ[2,1]) = ", Valuation(RatZ[2,1], P1), "  expected  ", 1;
        print "val(RatZ[3,1]) = ", Valuation(RatZ[3,1], P2), "  expected  ", l+1;
        fprintf filename, "val(RatZ[1,1]) = %o expected %o\n", Valuation(Rationals()!RatZ[1,1], l), 0;
        fprintf filename, "val(RatZ[2,1]) = %o expected %o\n", Valuation(RatZ[2,1], P1), 1;
        fprintf filename, "val(RatZ[3,1]) = %o expected %o\n", Valuation(RatZ[3,1], P2), l+1;
        return l+1 eq Valuation(RatZ[3,1], P2) and 1 eq Valuation(RatZ[2,1], P1) and 0 eq Valuation(Rationals()!RatZ[1,1], l);
    end if;
end intrinsic;


/* A Perfect example is a semisimple elliptic curve over the rationals and a cyclic extension K/Q
   of odd prime degree l such that E(K)_l is a permutation module of the form \ZlG].
   We also require all the assumptions used in the paper etncIV (essentially checked in ConditionB).

   The intrinsic DoperfectExample produces all the data and writes it on the file filename that is
   needed to check the ETNC at l over the maximal order (using the results of etncIV).

   These are the examples described in Sec.5.1 of the paper.
*/

intrinsic DoPerfectExample(filename, data) -> BoolElt, BoolElt
{}
    local i;

    E := data[1]; p := data[2]; l := data[3]; n := data[4]; prec := data[5]; Bound := data[6]; CheckSaturation := data[7];
    fprintf filename, "****************************************************************************************\n";

    SetDefaultRealFieldPrecision(prec);

    C<zeta> := CyclotomicField(p);
    w := PrimitiveRoot(p);
    w := w^(l^n) mod p;
    delta1 := &+[zeta^(w^i) : i in [0..((p-1)/l^n)-1]];
    K := sub<C | delta1>;

    /* Compute Artin representations etc. */
    print "Compute the Artin representations";
    A := ArtinRepresentations(K);
    G := K`artinrepdata`G;
    h := CompatibleIso(K, G);
    OK := MaximalOrder(K);
    QG := InitGroupAlgebra(G);
    delta := [h(g)(delta1) : g in G | Abs( Conjugates(h(g)(K!delta1))[1] - Conjugates(delta1)[1]  ) le 10^-10 ][1];

    print "Condition B is ", ConditionB(E, OK, l);
    fprintf filename, "Condition B is %o\n", ConditionB(E, OK, l);

    EK := BaseChange(E, K);
    MW, iota := MordellWeilGroup(E);

    fprintf filename, "%o = %o\np =%o, l = %o n = %o\n", CremonaReference(E), E, p, l, n;
    fprintf filename, "K defined by %o\n", MinimalPolynomial(K.1);
    fprintf filename, "data = %o\n\n", data;

    print "E(Q) = ", MW, " ---- ", [iota(g) : g in Generators(MW)];
    fprintf filename, "E(Q) = %o\n ------- %o\n", MW, [iota(g) : g in Generators(MW)];
    Qpts := [ iota(g) : g in Generators(MW) | Order(g) eq 0 ];
    pts := [ EK ! [K!x : x in ElementToSequence(iota(g))] : g in Generators(MW) | Order(g) eq 0 ];
    reg := MyRegulator(pts);
    r := #pts;
    print "The rank of E(Q) is ", r;

    print "Compute resolvents, periods and L-values";
    /* Compute the resolvents */
    Res := Resolvents(QG, h, delta);
    fprintf filename, "Resolvents =\n%o\n", Res;

    /* Compute the real and purely imaginary period of E/Q */
    Ovals := CompOvals(QG, E, OK, h);
    fprintf filename, "Periods =\n%o\n", Ovals;


    /* Compute the leading terms of the L-series and the order of vanishing */
    SetVerbose("LSeries", 1);
    Lvec := InitLSeries(QG, E, K, prec);
    Lvals,orders := Evaluate(Lvec);
    E`Lvals := Lvals; E`orders := orders; E`OK := OK; E`h := h;
    print "The analytic rank conjecturally is ", orders;
    fprintf filename, "The analytic rank conjecturally is %o\n", orders;
    fprintf filename, "Leading terms =\n%o\n", Lvals;

    print "Search for a point";
    Kpts := MyPoints(E, EK, Bound);
    Kpts := [P : P in Kpts | not P[1] in Rationals()];
    if #Kpts eq 0 then
        fprintf filename, "No point found\n";
        print "No point found";
        return false, false;
    end if;
    found := false; i := 1;
    while not found and i lt #Kpts do
        R := Kpts[i];
        Q1 := E ! Trace(R, h);
        log := NaiveLogModTorsion(Qpts,Q1, 50);  /* The bound 50 has to vary according to the example */
        print "log = ", log;
        if not IsDivisibleBy(log[1], l) then found := true; end if;
        i := i+1;
    end while;
    if not found then
        fprintf filename, "No point with prime to l log found. Further results not correct\n";
        print "No point with prime to l log found. Further results not correct\n";
        R := Kpts[1];
    end if;

   fprintf filename, "R = %o\n", R;
   fprintf filename, "Tr(R) = %o\nlog = %o\n", Q1, log;

    d, i0 := Max([Abs(l) : l in log]);
    for i:=1 to r do
        if Abs(log[i]) lt d and log[i] ne 0 then
            i0 := i; d := Abs(log[i]);
        end if;
    end for;
    Qs := [EK ! Q1] cat [EK ! Qpts[i] : i in [1..r] | i ne i0];


    reg2 := MyRegulator(Qs);
    print "d = ", d, " reg2 / reg = ", reg2/reg;
    fprintf filename, "d = %o  reg2/reg = %o (we should have d^2 = reg2/reg)\n", d, reg2/reg;


    sigma := G.1;
    gens := [Qs[i] : i in [2..r]] cat [GAction(g, R, h) : g in G];
    if CheckSaturation then
        fprintf filename, "V is %o-saturated : %o\n", l, IslSaturated(gens, l);
    else
        fprintf filename, "l-Saturation of V not checked\n";
    end if;
    fprintf filename, "V = %o\n", gens;
    Kreg := MyRegulator(gens);

    RR := [ [Kreg] ] cat [ [1 : chi in QG`X[j]] : j in [2..#QG`X] ];
    ts := AnalyticTateShafarevicGroup(E, QG, Lvals, RR, OK);
    fprintf filename, "Analytic order of Sha times [E(K) : <gens>]^2 = %o\n", ts;


    Rvals := lCyclicRegulator(QG, h, R , prec);

    Z := ComputeZ(QG, E, Lvals, Res, Ovals, Rvals);
    print "Z = ", Z;
    fprintf filename, "Z = %o\n", Z;
    RatZ := MakeRational(Z, QG);
    print "RatZ = ", RatZ;
    fprintf filename, "RatZ = %o\n", RatZ;

    /* We check the full ETNC, hence we must also consider the Euler factor at p (which is the conductor of K). */
    eul := DualLocalEulerFactors(QG, E, E`OK, E`h, p);
    print "eul = ", eul;
    fprintf filename, "eul = %o\n", eul;

    Lstar := [ [ Z[i,j]*Conjugates(eul[i,j])[1] : j in [1..#Z[i]] ] : i in [1..#Z] ];
    print "Lstar = ", Lstar;
    fprintf filename, "Lstar = %o\n", Lstar;

    RatLstar := MakeRational(Lstar, QG);
    print "RatLstar = ", RatLstar;
    fprintf filename, "RatLstar = %o\n", RatLstar;

    u := [* RatLstar[i,1] : i in [1..#RatZ] *];
    print "u = ", u;
    fprintf filename, "u = %o\n", u;

    /* We have to check that (a) u[1] and u[2] are l-units and (b) u[1] equiv u[2] (mod (1-zeta)) */ 

    C := Parent(RatZ[2,1]);
    OC := MaximalOrder( C );
    P := Factorization(l*OC)[1,1];
    print "val(u[1]) = ", Valuation(Rationals()!u[1], l), "   expected ", 0;
    print "val(u[2,1]) = ", Valuation(C ! u[2], P), "   expected ", 0;
    fprintf filename, "val(u[1]) = %o  expected %o\n", Valuation(Rationals()!u[1], l), 0;
    fprintf filename, "val(u[2]) = %o  expected %o\n", Valuation(C ! u[2], P), 0;

    print "val(u[1] - u[2]) = ", Valuation(C ! (u[1] - u[2]), P), "   expected >= ", 1; 
    fprintf filename, "val(u[1] - u[2]) = %o expected >= 1\n", Valuation(C ! (u[1] - u[2]), P);

    return true, (Valuation(C ! (u[1] - u[2]), P) gt 0) and (Valuation(Rationals()!u[1], l) eq 0) and (Valuation(C ! u[2], P) eq 0);

end intrinsic;


/* A NonPerfectAnnihilator example is a an example of type E(K)_l \simeq \Zl^{n_0}. Here we check whether part (i) of
   Cor.2.11 or Th. 2.12 holds.
   As explained in Sec.5.2 we are particularly interested in cases where [K:Q] = l^2 and n_0 ge 2.

   We produce all the date, and write it to the file fielname, which is needed to verify part (i) of Cor.2.11 or Th. 2.12 holds.
*/

intrinsic DoNonPerfectAnnihilator(filename, data) -> BoolElt, BoolElt
{}
   local i;

    E := data[1]; p := data[2]; l := data[3]; n := data[4]; prec := data[5]; Bound := data[6]; CheckSaturation := data[7];
    fprintf filename, "****************************************************************************************\n";
    fprintf filename, "data = %o\n", data;
    fprintf filename, "E = %o\n\n", CremonaReference(E);
    print "data = ", data, "E = ", CremonaReference(E);

    SetDefaultRealFieldPrecision(prec);

    if not n le 2 then
        fprintf filename, "n = %o not allowed\n", n;
        print "n = ", n, " not allowed";
        return false, false;
    end if;

    if n eq 2 then
        /* Ad hoc; could and should be done in a more efficient way. */
        fprintf filename, "--------------->>>>>  We first do the n = 1 case in order to check the order of Sha in the intermediate field\n";
        print "--------------->>>>> We first do the n = 1 case in order to check the order of Sha in the intermediate field";
        inter_data := data; inter_data[4] := 1;

        computable, ETNCvalid := DoNonPerfectAnnihilator(filename, inter_data); 
        if not computable or not ETNCvalid then 
            return computable, ETNCvalid; 
        end if;

        fprintf filename, "--------------->>>>>  Now we do  the n = 2 case\n";
        print "--------------->>>>> Now we do the n = 2 case";
    end if;

    C<zeta> := CyclotomicField(p);
    w := PrimitiveRoot(p);
    w := w^(l^n) mod p;
    delta1 := &+[zeta^(w^i) : i in [0..((p-1)/l^n)-1]];
    K := sub<C | delta1>;

    /* Compute Artin representations etc. */
    print "Compute the Artin representations";
    A := ArtinRepresentations(K);
    G := K`artinrepdata`G;
    h := CompatibleIso(K, G);
    OK := MaximalOrder(K);
    QG := InitGroupAlgebra(G);
    delta := [h(g)(delta1) : g in G | Abs( Conjugates(h(g)(K!delta1))[1] - Conjugates(delta1)[1]  ) le 10^-10 ][1];

    b := ConditionB(E, OK, l);
    print "Condition B is ", b;
    fprintf filename, "Condition B is %o\n", b;
    if not b then return false, false; end if;

    EK := BaseChange(E, K);
    MW, iota := MordellWeilGroup(E);

    fprintf filename, "E = %o\n", E;
    fprintf filename, "K defined by %o\n", MinimalPolynomial(K.1);
    fprintf filename, "****************************************************************************************\n";
    print "E(Q) = ", MW, " ---- ", [iota(g) : g in Generators(MW)];
    fprintf filename, "E(Q) = %o\n ------- %o\n", MW, [iota(g) : g in Generators(MW)];
    Qpts := [ iota(g) : g in Generators(MW) | Order(g) eq 0 ];
    pts := [ EK ! [K!x : x in ElementToSequence(iota(g))] : g in Generators(MW) | Order(g) eq 0 ];
    reg := MyRegulator(pts);
    r := #pts;
    print "The rank of E(Q) is ", r;

    print "Compute resolvents, periods and L-values";
    /* Compute the resolvents */
    Res := Resolvents(QG, h, delta);
    fprintf filename, "Resolvents =\n%o\n", Res;

    /* Compute the real and purely imaginary period of E/Q */
    Ovals := CompOvals(QG, E, OK, h);
    fprintf filename, "Periods =\n%o\n", Ovals;


    /* Compute the leading terms of the L-series and the order of vanishing */
    SetVerbose("LSeries", 1);
    Lvec := InitLSeries(QG, E, K, prec);
    Lvals,orders := Evaluate(Lvec);
    E`Lvals := Lvals; E`orders := orders; E`OK := OK; E`h := h;
    print "The analytic rank conjecturally is ", orders;
    fprintf filename, "The analytic rank conjecturally is %o\n", orders;
    fprintf filename, "Leading terms =\n%o\n", Lvals;

    if not AllNonTrivialOrdersAreZero(orders) then
        print "Not all analytic orders for non-trivial characters are 0";
        fprintf filename, "Not all analytic orders for non-trivial characters are 0\n";
        return false, false;
    end if;

    if CheckSaturation then
        fprintf filename, "V is %o-saturated : %o\n", l, IslSaturated(pts, l);
    else
        fprintf filename, "l-Saturation of V not checked\n";
    end if;


    RR := [ [reg] ] cat [ [1 : chi in QG`X[j]] : j in [2..#QG`X] ];
    print "Compute the order of Sha(E/Q)";

    ts := AnalyticTateShafarevicGroup(E, QG, Lvals, RR, OK);
    ts_gerundet := Runde(ts, Floor(2*prec/3));

    fprintf filename, "Analytic order of Sha times [E(K) : <gens>]^2 = %o = %o\n", ts, Factorization(ts_gerundet);
    print "Conjectural order of Sha = ", ts, " = ", Factorization(ts_gerundet);
    if GCD(l, ts_gerundet) ne 1 then
        print "The order of Sha is conjecturally divisible by l => CANNOT check the precise order of vanishing";
        fprintf filename, "The order of Sha is conjecturally divisible by l  => CANNOT check the precise order of vanishing";
        check_precise_ord := false;
    else
        print "The order of Sha is conjecturally NOT divisible by l => check the precise order of vanishing";
        fprintf filename, "The order of Sha is conjecturally NOT divisible by l  => check the precise order of vanishing\n";
        check_precise_ord := true;
    end if;

    i := 1;
    g0 := G.1;
    /* Note that the delta_psi is included in the regulator */
    Rvals := [ [i*reg/l^(r*n)] ] cat
             [ [(Conjugates(chi(g0))[1] - 1)^r : chi in QG`X[j]] : j in [2..#QG`X] ];
    Z := ComputeZ(QG, E, Lvals, Res, Ovals, Rvals);
    print "Z = ", Z;
    fprintf filename, "Z = %o\n", Z;

    RatZ := MakeRational(Z, QG);
    print "RatZ = ", RatZ;
    fprintf filename, "RatZ = %o\n", RatZ;

    eul := DualLocalEulerFactors(QG, E, E`OK, E`h, p);
    print "eul = ", eul;
    fprintf filename, "eul = %o\n", eul;

    Lstar := [ [ Z[i,j]*Conjugates(eul[i,j])[1] : j in [1..#Z[i]] ] : i in [1..#Z] ];
    print "Lstar = ", Lstar;
    fprintf filename, "Lstar = %o\n", Lstar;

    RatLstar := MakeRational(Lstar, QG);
    print "RatLstar = ", RatLstar;
    fprintf filename, "RatLstar = %o\n", RatLstar;

    u := RatLstar;

    /* The element u is pricisely the u as in the proof of Corollary 2.11 (whereas in the
       intrinsics DoWuthrichTypeExample and DolSquaredExample the u's differ by the delta_psi
    */

    sig := [r] cat [0 : i in [1..n]];
    a := ShaAnnihilator(QG, sig, u);
    print "Annihilator = ", a;
    g0 := G.1;
    fprintf filename, "Ann = ";
    for i:=0 to #G-1 do
        print Coefficient(a, g0^i), "*S^", i, " + ";
        if i eq #G-1 then
             fprintf filename, "%o*S^%o\n", Coefficient(a, g0^i), i;
        else
            fprintf filename, "%o*S^%o + ", Coefficient(a, g0^i), i;
        end if;
    end for;

    h := &+sig;;
    sigma := QG ! G.1;
    ZG := ModuleInit( IdentityMatrix(Rationals(), #G) );
    I := ModuleScalarProd(QG, ZG, (sigma - 1)^h);
    J := ModuleScalarProd(QG, ZG, (sigma - 1)^(h+1));
    a_in_I := IsModuleElt(QG, a, I);
    a_in_J := IsModuleElt(QG, a, J);
    fprintf filename, "Ann in I_{G,p}^%o ---------> %o\n", h, a_in_I;
    fprintf filename, "Ann in I_{G,p}^%o ---------> %o\n", h+1, a_in_J;
    printf "Ann in I_{G,p}^%o ---------> %o\n", h, a_in_I;
    printf "Ann in I_{G,p}^%o ---------> %o\n", h+1, a_in_J;

    if check_precise_ord then
        ETNCvalid := a_in_I and not a_in_J;
    else
        ETNCvalid := a_in_I;
    end if;

    /* We also check the part (iii) of Theorem 2.8, also we do not include these computations in the paper. */
    if n eq 2 then
        OC2 := MaximalOrder( Parent(RatZ[3,1]) );
        P2 := Factorization(l*OC2)[1,1];

        OC1 := MaximalOrder( Parent(RatZ[2,1]) );
        P1 := Factorization(l*OC1)[1,1];

        print "val(RatZ[1,1]) = ", Valuation(Rationals()!RatZ[1,1], l);
        print "val(RatZ[2,1]) = ", Valuation(RatZ[2,1], P1);
        print "val(RatZ[3,1]) = ", Valuation(RatZ[3,1], P2);
        fprintf filename, "val(RatZ[1,1]) = %o\n", Valuation(Rationals()!RatZ[1,1], l);
        fprintf filename, "val(RatZ[2,1]) = %o\n", Valuation(RatZ[2,1], P1);
        fprintf filename, "val(RatZ[3,1]) = %o\n", Valuation(RatZ[3,1], P2);

        if check_precise_ord then
            ETNCvalid := ETNCvalid and Valuation(Rationals()!RatZ[1,1], l) eq 0  
                                   and Valuation(RatZ[2,1], P1) eq 0 
                                   and Valuation(RatZ[3,1], P2) eq 0;
        end if;
    end if;
    
    fprintf filename, "\n";
    return true, ETNCvalid;
end intrinsic;



/* This is essentially the MAGMA command Points. However, we first take the first curve, say A, in the isogeny class (which should be the
   Strong Weil curve and shound also behave better algorithmically) and compute the points in A(K). Then we use the isogeny to compute
   points in E(K).
*/
intrinsic MyPoints(E :: CrvEll, EK :: CrvEll, Bound :: RngIntElt) -> SeqEnum
{}
    local a, f, AK, Kpts, phi, psi, omega, E_Kpts;

    A := IsogenousCurves(E)[1];
    if A eq E then
        Kpts := Points(EK : Bound:=Bound);
	return Kpts;
    else
        _, f := IsIsogenous(A, E);
        AK := BaseChange(A, BaseField(EK));
        Kpts := Points(AK : Bound:=Bound);

        phi := IsogenyMapPhi(f);
        psi := IsogenyMapPsi(f);
        psi2 := IsogenyMapPsiSquared(f); 
        xNum := psi2;
        yNum := psi*psi2;
        omega := IsogenyMapOmega(f);

        E_Kpts := [ EK ! [ Evaluate(phi, P[1]) / Evaluate(xNum, P[1]), Evaluate(omega, [P[1], P[2]]) / Evaluate(yNum, P[1]) ] : 
                    P in Kpts | Evaluate(psi, P[1]) ne 0 and P ne AK!0 ];

        return E_Kpts;
    end if;
end intrinsic;





/* These are naive routines that check l-saturation for a prime l.
   We assume that the order of the torsion subgroup is prime to l.
   gens is a list of generators, each of infinite order.

   This works only for l in {3,5}.
*/
intrinsic IslSaturated(gens, l) -> BoolElt
{}
    local i;

    S := CartesianProduct( [ [0..(l-1)] : i in [1..#gens]  ]);
    for s in S do
        P := &+[s[i]*gens[i] : i in [1..#gens]];
        if P ne Parent(gens[1])!0 then
            bool, Q := IsDivisibleBy(P, l);
            print "s = ", s, "    bool = ", bool;
            if bool then
                return false;
            end if;
        end if;
    end for;
    return true;
end intrinsic;


intrinsic IslSaturatedFast(gens, l) -> BoolElt
{}
    local i;

    bool := false;
    if IsDivisibleBy(gens[#gens], l) then return false; end if;
    j := #gens - 1;
    while j gt 1 do
        print "j = ", j;
        S := CartesianProduct( [ [0..(l-1)] : i in [j+1..#gens]  ]);
        c := [0 : g in gens ];
        c[j] := 1;
        for s in S do
            for i:=1 to #s do c[j+i]:=s[i]; end for;
            P := &+[c[i]*gens[i] : i in [1..#gens]];
            if P ne Parent(gens[1])!0 then
                bool, Q := IsDivisibleBy(P, l);
                print "c = ", c, "    bool = ", bool;
                if bool then
                    return false;
                end if;
            end if;
        end for;
        j := j-1;
    end while;
    return true;
end intrinsic;


/* Again a very naive routine that saturates at l. */
intrinsic lSaturation(gens, l) -> SeqEnum
{}
    local i;

    Z := gens;
    S := CartesianProduct( [ [0..(l-1)] : i in [1..#gens]  ]);
    Saturated := false;
    NrOfLoops := 0;
    while not Saturated do
        NrOfLoops := NrOfLoops+1;
        print NrOfLoops, ".loop";
        Saturated := true;
        for s in S do
            P := &+[s[i]*Z[i] : i in [1..#gens]];
            if P ne Parent(gens[1])!0 then
                bool, Q := IsDivisibleBy(P, l);
                if bool then
                    print "s = ", s, "  not yet saturated";
                    d, i0 := Max([Abs(l) : l in s ]);
                    for i:=1 to #gens do
                        if Abs(s[i]) le d and GCD(s[i], l) eq 1 then
                            i0 := i; d := Abs(s[i]);
                        end if;
                    end for;
                    // i0 := [i : i in [1..#gens] | GCD(s[i], l) eq 1][1];
                    Z[i0] := Q;
                    Saturated := false;
                    break;
                end if;
            end if;
        end for;
    end while;
    return Z;
end intrinsic;


/* The next two functions compute the element calL of Cor. 2.11 and Th. 2.12. Recall that calL is expected to annihilate Sha. 
   u is essentially leading terms divided by periods and regulator. Depending on the context of the computation, 
   u_psi is already divided by delta_psi or not.
*/

/* Use this function if the element u is not yet divided by the delta_psi 
   This does not coincide with the notation used in the paper.
*/
intrinsic ShaAnnihilator(QG :: AlgGrp, t0 :: RngIntElt, u :: List) -> AlgGrpElt
{}
    local a, G, c, i, j, chi;

    a := QG ! 0;
    G := Group(QG);
    for g in G do
        c := 0;
        for i in [(t0+1)..#QG`X] do
             for j:=1 to #QG`X[i] do
                 chi := QG`X[i,j];
                 c := c + u[i,j] * chi(g^(-1));
             end for;
        end for;
        a := a + (Rationals() ! (c / #G)) * (QG ! g);
    end for;
    print "a = ", a;
    return a;
end intrinsic;

/* Use this function if u is as in the proof of Corollary 2.11. sig is the
   list [m_0, ..., m_n].
*/
intrinsic ShaAnnihilator(QG :: AlgGrp, sig :: SeqEnum, u :: List) -> AlgGrpElt
{}
    local a, G, c, i, j, chi, delta, t0, g0, p;

    a := QG ! 0;
    G := Group(QG);
    p := PrimeDivisors(#G)[1];
    for g in G do
        c := 0;
        for i in [1..#QG`X] do
             for j:=1 to #QG`X[i] do
                 chi := QG`X[i,j];
                 c := c + u[i,j] * chi(g^(-1));
             end for;
        end for;
        a := a + (Rationals() ! (c / #G)) * (QG ! g);
    end for;
    t0 := #sig - 1;
    while sig[t0] eq 0 do t0 := t0 - 1; end while;
    g0 := QG ! G.1;
    print "sig = ", sig, "  t0 = ", t0;
    delta := &*[ (g0^(p^(j-1)) - 1)^sig[j] : j in [1..t0] ]; 
    a := a * delta;
    print "a = ", a;

    return a;
end intrinsic;







/******************************** WuthrichType  ***************************************
 *
 * WuthrichType only checks whether the pair (E, K) with K the extension of degree l in Q(\zeta_p) is such that
 * E(K)_p \simeq \Zl^{n_0} \oplus \ZlG. To that end we assume that rank(E(Q)) > 1 (not checked) so that we only need to check
 * that the analytic order of Lvec[2,1] is one.
 * If the relevant conjectures are true, this should be enough (and in all examples it is enough).
 * We also check whether the other assumptions of etncIV are satisfied.
 * We only consider the case n = 1, i.e., [K:Q] = l.
 *
 * This is used to find examples for etncIV and very preliminary. We assume from the beginning that rank(E(Q)) > 1.
 *
 *******************************************************************************************/

intrinsic WuthrichType(E :: CrvEll, l :: RngIntElt, n :: RngIntElt, p :: RngIntElt, prec :: RngIntElt) -> BoolElt, BoolElt, CrvEll, FldNum
{}
    local C, zeta, w, delta, f, L, A, K, OK, G, QG, act, h, EK, MW, iota, P, reg,
          R, Ovals, Lvec, Lvals, orders, Krel, g0, idem, e0, calE, calU, i, Z, d, RatZ, u, log,
          lambda;

    delete E`ts;           /* just to be sure */

    /* Compute the unique subfield L of Q(zeta_p) of degree l^n together with a integral normal
       basis element delta
    */
    print "Compute the field K and an integral normal basis element";
    C<zeta> := CyclotomicField(p);
    w := PrimitiveRoot(p);
    w := w^(l^n) mod p;
    delta1 := &+[zeta^(w^i) : i in [0..((p-1)/l^n)-1]];
    K := sub<C | delta1>;

    /* Compute Artin representations etc. */
    print "Compute the Artin representations";
    SetDefaultRealFieldPrecision(prec);
    A := ArtinRepresentations(K);
    G := K`artinrepdata`G;
    h := CompatibleIso(K, G);
    OK := MaximalOrder(K);
    QG := InitGroupAlgebra(G);
    delta := [h(g)(delta1) : g in G | Abs( Conjugates(h(g)(K!delta1))[1] - Conjugates(delta1)[1]  ) le 10^-10 ][1];

    print "Check condition Burns' (B)";
    ETNCIsValid := false;
    ESatisfiesB := ConditionB(E, OK, l);
    if not ESatisfiesB then
        print "Condition B not satisfied";
        return ESatisfiesB, ETNCIsValid, E, K;
    end if;

    /* Compute the Mordell-Weil group and the regulator.
    */
    print "Compute the Modell-Weil group E(Q)";
    EK := BaseChange(E, K);
    MW, iota := MordellWeilGroup(E);
    pts := [ EK ! [K!x : x in ElementToSequence(iota(g))] : g in Generators(MW) | Order(g) eq 0 ];
    reg := MyRegulator(pts);
    r := #pts;
    print "The rank of E(Q) is ", r;

    print "Compute resolvents, periods and L-values";
    /* Compute the resolvents */
    R := Resolvents(QG, h, delta);

    /* Compute the real and purely imaginary period of E/Q */
    Ovals := CompOvals(QG, E, OK, h);

    /* Compute the leading terms of the L-series and the order of vanishing */
    SetVerbose("LSeries", 1);
    Lvec := InitLSeries(QG, E, K, prec);

    val, o := LeadingTerm(Lvec[#Lvec, 1]);
    if o eq 1 then
        print "Example found    p = ", p, "   l = ", l, " n = ", n;
        return true, true, E, K;
    else
        print "Not all non trivial analytic orders are 1";
        return false, false, E, K;
    end if;
end intrinsic;

/******************************** lSquaredExamples  ***************************************
 *
 * lSquaredExamples only checks whether the pair (E, K) with K the extension of degree l^2 in Q(\zeta_p) is such that
 * E(K)_p \simeq \Zl^{n_0} \oplus \Zl[G/H]^{n_1} \oplus \ZlG^{n_2}. 
 * To that end we compute the order of Lvec[2,1] and Lvec[3,1]. By standard conjectures this should be enough.
 * We also check whether the other assumptions of etncIV are satisfied.
 * We only consider the case n = 2, i.e., [K:Q] = l^2.
 *
 * This is used to find examples for etncIV and very preliminary. We assume from the beginning that rank(E(Q)) >= 1.
 *
 *******************************************************************************************/

intrinsic lSquaredExamples(E :: CrvEll, l :: RngIntElt, n :: RngIntElt, p :: RngIntElt, prec :: RngIntElt) -> List, BoolElt
{}
    local C, zeta, w, delta, f, L, A, K, OK, G, QG, act, h, EK, MW, iota, P, reg,
          R, Ovals, Lvec, Lvals, orders, Krel, g0, idem, e0, calE, calU, i, Z, d, RatZ, u, log,
          lambda;

    delete E`ts;           /* just to be sure */

    l := 3; n := 2;
    /* Compute the unique subfield L of Q(zeta_p) of degree l^n together with a integral normal
       basis element delta
    */
    SetDefaultRealFieldPrecision(prec);
    print "Compute the field K and an integral normal basis element";
    C<zeta> := CyclotomicField(p);
    w := PrimitiveRoot(p);
    w := w^(l^n) mod p;
    delta1 := &+[zeta^(w^i) : i in [0..((p-1)/l^n)-1]];
    K := sub<C | delta1>;

    /* Compute Artin representations etc. */
    print "Compute the Artin representations";
    A := ArtinRepresentations(K);
    G := K`artinrepdata`G;
    h := CompatibleIso(K, G);
    OK := MaximalOrder(K);
    QG := InitGroupAlgebra(G);
    delta := [h(g)(delta1) : g in G | Abs( Conjugates(h(g)(K!delta1))[1] - Conjugates(delta1)[1]  ) le 10^-10 ][1];

    print "Check condition Burns' (B)";
    ETNCIsValid := false;
    ESatisfiesB := ConditionB(E, OK, l);
    if not ESatisfiesB then
        print "Condition B not satisfied";
        return [* E, p,l,n,[] *], false;
    end if;

    /* Compute the Mordell-Weil group and the regulator.
    */
    print "Compute the Modell-Weil group E(Q)";
    EK := BaseChange(E, K);
    MW, iota := MordellWeilGroup(E);
    pts := [ EK ! [K!x : x in ElementToSequence(iota(g))] : g in Generators(MW) | Order(g) eq 0 ];
    reg := MyRegulator(pts);
    r := #pts;
    print "The rank of E(Q) is ", r;

    /* Compute the leading terms of the L-series and the order of vanishing */
    SetVerbose("LSeries", 1);
    Lvec := InitLSeries(QG, E, K, prec);

    orders := [r];
    val, o := LeadingTerm(Lvec[2, 1]);
    if o eq 0 then
        return [* E, p, l, n, [r,0,0] *], false;
    else
        Append(~orders, o);
        val, o := LeadingTerm(Lvec[3, 1]);
        Append(~orders, o);
        return [* E, p, l, n, orders *],true;
    end if;
end intrinsic;

intrinsic NaiveLogModTorsion(Qpts :: SeqEnum, P :: PtEll, B :: RngIntElt) -> List
{}

    local i;

    S := CartesianProduct( [ [-B..B] : i in [1..#Qpts]  ]);

    for s in S do
        Q := &+[ s[i]*Qpts[i] : i in [1..#s] ];
        if Order(Q - P) ne 0 then
            return [t : t in s];
            // print [t : t in s];
        end if;
    end for;
    return [];
end intrinsic;

intrinsic Trace(P ::  PtEll, h :: Map) -> PtEll
{}

    local G;

    return &+[ GAction(g, P, h) : g in Domain(h) ];

end intrinsic;


