declare attributes CrvEll : Lvals, Rvals, orders, RatZ, EK, OK, K, h, ts;

/******************************** etnc ***************************************
 *
 * pol is an irreducible normalized polynomial over Z which will define a splitting field K
 * E is an elliptic curve defined over Q
 * prec is a positive integer which gives the precision for the L-series computations
 *
 * The function tries to verify ETNC numerically. 
 *
 *******************************************************************************************/

intrinsic etnc(pol :: RngUPolElt, E :: CrvEll, prec :: RngIntElt) -> CrvEll, SetEnum, SetEnum, AlgGrp
{}
    local K, OK, f, X, theta, Atheta, AssOrd, nb, h, QH, Idem, G, lambda, delta, e, F,
          IrrMod, IrrRep, IrrChar, R, Ovals, Lvec, Lvals, Z, S;

    delete E`ts;           /* just to be sure */
    
    /* Compute the splitting field K and the Artin representations of K/Q */
    print "Compute the splitting field K and the Artin representations of K/Q";
/*  Does not work this way in newer MAGMA versions
    L := NumberField(pol);
    A := ArtinRepresentations(L);
    K :=L`artinrepdata`K;
*/    
    K := SplittingField(NumberField(pol));
    pol := MinimalPolynomial(K.1);
    L := NumberField(pol);
    A := ArtinRepresentations(L);
    K :=L`artinrepdata`K;
    assert IsTamelyRamified(K);


    /* Compute an integral normal basis element. We always assume that K/Q is tame.*/
    print "Compute an integral normal basis element delta";
    X, isfree := INB(K);
    K := X[1]; theta := X[2]; Atheta := X[3]; AssOrd := X[4]; nb := X[5]; h := X[6]; QG := X[7];
    OK := MaximalOrder(K);
    G := Domain(h);
    lambda := QG ! ElementToSequence(nb);
    delta := QGAction(lambda, theta, h);
    print "delta = ", delta;

    QG := InitGroupAlgebra(G);

    print "Compute resolvents, periods and L-values";
    /* Compute the resolvents of the integral normal basis element delta */
    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);
    Lvals,orders := Evaluate(Lvec);
    E`Lvals := Lvals; E`orders := orders; E`OK := OK; E`h := h;
    print "The analytic rank conjecturally is ", orders;

    Rvals := [ [ComplexField()!1 : i in [1..#z]] : z in Lvals];
    for i:=1 to #orders do
        for j:=1 to #orders[i] do
            if orders[i, j] gt 0 then
	        print "Rvals[", i, " ,", j, "] = ";
	        Rvals[i, j] := ReadReal("Regulator = ");
	    end if;
	end for;
    end for;


    /* Compute the ratios of L-values, periods and resolvents and check rationality */
    print "Compute approximations to the twisted BSD quotients";
    Z := ComputeZ(QG, E, Lvals, R, Ovals, Rvals);
    print "Z = ", Z;
    RatZ, fehler := MakeRational(Z, QG);
    print "Error = ", fehler;
    print "The rounded values are conjecturally: ", RatZ;

    E`Rvals := Rvals; E`K := NumberField(E`OK); E`RatZ := RatZ; E`EK := BaseChange(E, E`K); 
    S := ComputeS(E, OK); /* the set S of bad primes; those where the representation is ramified */
    print "Compute a conjectural value for the order of Sha(E/K)";
    E`ts := TateShafarevicGroup(E, QG, Lvals, Rvals, OK);
    print "Conjecturally #Sha(E/K) = ", E`ts;
    HP := S join HardPrimes(E, QG, Lvals, Rvals, OK);
    print "HP = ", HP;

    /* Check primes which are not in HP */
    if CheckEasyPrimes(RatZ, HP) then
        print "ETNC conjecturally true outside ", HP;
    else
        print "ETNC is false !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!";
    end if;

    /* Try to check the hard primes */
    print "Now check the hard primes";
    PrimesThatCannotBeChecked := [];
    for l in HP do
         print "\n\n******************* l = ", l, " *************************************";
         IsComputable, ETNCIsValid := CheckHardPrime(QG, E, S, l);
         if IsComputable then
             if ETNCIsValid then
	         print "ETNC is valid for l = ", l;
	     else
	         print "Counter example to ETNC for E = ", E, "  l = ", l, " K = ", K;
                 assert false;
	     end if;
	 else
	     print "Cannot check ETNC for l = ", l;
             Append(~PrimesThatCannotBeChecked, l); 
         end if;
    end for;
    print "ETNC is checked for all primes not in ", Set(PrimesThatCannotBeChecked);
    return E, S, HP, QG;
end intrinsic;


intrinsic NonPerfectEtnc(E :: CrvEll, l :: RngIntElt, n :: RngIntElt, p :: RngIntElt, prec :: RngIntElt) -> BoolElt, BoolElt, SetEnum, SetEnum, 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, {1}, {}, K;
    end if;

    /* Compute the Mordell-Weil group and the regulator. 
        irint c;
    */
    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);
    Lvals,orders := Evaluate(Lvec);
    E`Lvals := Lvals; E`orders := orders; E`OK := OK; E`h := h;
    print "The analytic rank conjecturally is ", orders;

    if not AllNonTrivialOrdersAreZero(orders) then
        print "Not all analytic orders for non-trivial characters are 0";
        return false, false, {1}, {}, K;
    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 := TateShafarevicGroup(E, QG, Lvals, RR, OK);
    print "Conjectural order of Sha = ", ts;
    if GCD(l, ts) ne 1 then 
        print "The order of Sha is conjecturally divisible by l Hurra, Hurrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrra";
        return false, false, {1}, {}, K;
    end if;
    EKtor := TorsionSubgroup(EK); Etor := TorsionSubgroup(E);
    if #EKtor ne #Etor then 
        print "E(Q)_tors is not equal to E(K)_tors";
        return false, false, {1}, {}, K;
    end if;

    print "Compute the relative group";
    Krel := RelativeGroup(QG, l);
    g0 := G.1;
    idem := QIdempotents(QG, QG`X, QG`H);
    e0 := idem[1];
    calU := []; calE := [];

    if r eq 0 then M := 1; else M := l^n-1; end if;
    eul := DualLocalEulerFactors(QG, E, E`OK, E`h, p);
    u := [* z[1] : z in eul *];
    log := K0RelLog(Krel, u);
    print "Log of xi:   ", log;

    print "Compute the sets calE and calC";
    for i:=1 to M do
        if GCD(i, l)  eq 1 then
            // print "************************* i = ", i;
            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, R, Ovals, Rvals);
            d := Denominator( BestApproximation(Z[1,1], 100000) );
            RatZ := MakeRational(Z, QG, d, Round(30));
            u := [* RatZ[i,1] * eul[i,1] : i in [1..#RatZ] *];
            log := K0RelLog(Krel, u);
            if not IdealPartIsTrivial(log) then
	        print "Non trivial ideal part --> ETNC not valid mod torsion !!!!!!!!!!!!!!";
		return true, false, {0}, {}, K;
	    end if;
            // print "i = ", i, "   log = ", log;
            Append(~calU, log[2]);
            // lambda:=e0 + (One(QG) - e0) * &+[(QG!g0_L)^j : j in [0..i-1]];
            lambda:=i*e0 + (One(QG) - e0);
            log := K0RelLog(Krel, ReducedNorm(lambda));
            Append(~calE, log[2]);
         end if;
    end for;

    calU := Set(calU); calE := Set(calE);
    return true, calU eq calE, calU, calE, K;

end intrinsic;

/* Only for testing; not used at the moment, see DoNonPerfectAnnihilator instead. */

intrinsic NonPerfectAnnihilator(E :: CrvEll, l :: RngIntElt, n :: RngIntElt, p :: RngIntElt, prec :: RngIntElt) -> BoolElt, AlgGrpElt
{}
    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 */


    print "E = ", CremonaReference(E), " l = ", l, " n = ", n, " p = ", p, " prec = ", prec;
    /* 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, {1}, {}, K;
    end if;

    /* Compute the Mordell-Weil group and the regulator.
        irint c;
    */
    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);
    Lvals,orders := Evaluate(Lvec);
    E`Lvals := Lvals; E`orders := orders; E`OK := OK; E`h := h;
    print "The analytic rank conjecturally is ", orders;

    if not AllNonTrivialOrdersAreZero(orders) then
        print "Not all analytic orders for non-trivial characters are 0";
        return false, false, {1}, {}, K;
    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 := TateShafarevicGroup(E, QG, Lvals, RR, OK);
    print "Conjectural order of Sha = ", ts;
    if GCD(l, ts) ne 1 then
        print "The order of Sha is conjecturally divisible by l Hurra, Hurrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrrra";
    end if;
    EKtor := TorsionSubgroup(EK); Etor := TorsionSubgroup(E);
    if #EKtor ne #Etor then
        print "E(Q)_tors is not equal to E(K)_tors";
        return false, false, {1}, {}, K;
    end if;

    eul := DualLocalEulerFactors(QG, E, E`OK, E`h, p);
    i := 1;
    g0 := G.1;
    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, R, Ovals, Rvals);
    d := Denominator( BestApproximation(Z[1,1], 100000) );
    RatZ := MakeRational(Z, QG, d, Round(10));
    u := [* [ RatZ[i,j] * eul[i,j] : j in [1..#RatZ[i]] ] : i in [1..#RatZ] *];
    a := ShaAnnihilator(QG, [1,0], u);

    print "Annihilator = ", Eltseq( a );
    return true, a;

end intrinsic;

/* No longer needed and not used in the moment. */
intrinsic AnnihilationExample(E :: CrvEll, l :: RngIntElt, n :: RngIntElt, p :: RngIntElt, prec :: RngIntElt) -> 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 */

    /* 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, {1}, {}, K;
    end if;

    /* Compute the Mordell-Weil group and the regulator.
    */
    print "Compute the Modell-Weil group E(Q)";

    ZX<x> := PolynomialRing(Integers());
    Q := NumberField(x-1 : DoLinearExtension:=true);
    EQ := BaseChange(E, Q);
    EK := BaseChange(E, K);
    MW, iota := MordellWeilGroup(E);
    Kpts := [ EK ! [K!x : x in ElementToSequence(iota(g))] : g in Generators(MW) | Order(g) eq 0 ];
    Qpts := [ EQ ! [K!x : x in ElementToSequence(iota(g))] : g in Generators(MW) | Order(g) eq 0 ];

    Q_ts := Integers() ! Round( ConjecturalSha(EQ, Set(Qpts)) ); 
    K_ts := Integers() ! Round( ConjecturalSha(EK, Set(Kpts)) );

    print "l = ", l, "     K_ts = ", K_ts, "Q_ts = ", Q_ts;
    if GCD(l, Q_ts) eq 1 and GCD(l, K_ts) gt 1 then
        print "Example found ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++";
        return [* E, p, l, n, [] *],true;
    else
        return [* *], false; 
    end if;
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;






/******************************** PerfectEtnc  ***************************************
 *
 * PerfectEtnc 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 \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 enouhg (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 PerfectEtnc(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;


/******************************** 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;



/* This routine is used to check the full ETNC at all primes for the examples as described in etncII, Section 6.
   We can only check the ETNC at l if E(K)_l is perfect, which in the context of etncII means that rk(E(K)) = 0. 
   We first check ETNC at l and then use the methods of etncI to deal with the primes different from l.
*/
intrinsic etncII(E :: CrvEll, l :: RngIntElt, n :: RngIntElt, p :: RngIntElt, prec :: RngIntElt) -> BoolElt, BoolElt, SeqEnum
{}
    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 K and an integral normal basis element";
    C<zeta> := CyclotomicField(p);
    w := PrimitiveRoot(p);
    w := w^(l^n) mod p;
    delta := &+[zeta^(w^i) : i in [0..((p-1)/l^n)-1]];
    f := MinimalPolynomial(delta);
    L := NumberField(f);


    /* Compute Artin representations etc. */
    print "Compute the Artin representations";
    SetDefaultRealFieldPrecision(prec);
    A := ArtinRepresentations(L);
    K :=L`artinrepdata`K;
    OK := MaximalOrder(K);
    G := L`artinrepdata`G;
    QG := InitGroupAlgebra(G);
    act := L`artinrepdata`act;
    h := map<G->Parent(act(G.1)) | g:->act(g^-1)>;
    delta := K.1;

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

    /* Compute the Mordell-Weil group and the regulator. 
    */
    print "Compute the Mordell-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;
    assert r eq 0;
	 
    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);
    Lvals,orders := Evaluate(Lvec);
    E`Lvals := Lvals; E`orders := orders; E`OK := OK; E`h := h;
    print "The analytic rank conjecturally is ", orders;

    if not AllNonTrivialOrdersAreZero(orders) then
        print "Not all analytic orders for non-trivial characters are 0";
        return false, false, {1}, {}, K;
    end if;
    Rvals := [ [ComplexField()!1 : i in [1..#z]] : z in Lvals];

    print "Compute a conjectural value for the order of Sha";
    E`ts := TateShafarevicGroup(E, QG, Lvals, Rvals, OK);
    print "Conjecturally #Sha(E/K) = ", E`ts;

    if GCD(l, E`ts) ne 1 then 
        print "The order of Sha is conjecturally divisible by l";
        return false, false, {1}, {}, K;
    end if;
    EKtor := TorsionSubgroup(EK); Etor := TorsionSubgroup(E);
    if #EKtor ne #Etor then 
        print "E(Q)_tors is not equal to E(K)_tors";
        return false, false, {1}, {}, K;
    end if;

    print "Compute the relative group";
    Krel := RelativeGroup(QG, l);
    g0 := G.1;
    idem := QIdempotents(QG, QG`X, QG`H);
    e0 := idem[1];
    calU := []; calE := [];

    eul := DualLocalEulerFactors(QG, E, E`OK, E`h, p);
    u := [* z[1] : z in eul *];
    log := K0RelLog(Krel, u);
    print "Log of xi:   ", log;
 
    Z := ComputeZ(QG, E, Lvals, R, Ovals, Rvals);
    print "Z = ", Z;
    RatZ, fehler := MakeRational(Z, QG);
    print "Error = ", fehler;
    print "The rounded values are conjecturally: ", RatZ;

    print "Check ETNC at l = ", l;
    u := [* RatZ[i,1] * eul[i,1] : i in [1..#RatZ] *];
    log := K0RelLog(Krel, u);
    ETNCIsValid := IdealPartIsTrivial(log) and (log[2] eq Id(Parent(log[2])));
    if not ETNCIsValid then
        print "Counter example to ETNC !!!!!!!!!!!!!!!!!!!!!!!!!!!";
	return true, false, [];
    end if;

    E`Rvals := Rvals; E`K := NumberField(E`OK); E`RatZ := RatZ; E`EK := BaseChange(E, E`K);
    S := ComputeS(E, OK); /* the set S of bad primes; those where the representation is ramified */
    HP := S join HardPrimes(E, QG, Lvals, Rvals, OK);
    print "HP = ", HP;

    /* Check primes which are not in HP */
    if CheckEasyPrimes(RatZ, HP) then
        print "ETNC conjecturally true outside ", HP;
    else
        print "ETNC is false !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!";
        return true, false, [];
    end if;

    /* Try to check the hard primes */
    print "Now check the hard primes";
    PrimesThatCannotBeChecked := [];
    for q in HP diff {l} do
         print "\n\n******************* q = ", q, " *************************************";
         IsComputable, ETNCIsValid := CheckHardPrime(QG, E, S, q);
         if IsComputable then
             if ETNCIsValid then
                 print "ETNC is valid for q = ", q;
             else
                 print "Counter example to ETNC for E = ", E, "  q = ", q, " K = ", K;
                 assert false;
             end if;
         else
             print "Cannot check ETNC for q = ", q;
             Append(~PrimesThatCannotBeChecked, q);
         end if;
    end for;
    print "ETNC is checked for all primes not in ", Set(PrimesThatCannotBeChecked);

    return true, true, PrimesThatCannotBeChecked;
end intrinsic;

/* Just for experiments, not used and superseded by the results in BMW. */    
intrinsic HeegnerEtnc(E :: CrvEll, D :: RngIntElt, prec :: RngIntElt) -> SeqEnum, SeqEnum
{}
    local i, j;

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

    SetDefaultRealFieldPrecision(prec);

    K := QuadraticField(D);
    P := HeegnerPoints(E,D);
    L<u> := NumberField(P[1]);
    PT := Points(ChangeRing(E,L),u)[1];
    H := AbsoluteField( HilbertClassField(K) );

    // X := INB(H : WeakINB:=true);
    X := INB(H);
    H := X[1]; theta := X[2]; Atheta := X[3]; AssOrd := X[4]; nb := X[5]; h := X[6]; QG := X[7];
    OH := MaximalOrder(H);
    // Idem := QIdempotents(QG, QG`X, QG`H);
    G := Domain(h);
    lambda := QG ! ElementToSequence(nb);
    delta := QGAction(lambda, theta, h);

    QG := InitGroupAlgebra(G);

    /* Compute the resolvents of the integral normal basis element delta */
    R := Resolvents(QG, h, delta);
    R;

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

    Lvec := InitLSeries(QG, E, H, prec);
    // Lvec := OldInitLSeries(QG, E, OH, h, prec);
    Lvals,orders := Evaluate(Lvec);
    E`Lvals := Lvals; E`orders := orders; E`OK := OH; E`h := h;
    print "orders = ", orders;

    _, j := IsSubfield(L, H);
    EH := BaseChange(E, H);
    Q := EH ! [j(x) : x in ElementToSequence(PT)];

    Rvals := [ [ComplexField()!1 : i in [1..#z]] : z in Lvals];
    for i:=1 to #orders do
        for j:=1 to #orders[i] do
	    if Degree(QG`X[i,j]) eq 1 then
	       K1 := FixedField(H, [h(g) : g in Kernel(QG`X[i,j])]);
               D1 := Discriminant( MaximalOrder( K1 ));
               Et := QuadraticTwist(E, D1);
	       print "Compute the Mordell-Weil group for the ", D1, "-twist";
               MW_t, iota_t:=MordellWeilGroup(Et);
               EK1 := BaseChange(E, K1);
               EtK1 := BaseChange(Et, K1);
               bool, f := IsIsomorphic(EtK1, EK1);
               pts := [ EH ! f(iota_t(g)) : g in Generators(MW_t) | Order(g) eq 0 ];
               print "i = ", i, " j = ", j, " Compute the regulator (deg(chi = 1)";
               Rvals[i,j] := MyRegulator(pts);
            end if;
            if Degree(QG`X[i,j]) eq 2 then
                print "i = ", i, " j = ", j, " Compute the regulator (deg(chi) = 2)";
                chi := QG`X[i, j]; 
                Rvals[i,j] := &+[HeightPairing(Q, GAction(Q, g, h), H) * Conjugates(chi(g^-1))[1] : g in G] / #G;
            end if;
        end for;
    end for;
     
    Z := ComputeZ(QG, E, Lvals, R, Ovals, Rvals);
    print "Z = ", Z;
    RatZ, fehler := MakeRational(Z, QG);
    print "Error = ", fehler;
    print "The rounded values are conjecturally: ", RatZ;

    E`Rvals := Rvals; E`K := NumberField(E`OK); E`RatZ := RatZ; E`EK := BaseChange(E, E`K);
    S := ComputeS(E, OH); /* the set S of bad primes; those where the representation is ramified */
    print "Compute a conjectural value for the order of Sha(E/K)";
    E`ts := TateShafarevicGroup(E, QG, Lvals, Rvals, OH);
    print "Conjecturally #Sha(E/K) = ", E`ts;
    HP := S join HardPrimes(E, QG, Lvals, Rvals, OH);
    print "HP = ", HP;

    /* Check primes which are not in HP */
    if CheckEasyPrimes(RatZ, HP) then
        print "ETNC conjecturally true outside ", HP;
    else
        print "ETNC is false !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!";
    end if;

    print "Now check the hard primes";
    PrimesThatCannotBeChecked := [];
    for l in HP do
         print "\n\n******************* l = ", l, " *************************************";
         IsComputable, ETNCIsValid := CheckHardPrime(QG, E, S, l);
         if IsComputable then
             if ETNCIsValid then
                 print "ETNC is valid for l = ", l;
             else
                 print "Counter example to ETNC for E = ", E, "  l = ", l, " H = ", H;
                 assert false;
             end if;
         else
             print "Cannot check ETNC for l = ", l;
             Append(~PrimesThatCannotBeChecked, l);
         end if;
    end for;
    print "ETNC is checked for all primes not in ", Set(PrimesThatCannotBeChecked);
    return E, S, HP, QG;
end intrinsic;


/******************************** AllNonTrivialOrdersAreZero ***************************************
 *
 * Self explanatory.
 *
 *******************************************************************************************/

intrinsic AllNonTrivialOrdersAreZero(orders :: SeqEnum) -> BoolElt
{}
    local i, j;

    for i := 2 to #orders do
        for j:= 1 to #orders[i] do
            if orders[i, j] ne 0 then
	        return false;
	    end if;
	end for;
    end for;
    return true;
end intrinsic;

/******************************** AllNonTrivialOrdersAreZero ***************************************
 *
 * Self explanatory.
 *
 *******************************************************************************************/

intrinsic AllNonTrivialOrdersAreOne(orders :: SeqEnum) -> BoolElt
{}
    local i, j;

    for i := 2 to #orders do
        for j:= 1 to #orders[i] do
            if orders[i, j] ne 1 then
                return false;
            end if;
        end for;
    end for;
    return true;
end intrinsic;


/******************************** AllNonTrivialOrdersAreZero ***************************************
 *
 * Self explanatory.
 *
 *******************************************************************************************/

intrinsic AllPrimitiveOrdersAreOne(orders :: SeqEnum) -> BoolElt
{}
    local i, j;

        last := #orders;
        for j:= 1 to #orders[last] do
            if orders[last, j] ne 1 then
                return false;
            end if;
        end for;
    return true;
end intrinsic;

/******************************** AllOrdersAreOne ***************************************
 *
 * Self explanatory.
 *
 *******************************************************************************************/

intrinsic AllOrdersAreOne(orders :: SeqEnum) -> BoolElt
{}
    local i, j;

    for i := 1 to #orders do
        for j:= 1 to #orders[i] do
            if orders[i, j] ne 1 then
                return false;
            end if;
        end for;
    end for;
    return true;
end intrinsic;

/******************************** SemistableEllCurves ***************************************
 *
 * Computes a list of [ L0, L1, L2, L3, L4, L5, L6 ], where each Lr is a list of
 * semistable elliptic curves of algebraic rank r. Each curve is given by a list of
 * the form [* CremonaReference(E), E *].
 *
 *******************************************************************************************/

intrinsic SemistableEllCurves(low :: RngIntElt, high :: RngIntElt) -> SeqEnum
{}
    local EllCurvesList, DB, N, NrOfCrvs, m, i, j, E, MW, iota, r;

    EllCurvesList := [[],[],[],[], [], [], []];
    DB := CremonaDatabase();

    for N:=low to high do
        NrOfCrvs := NumberOfCurves(DB, N);
        // print "N = ", N, " NumberOfCurves = ", NrOfCrvs;
        if NrOfCrvs gt 0 then
            m := NumberOfIsogenyClasses(DB, N);
            for i:=1 to m do
                for j:=1 to NumberOfCurves(DB,N, i) do
                    E := EllipticCurve(DB, N, i, j);
		    if IsSemistable(E) then
		        r := MordellWeilRank(E);
                        Append(~(EllCurvesList[r+1]), [* CremonaReference(E), E *]);
		    end if;
		end for;
	    end for;
	end if;
    end for;
    return EllCurvesList;
end intrinsic;

/******************************** SplitMultEllCurves ***************************************
 *
 * Computes a list of [ L0, L1, L2, L3, L4, L5, L6 ], where each Lr is a list of
 * split multiplicative  elliptic curves of algebraic rank r. Each curve is given by a list of
 * the form [* CremonaReference(E), E *].
 *
 *******************************************************************************************/

intrinsic SplitMultEllCurves(low :: RngIntElt, high :: RngIntElt) -> SeqEnum
{}
    local EllCurvesList, DB, N, NrOfCrvs, m, i, j, E, MW, iota, r;

    EllCurvesList := [[],[],[],[], [], [], []];
    DB := CremonaDatabase();

    for N:=low to high do
        NrOfCrvs := NumberOfCurves(DB, N);
        print "N = ", N, " NumberOfCurves = ", NrOfCrvs;
        if NrOfCrvs gt 0 then
            m := NumberOfIsogenyClasses(DB, N);
            for i:=1 to m do
                for j:=1 to NumberOfCurves(DB,N, i) do
                    E := EllipticCurve(DB, N, i, j);
                    if IsSplitMultiplicative(E) then
                        r := MordellWeilRank(E);
                        Append(~(EllCurvesList[r+1]), [* CremonaReference(E), E *]);
                    end if;
                end for;
            end for;
        end if;
    end for;
    return EllCurvesList;
end intrinsic;

/******************************** ExtClasses ***************************************
 *
 * No longer needed.
 *
 *******************************************************************************************/

intrinsic ExtClasses(l :: RngIntElt, n :: RngIntElt, p :: RngIntElt, prec :: RngIntElt, low :: RngIntElt, high :: RngIntElt) -> List
{}
    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, S1, S2, i, Z, d, RatZ, u, log,
          lambda, eul, Rat_eul;

    /* Compute the unique subfield L of Q(zeta_p) of degree l^n together with a integral normal
       basis element delta 
    */
    C<zeta> := CyclotomicField(p);
    w := PrimitiveRoot(p);
    w := w^(l^n) mod p;
    delta := &+[zeta^(w^i) : i in [0..((p-1)/l^n)-1]];
    f := MinimalPolynomial(delta);
    L := NumberField(f);

    /* Compute Artin representations etc. */
    SetDefaultRealFieldPrecision(prec);
    A := ArtinRepresentations(L);
    K :=L`artinrepdata`K;
    OK := MaximalOrder(K);
    G := L`artinrepdata`G;
    QG := InitGroupAlgebra(G);
    act := L`artinrepdata`act;
    h := map<G->Parent(act(G.1)) | g:->act(g^-1)>;
    delta := K.1;

    Krel := RelativeGroup(QG, l);
    g0 := G.1;
    idem := QIdempotents(QG, QG`X, QG`H);
    e0 := idem[1];

    ClassList := [**];
    fixed_rank := 1;
    DB := CremonaDatabase();

    for N:=low to high do
        NrOfCrvs := NumberOfCurves(DB, N);
        print "N = ", N, " NumberOfCurves = ", NrOfCrvs;
        if NrOfCrvs gt 0 then
            m := NumberOfIsogenyClasses(DB, N);
            for i:=1 to m do
                for j:=1 to NumberOfCurves(DB,N, i) do
                    E := EllipticCurve(DB, N, i, j);
                    if not IsSemistable(E) or not ConditionB(E, OK, l) then break; end if;
                    print "E = ", CremonaReference(E);
                    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 ];
                    r := #pts;
                    if r eq fixed_rank then
                        print "r = ", r;
			reg := MyRegulator(pts);
                        /* 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;
                        if AllNonTrivialOrdersAreZero(orders) then
                            print "The analytic rank is ", orders;

                            /* Compute the resolvents */
                            R := Resolvents(QG, h, delta);

                            /* Compute the real and purely imaginary period of E/Q */
                            Ovals := CompOvals(QG, E, OK, h);
                            RR := [ [reg] ] cat [ [1 : chi in QG`X[j]] : j in [2..#QG`X] ];
                            ts := TateShafarevicGroup(E, QG, Lvals, RR, OK);
                            if GCD(l, ts) ne 1 then break; end if;
                            EKtor := TorsionSubgroup(EK); Etor := TorsionSubgroup(E);
                            if #EKtor ne #Etor then break; end if;

                            for q:=1 to l^n-1 do
                                if GCD(q, l)  eq 1 then
                                    Rvals := [ [q*reg/l^(r*n)] ] cat
                                             [ [(Conjugates(chi(g0))[1] - 1)^r : chi in QG`X[j]] : j in [2..#QG`X] ];
                                    RR := [ [reg] ] cat
                                             [ [1 : chi in QG`X[j]] : j in [2..#QG`X] ];

                                    Z := ComputeZ(QG, E, Lvals, R, Ovals, Rvals);
                                    // print "Z = ", Z;
                                    // readi d, "Input a denominator";
                                    d := Denominator( BestApproximation(Z[1,1], 100000) );
                                    RatZ := MakeRational(Z, QG, d, Round(4));
                                    eul := DualLocalEulerFactors(QG, E, E`OK, E`h, p);
				    // u := [* z[1] : z in RatZ *];
				    u := [* RatZ[i,1] * eul[i,1] : i in [1..#RatZ] *]; 
                                    log := K0RelLog(Krel, u);
				    print "log = ", log;
                                    if log[2] eq Id(Krel`DT) then
                                        Append(~ClassList,  [* CremonaReference(E), q, LocalInformation(E) , MW, iota, Conductor(E), Discriminant(E) *]);
                                        break;
                                    end if;
                                 end if;
                             end for;
                         end if;
                     end if;
                end for;
            end for;
        end if;
    end for;
    
    return ClassList;

end intrinsic;

/******************************** ConditionB ***************************************
 *
 * This routine checks the hypothesis imposed in the introduction of etncII. We already assume that K/Q is abelian of l-power order.
 * Then we verify all the other hypothesis but the last one which requires l \nmid #Sha. A numerical verification (modulo the validity of BSD)
 * of this condition should be done later when all the periods, leading terms and regulators are computed by using BSD for E/K. 
 *
 *******************************************************************************************/

intrinsic ConditionB(E :: CrvEll, OK :: RngOrd, l :: RngIntElt) -> BoolElt
{}
    local i;

    N := Conductor(E);
    D := Discriminant(OK);
    if GCD(l*N, Discriminant(OK)) ne 1 then return false; end if;
    Etors := TorsionSubgroup(E);
    m := #Etors * &*[ NrOfNsPoints(E, q) : q in PrimeDivisors(D) ];
    if GCD(l, m) ne 1 then return false; end if;
    m := &*TamagawaNumbers(E);
    if GCD(l, m) ne 1 then return false; end if;
  
    return true;
end intrinsic;


/******************************** CompLvals ***************************************
 *
 * No longer used. If needed this function must be adapted as the beginning of the function 
 * etnc in order to use ArtinRepresentations etc.
 * !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 *
 *******************************************************************************************/

intrinsic CompLvals(pol :: RngUPolElt, E :: CrvEll, prec :: RngIntElt) -> List
{}
    local K, OK, f, X, theta, Atheta, AssOrd, nb, h, QH, Idem, G, lambda, delta, e, F,
          IrrMod, IrrRep, IrrChar, R, Ovals, Lvec, Lvals, Z, S;

    K := SplittingField(NumberField(pol));
    OK := MaximalOrder(K);


    K := SplittingField(NumberField(pol));
    pol := MinimalPolynomial(K.1);
    L := NumberField(pol);
    A := ArtinRepresentations(L);
    K :=L`artinrepdata`K;
    // assert IsTamelyRamified(K);

    print "Compute an integral normal basis element delta";
    X, isfree := INB(K);
    K := X[1]; theta := X[2]; Atheta := X[3]; AssOrd := X[4]; nb := X[5]; h := X[6]; QG := X[7];
    OK := MaximalOrder(K);
    G := Domain(h);
    lambda := QG ! ElementToSequence(nb);
    delta := QGAction(lambda, theta, h);
    print "delta = ", delta;

    QG := InitGroupAlgebra(G);

    /* Compute the resolvents of the integral normal basis element delta */
    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);
    Lvals,orders := Evaluate(Lvec);
    E`Lvals := Lvals; E`orders := orders; E`OK := OK; E`h := h;
    print "The analytic rank is ", orders;

    return [* E, pol, Lvals, orders *];
end intrinsic;


/******************************** AnalyticRank ***************************************
 *
 * Returns the conjectural order of the Hasse-Weil L-series L(E/K, s) at s=1.
 * orders are the conjectural analytic orders of the twisted Hasse-Weil-L-functions at s=1.
 *
 *******************************************************************************************/

intrinsic AnalyticRank(orders :: SeqEnum) -> RngIntElt
{}
    return &+[ &+o : o in orders];
end intrinsic;



/******************************** Frobenius ***************************************
 *
 * OK ring of ints in K, h map provided by AutomorphismGroup or ArtinRepresentations
 * h:G -> Aut(K/Q), p is a rational prime
 *
 * Computes  a lift of Frobenius of \frp / p in the Galois group G, where \frp is a chosen prime ideal
 * in K above p.
 *
 *******************************************************************************************/

intrinsic Frobenius(OK :: RngOrd, h :: Map, p :: RngIntElt) -> GrpPermElt
{}
    local G, b, is_frob, vals, i;

    G := Domain(h);

    P := Factorization(p*OK)[1,1];
    for g in G do
        vals := [Valuation( h(g)(b) - b^p, P ) : b in Basis(OK)];
	is_frob := true;
	i := 2;
	while i le #vals and is_frob do
	    if vals[i] lt 1 then 
		is_frob := false;
	    end if;
	    i :=i+1;
	end while;
	if is_frob then
	    return(g);
	end if;
    end for;
end intrinsic;


/******************************** ComplexConjugation ***************************************
 *
 * OK ring of ints in K, h map provided by AutomorphismGroup or ArtinRepresentations
 *
 * Computes  complex conjugation as an element of the Galois group G. For this we always assume
 * that an embedding of K into C was chosen once and for all.
 *
 *******************************************************************************************/

intrinsic ComplexConjugation(OK :: RngOrd, h :: Map) -> GrpPermElt
{}
    local G, b, is_c, a, i, prec;

    G := Domain(h);

    prec := 10;   /* very preliminary */
    for g in G do
        a := [Abs( Conjugates(h(g)(b))[1] - Conjugate(Conjugates(b)[1]) ) : b in Basis(OK)];
        is_c := true;
        i := 2;
        while i le #a and is_c do
            if a[i] gt 10^(-prec) then
                is_c := false;
            end if;
            i :=i+1;
        end while;
        if is_c then
            return(g);
        end if;
    end for;
end intrinsic;


/******************************** InertiaGroup ***************************************
 *
 * OK ring of ints in K, h map provided by AutomorphismGroup or ArtinRepresentations, p a rational prime
 *
 * Computes the inertia group of a chosen prime \frp above p.
 *
 *******************************************************************************************/

intrinsic InertiaGroup(OK :: RngOrd, h :: Map, p :: RngIntElt) -> GrpPerm, Map
{}
    local G, P, I, g, vals, is_in_I, i; 

    G := Domain(h);

    P := Factorization(p*OK)[1,1];
    I := [];
    for g in G do
        vals := [Valuation( h(g)(b) - b, P ) : b in Basis(OK)];
        is_in_I := true;
        i := 1;
        while i le #vals and is_in_I do
            if vals[i] lt 1 then
                is_in_I := false;
            end if;
            i :=i+1;
        end while;
        if is_in_I then
            Append(~I, g);
        end if;
    end for;
    return sub<G | I>;
end intrinsic;


/******************************** RamificationGroup ***************************************
 *
 * OK ring of ints in K, h map provided by AutomorphismGroup or ArtinRepresentation, 
 * p a rational prime, s a integer >= 0
 *
 * Computes the s-the ramification group G_s of a chosen prime \frp above p. (Of course, for
 * s = 0 this computes the inertia group.
 *
 *******************************************************************************************/

intrinsic RamificationGroup(OK :: RngOrd, h :: Map, p :: RngIntElt, s :: RngIntElt) -> GrpPerm, Map
{}
    local G, P, I, g, vals, is_in_I, i;

    G := Domain(h);

    P := Factorization(p*OK)[1,1];
    G_s := [];
    for g in G do
        vals := [Valuation( h(g)(b) - b, P ) : b in Basis(OK)];
        is_in_G_s := true;
        i := 1;
        while i le #vals and is_in_G_s do
            if vals[i] lt s+1 then
                is_in_G_s := false;
            end if;
            i :=i+1;
        end while;
        if is_in_G_s then
            Append(~G_s, g);
        end if;
    end for;
    return sub<G | G_s>;
end intrinsic;


/******************************** DecompositionGroup ***************************************
 *
 * OK ring of ints in K, h map provided by AutomorphismGroup or ArtinRepresentations, 
 * p a rational prime
 *
 * Computes the decomposition group of a chosen prime \frp above p. 
 *
 *******************************************************************************************/

intrinsic DecompositionGroup(OK :: RngOrd, h :: Map, p :: RngIntElt) -> GrpPerm, Map
{}
    local G, P, D, g, vals, is_in_D, i;

    G := Domain(h);

    P := Factorization(p*OK)[1,1];
    D := [];
    for g in G do
        Q := ideal<OK | [h(g)(b) : b in Basis(P)]>;
        if Q eq P then
            Append(~D, g);
        end if;
    end for;
    return sub<G | D>;
end intrinsic;


/******************************** FixedModule ***************************************
 *
 * V is a G-module and N a normal subgroup of G.
 *
 * Computes  V^N.
 *
 *******************************************************************************************/

intrinsic FixedModule(V :: ModGrp, N :: GrpPerm) -> ModGrp
{} 
    local ResV, VfixN, W;

    assert IsNormal(Group(V), N);

    ResV := Restriction(V, N);
    VfixN := Fix(ResV);
    W := sub<V | [V ! (ResV!b) : b in Basis(VfixN)]>;
    return W;
end intrinsic;



/******************************** ArtinConductor ***************************************
 *
 * V is a G-module, OK the ring of ints in K and h the map provided by AutomorphismGroup
 * or ArtinRepresentations.
 *
 * Computes the Artin conductor of the G-representation V.
 *
 *******************************************************************************************/

intrinsic ArtinConductor(V :: ModGrp,  OK :: RngOrd, h :: Map) -> RngIntElt
{}
    local F, p, Fp, D, ResV, H, g_0, s, W, i;

    F := [];
    for p in PrimeDivisors( Discriminant(OK) ) do
        Fp := [];
        D := DecompositionGroup(OK, h, p);
	ResV := Restriction(V, D);
	H := RamificationGroup(OK, h, p, 0);
	g_0 := #H;
	s := 0;
	while #H gt 1 do
	    W := FixedModule(ResV, H);
	    Append(~Fp, (#H / g_0 * (Dimension(V) - Dimension(W))) );
	    s := s+1;
            H := RamificationGroup(OK, h, p, s);
	end while;
	Append(~F, p^( Integers() ! (&+Fp)));
    end for;
    return &*F;
end intrinsic;


/******************************** Conductor ***************************************
 *
 * V is a G-module, E elliptic curve over Q, OK the ring of ints in K 
 * and h the map provided by AutomorphismGroup
 *
 * Computes the conductor of the system of compatible G-representations V \tensor T_p(E), where
 * T_p(E) is the usual Tate module of E/Q. We always assume that E and K have coprime 
 * ramification.
 *
 *******************************************************************************************/

intrinsic Conductor(V :: ModGrp, E :: CrvEll, OK :: RngOrd, h :: Map) -> RngIntElt
{}

    assert GCD(Discriminant(OK), Conductor(E)) eq 1;
    return ArtinConductor(V, OK, h)^2 * Conductor(E)^Dimension(V);
end intrinsic;


/******************************** OldInitLSeries ***************************************
 *
 * QG rational group ring, E elliptic curve over Q, OK the ring of ints in K (G = Gal(K/Q)),
 * h the map provided by AutomorphismGroup and prec the preciosion for L-series computations.
 *
 * Initializes the L-series [ [L(E, psi, s) : psi|chi] : chi \in Irr_Q(G) ].
 * No longer used. Dokchitser's implementation is much faster.
 *
 *******************************************************************************************/

intrinsic OldInitLSeries(QG :: AlgGrp, E :: CrvEll, OK :: RngOrd, h :: Map, prec :: RngIntElt) -> List
{}
    local L, orbit, M, f, gamma, Z;

    L := [* *];
    for orbit in QG`Mods do
        M := [* *];
	for V in orbit do
	    f := func< p, d | LocalPolynomial(V, E, OK, h, p) >;
	    gamma := [0 : i in [1..Dimension(V)]] cat [1 : i in [1..Dimension(V)]]; 
	    Z := LSeries(2, gamma, Conductor(V, E, OK, h), 0);
	    LSetPrecision(Z, prec);
	    LSetCoefficients(Z, f);
	    Append(~M, Z);
	end for;
	Append(~L, M);
    end for;
    return L;
end intrinsic;

/******************************** InitLSeries ***************************************
 *
 * QG rational group ring, E elliptic curve over Q, K the number field (G = Gal(K/Q)),
 * and prec the preciosion for L-series computations.
 *
 * Initializes the L-series [ [L(E,\bar\psi, s) : psi|chi] : chi \in Irr_Q(G) ].
 *
 * Attention: Note the \bar\psi !! This is exactly what is needed for ETNC.
 *
 *******************************************************************************************/
intrinsic InitLSeries(QG :: AlgGrp, E :: CrvEll, K :: FldNum, prec :: RngIntElt) -> List
{}
    local L, orbit, M, Z, A;

    L := [* *];
    for orbit in QG`Mods do
        M := [* *];
        for V in orbit do
            // A := K !! ElementToSequence( Character(Dual(V)) );
            // A := K !! Character(Dual(V) );
            A := K !! Character(V );

            Z := LSeries(E, A);
            LSetPrecision(Z, prec);
            Append(~M, Z);
        end for;
        Append(~L, M);
    end for;
    return L;
end intrinsic;


/******************************** InitLSeries2 ***************************************
 *
 * QG rational group ring, E elliptic curve over Q, K the number field (G = Gal(K/Q)),
 * and prec the preciosion for L-series computations.
 *
 * Initializes the L-series [ [L(E,psi, s) : psi|chi] : chi \in Irr_Q(G) ]. As InitLSeries,
 * but does not use the contragredient of psi. Note that for ETNC we have to work with
 * the contragredient !!!!
 *
 *******************************************************************************************/

intrinsic InitLSeries2(QG :: AlgGrp, E :: CrvEll, K :: FldNum, prec :: RngIntElt) -> List
{}
    local L, orbit, M, Z, A;

    L := [* *];
    for orbit in QG`X do
        M := [* *];
        for chi in orbit do
            A := K !! ElementToSequence( chi );
            Z := LSeries(E, A);
            LSetPrecision(Z, prec);
            Append(~M, Z);
        end for;
        Append(~L, M);
    end for;
    return L;
end intrinsic;

/******************************** Evaluate ***************************************
 *
 * Lvec is the list of L-series intitialized by InitLSeries. Computes the conjectural
 * analytic order and the leading term of the L-series.
 *
 *******************************************************************************************/

intrinsic Evaluate(Lvec :: List) -> SeqEnum, SeqEnum
{}
    local L, orbit, M, Z;

    L := [];
    orders := [];
    for i:=1 to #Lvec do
        orbit := Lvec[i];
        M := [];
        ord := [];
        for Z in orbit do
            val, o := LeadingTerm(Z);
            Append(~M, val);
            Append(~ord, o);
        end for;
        Append(~L, M);
        Append(~orders, ord);
    end for;
    return L, orders;
end intrinsic;

/******************************** AnalyticRankIsZero ***************************************
 *
 * Lvec is the list of L-series intitialized by InitLSeries. Returns true, if and only if
 * the analytic rank is zero for all twisted Hasse-Weil L-functions.
 *
 *******************************************************************************************/
intrinsic AnalyticRankIsZero(Lvec :: List) -> BoolElt
{}
    local i, val, prec;

    for i:=1 to #Lvec do
        val := Evaluate(Lvec[i, 1], 1);
	prec := Precision(val) - 3;  
	if Abs(val) lt 10^-prec then
	    return false;
	end if;
    end for;
    return true;
end intrinsic;


/******************************** LeadingTerm ***************************************
 *
 * Computes the conjectural order of vanishing and the leading coefficient in the Taylor expansion
 * of the L-series Z. Note that we are not able to decide whether a real complex number is actually
 * equal to 0.
 *
 *******************************************************************************************/

intrinsic LeadingTerm(Z :: LSer) -> FldComElt, RngIntElt
{}
    local k, val, prec;

    k := 0;
    val := Evaluate(Z, 1);
    prec := Floor( Precision(val) * 3/4);  /* this is very preliminary */
    while Abs(val) lt 10^-prec do
        print "In LeadingTerm k = ", k;
        k +:= 1;
	val :=  Evaluate(Z, 1 : Derivative:=k, Leading:=true);
    end while;
    return ComplexField() ! val/Factorial(k), k;
end intrinsic;


/******************************** LocalPolynomial ***************************************
 *
 * V G-module, E elliptic curve over Q, OK ring of ints, h as provided by AutomorphismGroup,
 * p a rational prime.
 *
 * Computes the local polynomial of the representation V \tensor T_p(E). The return value is
 * polynomial over the complex field.
 *
 *******************************************************************************************/

intrinsic LocalPolynomial(V :: ModGrp, E :: CrvEll, OK :: RngOrd, h :: Map, p :: RngIntElt) -> RngUPolElt
{}
    local ap, Lp, G, CX, alphas, betas, alpha, beta, Frob, W;

    assert GCD(Discriminant(OK), Conductor(E)) eq 1;
    G := Domain(h);
    RamifiedPrimes := PrimeDivisors( Discriminant(OK) );
    BadPrimes := PrimeDivisors( Conductor(E) );
    CX<x>:=PolynomialRing(ComplexField());

    if not(p in RamifiedPrimes) and not(p in BadPrimes)  then
        dim := 2 * Dimension(V);
        ap := FrobeniusTraceDirect(E,p);
        alphas := [<b[1]^-1, b[2]> : b in Roots( CX ! 1-ap*x + p*x^2 )];

        Frob := Frobenius(OK, h, p);
        d := Representation(V);
        W := d(Frob^-1);
        betas := [<ComplexField()!b[1], b[2]> : b in Eigenvalues(W)];
        Lp := &*[ (1 - a[1]*b[1]*x)^(b[2]*a[2]) : a in alphas, b in betas];
        // Lp := P ! [Round(a) : a in ElementToSequence(Lp)];
        return Lp;
    end if;

    if p in RamifiedPrimes  then
        ap := FrobeniusTraceDirect(E,p);
        alphas := [<b[1]^-1, b[2]> : b in Roots( CX ! 1-ap*x + p*x^2 )];
        D := DecompositionGroup(OK, h, p);
        I := InertiaGroup(OK, h, p);
        VfixI := FixedModule(Restriction(V, D), I);
        dim := 2 * Dimension(VfixI);
        if dim eq 0 then
            Lp := CX ! 1;
        else
            Frob := Frobenius(OK, h, p);
            d := Representation(VfixI);
            W := d(Frob^-1);
            betas := [<ComplexField()!b[1], b[2]> : b in Eigenvalues(W)];
            Lp := &*[ (1 - a[1]*b[1]*x)^(b[2]*a[2]) : a in alphas, b in betas];
            // Lp := P ! [Round(a) : a in ElementToSequence(Lp)];
            // Lp := Round( Integers() ! &*[ (p - a[1]*b[1])^(b[2]*a[2]) : a in alphas, b in betas] ) / p^dim;
        end if;
        return Lp;
    end if;

    if p in BadPrimes then
        if ReductionType(E, p) eq "Additive" then
            alphas := [];
            return CX ! 1;
        end if;
        if ReductionType(E, p) eq "Split multiplicative" then
            alphas := [1];
        end if;
        if ReductionType(E, p) eq "Nonsplit multiplicative" then
            alphas := [-1];
        end if;
        Frob := Frobenius(OK, h, p);
        d := Representation(V);
        dim := #alphas * Dimension(V);
        W := d(Frob^-1);
        betas := [<ComplexField()!b[1], b[2]> : b in Eigenvalues(W)];
        Lp :=  &*[ (1 - a*b[1]*x)^(b[2]) : a in alphas, b in betas];
        // Lp := P ! [Round(a) : a in ElementToSequence(Lp)];
        return Lp;
        // Lp := Round( Integers() ! &*[ (p - a*b[1])^(b[2]) : a in alphas, b in betas] ) / p^dim;
    end if;
end intrinsic;


/******************************** LocalEulerFactor ***************************************
 *
 * orbit is an orbit of conjugated absolutely irreducible representations, E elliptic curve over Q,
 * OK, h, p as usual.
 *
 * Evaluates the local polynomials at p^-1. This gives a SeqEnum of complex numbers which should
 * be conjugates in the character field K_i if orbit corresponds to the i-th component of QG.
 *
 * This must be changed so that one does not need any rounding process in the routine
 * LocalEulerFactors which builds on this function !!!
 *
 *******************************************************************************************/

intrinsic LocalEulerFactor(orbit :: List, E :: CrvEll, OK :: RngOrd, h :: Map, p :: RngIntElt) -> SeqEnum, RngIntElt
{}
    local V, maxdim, Lps, val;

    maxdim := 0;
    Lps := [];
    for V in orbit do
        f := LocalPolynomial(V, E, OK, h, p);
	val := Evaluate(f, p^-1);
	maxdim := Max(Degree(f), maxdim);
        Append(~Lps, val );
    end for;

    return Lps, p^maxdim;
end intrinsic;

/******************************** DualLocalEulerFactor ***************************************
 *
 * orbit is an orbit of conjugated absolutely irreducible representations, E elliptic curve over Q,
 * OK, h, p as usual.
 *
 * Evaluates the local polynomials at p^-1. This gives a SeqEnum of complex numbers which should
 * be conjugates in the character field K_i if orbit corresponds to the i-th component of QG.
 *
 * This must be changed so that one does not need any rounding process in the routine
 * DualLocalEulerFactors which builds on this function !!!
 *
 *******************************************************************************************/

intrinsic DualLocalEulerFactor(orbit :: List, E :: CrvEll, OK :: RngOrd, h :: Map, p :: RngIntElt) -> SeqEnum, RngIntElt
{}
    local V, maxdim, Lps, val;

    maxdim := 0;
    Lps := [];
    for V in orbit do
        f := LocalPolynomial(Dual(V), E, OK, h, p);
        val := Evaluate(f, p^-1);
        maxdim := Max(Degree(f), maxdim);
        Append(~Lps, val );
    end for;

    return Lps, p^maxdim;
end intrinsic;


/******************************** LocalEulerFactors ***************************************
 *
 * QG  is a group algebra with QG`Mods initialized, E elliptic curve over Q,
 * OK, h as usual, q a rational prime.
 *
 * Computes the value of the local polynomial at q evaluated at q^-1 as an element of the 
 * character field K_i, i = 1, ...,r.
 *
 * Should be changed so that we need no rounding process !!! For large primes q there are
 * actually problems with rounding.
 *
 *******************************************************************************************/

intrinsic LocalEulerFactors(QG :: AlgGrp, E :: CrvEll, OK :: RngOrd, h :: Map, q :: RngIntElt) -> List
{}

    local Z, ds, orbit, f, d, prec, RatZ;

    Z := [];
    ds := [];
    for orbit in QG`Mods do
        f, d := LocalEulerFactor(orbit, E, OK, h, q);
        Append(~Z,  f );
        Append(~ds, d);
    end for;
    prec := Round( Precision(Z[1][1]) / 2 );
    RatZ := MakeRational(Z, QG, LCM(ds), prec);
    return RatZ;
end intrinsic;


/******************************** DualLocalEulerFactors ***************************************
 *
 * QG  is a group algebra with QG`Mods initialized, E elliptic curve over Q,
 * OK, h as usual, q a rational prime.
 *
 * Computes the value of the local polynomial at q evaluated at q^-1 as an element of the
 * character field K_i, i = 1, ...,r.
 *
 * Should be changed so that we need no rounding process !!! For large primes q there are
 * actually problems with rounding.
 *
 *******************************************************************************************/

intrinsic DualLocalEulerFactors(QG :: AlgGrp, E :: CrvEll, OK :: RngOrd, h :: Map, q :: RngIntElt) -> List
{}

    local Z, ds, orbit, f, d, prec, RatZ;

    Z := [];
    ds := [];
    for orbit in QG`Mods do
        f, d := DualLocalEulerFactor(orbit, E, OK, h, q);
        Append(~Z,  f );
        Append(~ds, d);
    end for;
    prec := Round( Precision(Z[1][1]) / 2 );
    RatZ := MakeRational(Z, QG, LCM(ds), prec);
    return RatZ;
end intrinsic;


/******************************** Resolvents ***************************************
 *
 * QG  is a group algebra with QG`Mods initialized, 
 * h as usual, delta is a normal basis element for K/Q.
 *
 * Computes the resolvents of delta.  We return a List
 *             [ [ resolvent of delta wrt rep : rep in orbit] : orbits in QG`Rep ]
 *
 *******************************************************************************************/

intrinsic Resolvents(QG :: AlgGrp, h :: Map, delta :: FldNumElt) -> List
{}
    local G, comp, chi, T, d, c, S, k, t, s, i, j, res, R;

    G := Group(QG);
    R := [* *];

    for orbits in QG`Reps do
        M := [* *];
        for rep in orbits do
            T := [rep(g)  : g in G];
            d := Integers() ! Trace(rep(Id(G)));
            c := [Conjugates( h(g^-1)(delta) )[1] : g in G];
            // c := [Conjugates( h(g)(delta) )[1] : g in G]; /* Just for a test */

            S := [];
            for k:=1 to #T do
                t := T[k];
                s := ZeroMatrix(ComplexField(), d,d);

                for i:=1 to Nrows(t) do
                    for j:=1 to Ncols(t) do
                        // s[i,j] := (Rationals()!t[i,j])*c[k];
                        s[i,j] := Conjugates(t[i,j])[1] * c[k];
                    end for;
                end for;
                Append(~S, s);
            end for;

            // res := Conjugates( Determinant( &+S ) )[1];
            res := Determinant( &+S );
            Append(~M, res);
        end for;
        Append(~R, M);
    end for;
    return R;
end intrinsic;


/******************************** lCyclicRegulator ***************************************
 *
 * QG  is a group algebra with QG`Mods initialized,
 * h as usual, P is a point such that E(K)_p \simeq \Zl Q_1 \oplus \Zl Q_{n_0} \oplus \ZlG P.
 * Here K/Q is cyclic of order l.
 *
 * This function is used to compute the regulator at non-trivial characters in Wuthrich type examples.
 * For further details see the paper etncIV.
 *
 *******************************************************************************************/

intrinsic lCyclicRegulator(QG :: AlgGrp, h :: Map, P :: PtEll, prec :: RngIntElt) -> List
{}
    local G, R, orbits, M, chi, Rchi, g;

    G := Group(QG);
    K := Parent( P[1] );
    R := [ ];

    for orbits in QG`X do
        M := [];
        for chi in orbits do
            Rchi := &+[HeightPairing(P, GAction(P, g, h), K) * Conjugates(chi(g^-1))[1] : g in G];
            Append(~M, Rchi);
        end for;
        Append(~R, M);
    end for;
    // R[1,1] := 1;
    return R;
end intrinsic;


/******************************** lSquaredCyclicRegulator ***************************************
 *
 * QG  is a group algebra with QG`Mods initialized,
 * h as usual. 
 * We assume that K/Q is cyclic of order l^2 and that E(K)_p \simeq \Zl^{n_0} \oplus \Zl[G/H].
 * Here H is the subgroup of order l and we write F for the fixed field of H.
 * P in E(K) is a point which generates the \Zl[G/H]-component. So actually P in E(F).
 * sig = [u_0, u_1, u_2], and we assume that u_0 in {0,1}, u_1 = 1 and  u_2 = 0.
 *
 * This function is used to compute the regulator at non-trivial characters in lSquared type examples.
 * For further details see the paper etncIV.
 *
 *******************************************************************************************/

intrinsic lSquaredCyclicRegulator(QG :: AlgGrp, h :: Map, P :: SeqEnum, sig :: SeqEnum, prec :: RngIntElt) -> List
{}
    local G, R, orbits, M, chi, Rchi, g;

    G := Group(QG);
    K := Parent( P[1,1] );
    l := Integers() ! Sqrt(Degree(K));
    H := [S`subgroup : S in Subgroups(G) | S`order eq l][1];
    R := [ ];

    GmodH, pi := quo<G | H>;
    GmodH := [(pi^-1)(g) : g in GmodH];

    if sig[1] eq 0 and sig[2] eq 1 and sig[3] eq 0 then
        P11 := P[1];
    end if;
    if sig[1] eq 1 and sig[2] eq 1 and sig[3] eq 0 then
       P01 := P[1];     P11 := P[2];
    end if;

    for orbits in QG`X do
        M := [];
        if sig[1] eq 0 and sig[2] eq 1 and sig[3] eq 0 then
            for chi in orbits do
                Rchi := &+[1/l * HeightPairing(P11, GAction(P11, g, h), K) * Conjugates(chi(g^-1))[1] : g in GmodH];
                Append(~M, Rchi);
            end for;
        end if;
        if sig[1] eq 1 and sig[2] eq 1 and sig[3] eq 0 then
            if Order( orbits[1] ) eq 1 then
                RegMat := Matrix(RealField(), 2, 2, [1,2,3,4]);
		RegMat[1,1] := 1/l^2 * HeightPairing(P01, P01 , K);
		RegMat[1,2] := 1/l^2 * HeightPairing(P01, P11 , K);
		RegMat[2,1] := &+[1/l * HeightPairing(GAction(P11, g, h), P01, K) : g in GmodH];
		RegMat[2,2] := &+[1/l * HeightPairing(GAction(P11, g, h), P11, K)  : g in GmodH];
		Append(~M, Determinant(RegMat));
            else 
                for chi in orbits do
                    Rchi := &+[1/l * HeightPairing(P11, GAction(P11, g, h), K) * Conjugates(chi(g^-1))[1] : g in GmodH];
                    Append(~M, Rchi);
                end for;
            end if;
        end if;
        Append(~R, M);
    end for;
    for i:=1 to l*(l-1) do 
        R[3,i] := 1;
    end for;

    return R;
end intrinsic;


/******************************** ComputeZ ***************************************
 *
 * QG  is a group algebra, E elliptic curve over Q, Lvals is a list of complex L-values, more
 * precisely,
 * 		Lvals = [ [L(E, rep, 0) : rep in orbit] : orbit in QG`Reps
 * R is the list of resolvent values ordered as in Lvals,
 * Ovals is the list of period values, Rvals is the list of Regulators.
 * Note that in general we do not know how to compute E(K), hence in most cases where
 * E(K) is expected to be infinite we are not able to compute Rvals.
 *
 *
 * All values are complex numbers. This function computes for each absolutely irreducible character chi
 * the ratio
 *
 *                     L(E,chi,0) * (delta|chi) / ( Omega_+^{d^+(chi)} Omega_+^{d^-(chi)} )
 *
 * The elements in each component are expected to be conjugates in the correspoonding 
 * field K_i.
 *
 *******************************************************************************************/

intrinsic ComputeZ(QG :: AlgGrp, E :: CrvEll, Lvals :: SeqEnum, R:: List, Ovals :: List, Rvals :: SeqEnum) -> List
{}
    local comp, Z, z, i;

    Z := [];
    for comp:=1 to #Lvals do
        z := [];
	for i:=1 to #Lvals[comp] do
	    Append(~z, Lvals[comp][i]*R[comp][i]*Ovals[comp][i] / Rvals[comp][i]);
	end for;
	Append(~Z, z);
    end for;
    return  Z;
end intrinsic;
	    
        
/******************************** RealImagPeriod ***************************************
 *
 * E elliptic curve over Q
 *
 * Computes the real and purely imaginary period of E/Q. The definition of the purely imaginary period
 * is not completely clear to me.
 *
 * Attention: This has to be clarified as soon as we compute examples which are not totally real !!!!
 *
 *******************************************************************************************/
        
intrinsic RealImagPeriod(E :: CrvEll) -> SeqEnum
{}
    local om, k, prec, i;

    om := Periods(E);
    if  Sign( Discriminant(E) ) eq 1 then
        om[1] := 2*om[1];
        return om;
    else
        /* Attention: definition of the Omega_- is not clear. See also the function OldRealImagPeriod */
        return [om[1], -om[1] + 2*om[2]];
    end if;
end intrinsic;


intrinsic OldRealImagPeriod(E :: CrvEll) -> SeqEnum
{}
    local om, k, prec, i;

    om := Periods(E);
    k := 2*Real(om[2]) /  Real(om[1]) ;  /* The definiton of the imaginary period is not clear
                                            to me. The real period is probably correct, the
                                            imaginary period may well be wrong!!!!!!!!!!
                                         */
    prec := Precision(Parent(k))-1;
    assert Abs( k - Round(k) ) lt 10^(-prec);
    return [om[1], om[2]-k/2*om[1] ];
end intrinsic;


/******************************** CompOvals ***************************************
 *
 * QG the groupalgebra, E elliptic curve over Q, OK, h as usual.
 *
 * Computes for each absolutely irreducible representation V the value 
 * Omaga_+^{-dplus} * Omega_-^{-dminus}, where dplus (resp. dminus) is the dimension
 * (resp. codim) of V^{c=1}, c complex conjugation.
 * 
 * We return a list [ [value for rep : rep in orbit] : orbit in QG`Mods ]
 *
 * Attention: The imaginary periods may well be incorrect. This has to be clarified before
 * used with fields which are not totally real !!!!
 *
 *******************************************************************************************/

intrinsic CompOvals(QG :: AlgGrp, E :: CrvEll, OK :: RngOrd, h :: Map) -> List
{}
    local orbit, O, o, V, C, Vfix, dplus, dminus;

    O := [* *];
    om := RealImagPeriod(E);
    for orbit  in QG`Mods do
        o := [* *];
	for V in orbit do
	    C := sub< Group(QG) | ComplexConjugation(OK, h) >;
	    Vfix := FixedModule(Restriction(V, C), C);
	    dplus := Dimension(Vfix); dminus := Dimension(V) - Dimension(Vfix);
	    Append(~o, om[1]^(-dplus) * om[2]^(-dminus) );
	end for;
	Append(~O, o);
    end for;
    return O;
end intrinsic;


/******************************** IsSemistable ***************************************
 *
 * E elliptic curve over Q
 *
 * Checks whether E/Q is semistable.
 *
 *******************************************************************************************/

intrinsic IsSemistable(E :: CrvEll) -> BoolElt
{}
    local p, t;

    for p in PrimeDivisors(Conductor(E)) do
        t := ReductionType(E, p);
	if not( t eq "Split multiplicative" or t eq "Nonsplit multiplicative") then
	    return false;
	end if;
    end for;
    return true;
end intrinsic;


/******************************** IsSplitMultiplicative ***************************************
 *
 * E elliptic curve over Q
 *
 * Checks whether E/Q is split multiplicative.
 *
 *******************************************************************************************/

intrinsic IsSplitMultiplicative(E :: CrvEll) -> BoolElt
{}
    local p, t;

    for p in PrimeDivisors(Conductor(E)) do
        t := ReductionType(E, p);
        if not( t eq "Split multiplicative") then
            return false;
        end if;
    end for;
    return true;
end intrinsic;


/******************************** TateShafarevicGroup ***************************************
 *
 * E elliptic curve over Q, QG rational group ring, Lvals precomputed vector of leading terms,
 * Rvals precomputed vector of equivariant regulators, OK ring of integers.
 *
 * Uses the BSD for E/K to compute a conjectural value for the order of Sha(E/K). 
 *
 * Attention: One has to be carefull with the equivariant regulators. One has to make sure
 * that the product over the regulators according to Artin formalism is the regulator of $E(K)$.
 * In practice, this is difficult.
 *
 *******************************************************************************************/

intrinsic TateShafarevicGroup(E :: CrvEll, QG :: AlgGrp, Lvals :: SeqEnum, Rvals :: SeqEnum, OK :: RngOrd) -> RngIntElt, FldReElt
{}
    local K, EK, LK, om, r1, r2, c, w, d, ts, prec;

    if assigned E`ts then 
        print "E`ts already assigned";
        return E`ts; 
    end if;

    // assert IsSemistable(E);
    prec := Precision(RealField());

    K := NumberField(OK);
    EK := BaseChange(E, K);

    /* Use previously computed L-values and Artin formalism to compute L(E/K, 1) */
    LK := 1;
    for i:=1 to #Lvals do
        for j:=1 to #Lvals[i] do
            // LK *:= Lvals[i][j]^Degree(QG`X[i][j]) / Rvals[i][j];
	    LK *:= (Lvals[i][j] / Rvals[i][j])^Degree(QG`X[i][j]);
        end for;
    end for;
    om := RealImagPeriod(E);
    r1, r2 := Signature(OK);
    c := &*TamagawaNumbers(E, OK);
    w := #TorsionSubgroup(EK);
    d := Sqrt(Abs(Discriminant(OK)));

    // ts := (LK * d * w^2) / (om[1]^(r1+r2) * Abs(2*om[2])^r2  * c);
    ts := (LK * d * w^2) / (om[1]^(r1+r2) * Abs(om[2])^r2  * c);

    // print "Approximate conjectural real value for the order of Sha = ", ts;
    // readi OrderOfSha, "Input an order";
    // return OrderOfSha;

    OrderOfSha := Runde(ts, Round(prec / 2) );
    // assert IsSquare(OrderOfSha);
    return OrderOfSha, ts;
end intrinsic;

/******************************** AnalyticTateShafarevicGroup ***************************************
 *
 * E elliptic curve over Q, QG rational group ring, Lvals precomputed vector of leading terms,
 * Rvals precomputed vector of equivariant regulators, OK ring of integers.
 *
 * Uses the BSD for E/K to compute a conjectural value for the order of Sha(E/K).
 * essentially the same as TateShafarevicGroup but does not try to round the complex value.
 *
 * Attention: One has to be carefull with the equivariant regulators. One has to make sure
 * that the product over the regulators according to Artin formalism is the regulator of $E(K)$.
 *
 *******************************************************************************************/


intrinsic AnalyticTateShafarevicGroup(E :: CrvEll, QG :: AlgGrp, Lvals :: SeqEnum, Rvals :: SeqEnum, OK :: RngOrd) -> FldReElt
{}
    local K, EK, LK, om, r1, r2, c, w, d, ts, prec;

    if assigned E`ts then
        print "E`ts already assigned";
        return E`ts;
    end if;

    // assert IsSemistable(E);
    prec := Precision(RealField());

    K := NumberField(OK);
    EK := BaseChange(E, K);

    /* Use previously computed L-values and Artin formalism to compute L(E/K, 1) */
    LK := 1;
    for i:=1 to #Lvals do
        for j:=1 to #Lvals[i] do
            // LK *:= Lvals[i][j]^Degree(QG`X[i][j]) / Rvals[i][j];
            LK *:= (Lvals[i][j] / Rvals[i][j])^Degree(QG`X[i][j]);
        end for;
    end for;
    om := RealImagPeriod(E);
    r1, r2 := Signature(OK);
    c := &*TamagawaNumbers(E, OK);
    w := #TorsionSubgroup(EK);
    d := Sqrt(Abs(Discriminant(OK)));

    // ts := (LK * d * w^2) / (om[1]^(r1+r2) * Abs(2*om[2])^r2  * c);
    ts := (LK * d * w^2) / (om[1]^(r1+r2) * Abs(om[2])^r2  * c);

    return ts;
end intrinsic;

/******************************** TamagawaNumbers ***************************************
 *
 * E elliptic curve over Q, OK ring of ints of K.
 *
 * We assume that E/Q is semistable.
 * Computes the list of Tamagawa numbers for all bad primes. We use Silverman, Th. 6.1. Note that
 * if E/Q is semistable, then E/K is also semistable.
 * One should apply 6.1 directly, without using the magma function TamagawaNumber. If the
 * reduction is split multiplicative, this is easy by 6.1
 *
 * See the next function NewTamagawaNumbers.
 *
 *******************************************************************************************/

intrinsic TamagawaNumbers(E :: CrvEll, OK :: RngOrd) -> SeqEnum
{}
    local c, p, cp, fac, f;

    // assert IsSemistable(E);

    c := [];
    for p in PrimeDivisors( Conductor(E) ) do
        cp := TamagawaNumber(E, p);
	fac := Factorization(p*OK);
	c cat:= [f[2]*cp : f in fac];
    end for;
    return c;
end intrinsic;

/******************************** NewTamagawaNumbers ***************************************
 *
 * E elliptic curve over Q, OK ring of ints of K.
 *
 * We assume that E/Q is semistable.
 * We directly apply Silverman Th.6.1 without using the magma function TamagawaNumber. If the
 * reduction is split multiplicative, this is easy by 6.1
 *
 * See the next function NewTamagawaNumbers.
 *
 *******************************************************************************************/

intrinsic NewTamagawaNumbers(E :: CrvEll, OK :: RngOrd) -> SeqEnum
{}
    local c, p, cp, fac, f;

    assert IsSemistable(E);

    c := [];
    j := jInvariant(E);
    for p in PrimeDivisors( Conductor(E) ) do
        fac := Factorization(p*OK);
        if ReductionType(E, p) eq "Split multiplicative" then
            c cat:= [-Valuation(j*OK, f[1]) : f in fac];
        elif IsDivisibleBy( Valuation(j*OK, fac[1,1]), 2) then
            c cat:= [2 : f in fac];
        else
           c cat:= [1 : f in fac];
        end if;
    end for;
    return c;
end intrinsic;


/******************************** TamagawaNumber ***************************************
 *
 * E elliptic curve over Q, P is a prime ideal in OK
 *
 * We assume that E/Q is semistable. (It would be enough to assume semistability of E at P).
 * Again we use Silverman, Th.6.1, not directly,  but for E/Q.
 * This should be changed and 6.1 should be applied directly. See the next function NewTamagawaNumber.
 *
 *******************************************************************************************/

intrinsic TamagawaNumber(E :: CrvEll, P :: RngOrdIdl) -> RngIntElt
{}
    local p;

    assert GCD(Conductor(E), Discriminant(Order(P))) eq 1;
    assert IsSemistable(E);

    p := Characteristic( quo< Order(P) | P > );
    return TamagawaNumber(E, p);
end intrinsic;

/******************************** NewTamagawaNumber ***************************************
 *
 * E elliptic curve over Q, P is a prime ideal in OK
 *
 * We assume that E/Q is semistable. (It would be enough to assume semistability of E at P).
 * We apply Th.6.1 directly.
 *
 *******************************************************************************************/

intrinsic NewTamagawaNumber(E :: CrvEll, P :: RngOrdIdl) -> RngIntElt
{}
    local j, p, type;


    j := jInvariant(E);
    p := Characteristic( quo< Order(P) | P > );
    type := ReductionType(E, p);
    assert type eq "Split multiplicative" or type eq "Nonsplit multiplicative" or type eq "Good";

    if type eq "Good" then
       return 1;
    elif type eq "Split multiplicative" then
       return -Valuation(j*Order(P), P);
    elif IsDivisibleBy( Valuation(j*Order(P), P), 2) then
       return 2;
    else
       return 1;
    end if;
end intrinsic;

/******************************** EKtorRep ***************************************
 *
 * EKtor is the torsion subgroup of E(K), j the map EKtor ->E(K) given by the magma function
 * TorsionSubgroup, h the map given by AutomorphismGroup or ArtinRepresentations.
 *
 * Computes  the matrices A(g), g in G, which give the action of G in E(K)_tors with respect to
 * Generators(EKtor).
 *
 *******************************************************************************************/

intrinsic EKtorRep(EKtor :: GrpAb, j :: Map, h :: Map) -> SeqEnum
{}
    local A, G, gens, g, M;

    A := [];
    G := Domain(h);
    gens := [EKtor.i : i in [1..#Generators(EKtor)]];
    for g in G do
        M := [ EKtorLog(GAction(g, j(P), h), EKtor, j) : P in gens ];
	Append(~A, Transpose( Matrix(M) ));
    end for;
    return A;
end intrinsic;


/******************************** EKtorDualRep ***************************************
 *
 * EKtor is the torsion subgroup of E(K), j the map EKtor ->E(K) given by the magma function
 * TorsionSubgroup, h the map given by AutomorphismGroup or ArtinRepresentations.
 *
 * Computes  the matrices A(g), g in G, which give the action of G in E(K)_tors with respect to
 * Generators(EKtor).
 *
 *******************************************************************************************/

intrinsic EKtorDualRep(EKtor :: GrpAb, j :: Map, h :: Map) -> SeqEnum
{}
    local A, G, gens, g, M;

    A := [];
    G := Domain(h);
    gens := [EKtor.i : i in [1..#Generators(EKtor)]];
    for g in G do
        M := [ EKtorLog(GAction(g^-1, j(P), h), EKtor, j) : P in gens ];
        Append(~A, ( Matrix(M) ));
    end for;
    return A;
end intrinsic;


/******************************** ZGAction ***************************************
 *
 * lambda is a element of ZG,  P a point in E(K) and h the map given by 
 * AutomorphismGroup or ArtinRepresentations.
 *
 * Computes  P^lambda.
 *
 *******************************************************************************************/

intrinsic ZGAction(lambda :: AlgGrpElt, P :: PtEll, h :: Map) -> PtEll
{}
    local g;

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

end intrinsic;

/******************************** GAction ***************************************
 *
 * g is a element of G = Gal(K/Q),  P a point in E(K) and h the map given by AutomorphismGroup
 *
 * Computes  P^g.
 *
 *******************************************************************************************/

intrinsic GAction(g :: GrpPermElt, P :: PtEll, h :: Map) -> PtEll
{}
    return Parent(P) ! [h(g)(x) : x in ElementToSequence(P)];
end intrinsic;




/******************************** EKtorLog ***************************************
 *
 * P is a point in E(K)_tors, EKtor the abstract abelian group representing E(K)_tors,
 * j the map given by TorsionSubgroup.
 *
 * Computes  the coefficients of P with respect to the generators of EKtors.
 *
 *******************************************************************************************/

intrinsic EKtorLog(P :: PtEll, EKtor :: GrpAb, j :: Map) -> SeqEnum
{}
    local q;

    for q in EKtor do
        if j(q) eq P then 
	    return ElementToSequence(q);
	end if;
    end for;
end intrinsic;

/******************************** HardPrimes ***************************************
 *
 * E elliptic curve over Q, QG group algebra, Lvals the list of L-values ordered as
 * [ L(E, chi, 1) : chi in orbit] : orbit in QG`Mods ], Rvals the list of
 * equivariant regulators, OK ring of ints.
 *
 * Computes  the set of primes which either divide #E(K)_tors, #G, #Sha(E/K) or one of
 * the Tamagawa numbers for E/K.  We have a proposition which
 * allows to check ETNC (numerically) for all primes outside the set of hard primes.
 * 
 * Attention: This is conjectural because we use BSD for E/K in order to compute
 * a conjectural value for the order of Sha(E/K).
 *
 *******************************************************************************************/

intrinsic HardPrimes(E :: CrvEll, QG :: AlgGrp, Lvals :: SeqEnum, Rvals :: SeqEnum, OK :: RngOrd) -> SetEnum
{}
    local K, EK, b, c;

    K := NumberField(OK);
    EK := BaseChange(E, K);

    b := PrimeDivisors( #TorsionSubgroup(EK) );
    b cat:= [2];
    b cat:= PrimeDivisors( TateShafarevicGroup(E, QG, Lvals, Rvals, OK) );
    b cat:= PrimeDivisors( #Group(QG) );
    c := Set( TamagawaNumbers(E, OK) );
    for x in c do
        b cat:= PrimeDivisors(x);
    end for;
    return Set(b);
end intrinsic;


/******************************** ComputeS ***************************************
 *
 * E elliptic curve over Q, OK ring of ints.
 *
 * Computes  the set of bad primes for E/K.
 *
 *******************************************************************************************/

intrinsic ComputeS(E :: CrvEll, OK :: RngOrd) -> SetEnum
{}
    local b;

    b := PrimeDivisors( Conductor(E) );
    b cat:= PrimeDivisors( Discriminant(OK) );
    return Set(b);
end intrinsic;


    
    
/******************************** MakeRational ***************************************
 *
 * Z is a list of complex numbers ordered in the following way:
 *             Z = [ [z_\chi : chi in orbit] : orbit in QG`MOds ]
 * QG group algebra, d an integers such that all d*z_\chi should be integral, prec an integer such
 * that all rounding is done with an error less then 10^-prec.
 *
 * Computes z_\chi as an element of the character field K_i, i = 1, ..., r, and the maximal rounding error.
 * 
 * This version is still needed for the computation of the Euler factors. 
 * For other purposes use the function MakeRational below which uses the MAGMA functions
 * BestApproximation();
 *
 *******************************************************************************************/

intrinsic MakeRational(Z :: SeqEnum, QG :: AlgGrp, d :: RngIntElt, prec :: RngIntElt) -> List, FldReElt
{}
    local C, CX, ZX, RatZ, i, Ki, KiX, f, g, roots, found, j, conj, k, v, fehler;

    C := Parent(Z[1][1]);

    CX<T> := PolynomialRing(C);
    ZX := PolynomialRing( Integers());

    fehler := -1;
    RatZ := [**];
    for i:=1 to #Z do
        // print "Rationalize i = ", i;
        Ki := QG`H[i][1];              /* the character field */
	KiX := PolynomialRing(Ki);
	n := #Z[i];
	f := d^n * &*[ T - z : z in Z[i] ];
	g := KiX![Runde(c, prec)/d^n : c in ElementToSequence(f)];
	roots := Roots(g); /* This may take a long time */
        v := [];
        for k:=1 to #Z[i] do
	    found := false;
	    j:=1;
	    while not found and j le #roots do
	        conj := Conjugates(roots[j][1])[1];
	        if Abs(conj - Z[i][k]) lt 10^-prec then
	            found := true;
                    fehler := Max(fehler, Abs(conj - Z[i][k]));
                    // print "Error = ", Abs(conj - Z[i][k]) ;
		    Append(~v, roots[j][1]);
	        end if;
	        j := j+1;
	    end while;
	    assert found;
        end for;
        Append(~RatZ, v);
    end for;

    return RatZ, fehler;
end intrinsic;

/******************************** MakeRational ***************************************
 *
 * Z is a list of complex numbers ordered in the following way:
 *             Z = [ [z_\chi : chi in orbit] : orbit in QG`MOds ]
 * QG group algebra, d an integers such that all d*z_\chi should be integral, prec an integer such
 * that all rounding is done with an error less then 10^-prec.
 *
 * Computes z_\chi as an element of the character field K_i, i = 1, ..., r, and the maximal rounding error.
 *
 * This version uses BestApproximation(c, N) ); in order to find the denominators.
 *
 *******************************************************************************************/

intrinsic MakeRational(Z :: SeqEnum, QG :: AlgGrp) -> List, FldReElt
{}
    local C, CX, ZX, RatZ, i, Ki, KiX, f, g, roots, found, j, conj, k, v, fehler;

    N := 2^5 * 3^3 * 5^2 * 7 * 11 * 13; /* maximal denominator, this has to be improved */
    C := Parent(Z[1][1]);
    prec := Round( 2*Precision(C) / 3 );

    CX<T> := PolynomialRing(C);
    ZX := PolynomialRing( Integers());

    fehler := -1;
    RatZ := [**];
    for i:=1 to #Z do
        // print "Rationalize i = ", i;
        Ki := QG`H[i][1];              /* the character field */
        KiX := PolynomialRing(Ki);
        n := #Z[i];
        f := &*[ T - z : z in Z[i] ];
        
        fehler:= Max( [Abs(Imaginary(c)): c in ElementToSequence(f)] );
        assert fehler lt 10^-prec;
        g := KiX![BestApproximation(Real(c), N) : c in ElementToSequence(f)];
        roots := Roots(g); /* This may take a long time */
        v := [];
        for k:=1 to #Z[i] do
            found := false;
            j:=1;
            while not found and j le #roots do
                conj := Conjugates(roots[j][1])[1];
                if Abs(conj - Z[i][k]) lt 10^-prec then
                    found := true;
                    fehler := Max(fehler, Abs(conj - Z[i][k]));
                    // print "Error = ", Abs(conj - Z[i][k]) ;
                    Append(~v, roots[j][1]);
                end if;
                j := j+1;
            end while;
            assert found;
        end for;
        Append(~RatZ, v);
    end for;

    assert fehler lt 10^-prec;
    return RatZ, fehler;
end intrinsic;

/******************************** Runde ***************************************
 *
 * z is a complex number which is supposed to be a rational integer. prec gives the precision.
 *
 * We round z. The error in both real and imaginary part must be < 10^-prec.
 *
 *******************************************************************************************/

intrinsic Runde(z :: FldComElt, prec :: RngIntElt) -> RngIntElt
{}
    local x, y, N ;

    x := Round(Real(z)); y := Round(Imaginary(z));
    N := prec;  
    if not ( Abs(Imaginary(z)) lt 10^-N and Abs(Real(z) - x) lt 10^-N ) then
        print "Cannot round !!!    z = ", z;
        print "Error = ", Max(Abs(Imaginary(z)), Abs(Real(z) - x));
        readi x, "Input a conjectural integer"; 
    end if;
    // assert Abs(Imaginary(z)) lt 10^-N and Abs(Real(z) - x) lt 10^-N;
    return Integers() ! x;
end intrinsic;
	
	
/******************************** CheckEasyPrimes ***************************************
 *
 * RatZ is in each component of QG given by Lval*resolvent / Rval*(Omega_+^dplus * Omega_-^dminus)
 * as an element of the character field K_i, i = 1, ..., r. HardPrimes is the set of hard primes.
 *
 * We check that ETNC_p is true for all primes not in HardPrimes.This is based on our
 * Proposition.
 *
 *******************************************************************************************/

intrinsic CheckEasyPrimes(RatZ :: List, HardPrimes :: SetEnum) -> BoolElt
{}
    local i, n, S, p;

    for i:=1 to #RatZ do
           n := Norm(RatZ[i][1]);
           S := Set (PrimeDivisors( Numerator(n) ) cat PrimeDivisors( Denominator(n) ) );
           for p in S do
               if not(p in HardPrimes) then
	           return false;
	       end if;
	   end for;
    end for;
    return true;
end intrinsic;
       

/******************************** CheckHardPrime ***************************************
 *
 * QG group algebra, E elliptic curve over Q, S set of bad primes for E/K, l a `hard´ prime.
 *
 * At the moment we return two boolean values is_computable and ETNCIsValid;
 *
 * Because we always work with the cohomology modules (instead of the perfect complexes which
 * are used to state ETNC) we must assume that all cohomology modules are perfect. 
 * Unfortunately, even if we are in a regular situation where by
 * definition everything is perfect we may not be able to compute all Euler characteristics.
 *
 *******************************************************************************************/

intrinsic CheckHardPrime(QG :: AlgGrp, E :: CrvEll, S :: SetEnum, l :: RngIntElt) -> BoolElt, BoolElt
{}
    local Sl, eps_kl, EKtor, j, bool, eps_EKtor, eps_Sha, eps_Ens, q, eps, Prod_eps_Ens,
          eps_Tam, Prod_eps_Tam, euler, Prod_euler, xi, Krel,
	  is_regular, is_computable, ETNCIsValid;

    print "In CheckHardPrime l = ", l;

    if (l in S or l eq 2) and IsDivisibleBy(#Group(QG), l) then
        print "Cannot deal with l because (l in S or l eq 2) and IsDivisibleBy(#Group(QG), l)";
	return false, false;
    end if;
    if not(l in S or l eq 2) then
        tam := NewTamagawaNumbers(E, E`OK);
        if IsDivisibleBy( &*tam, l) then
            PrintList("Tamagawa numbers = ", tam);
            print "Cannot deal with l because l divides a Tamagawa number";
            return false, false;
        end if;
        for p in S do
	        Ip := InertiaGroup(E`OK, E`h, p);
		if IsDivisibleBy(#Ip, l) then
		    print "Cannot deal with l because l divides I_", p;
		    return false, false;
		end if;
	end for;
    end if;

    Sl := S join {l};
    xi := [* 1 : v in E`RatZ *];
    eta := [* 1 : v in E`RatZ *];

    /* Compute eps( E(K)_tors ) * eps( Pontrijagin dual of E(K)_tors */
    EKtor, j := TorsionSubgroup(E`EK);
    is_computable, eps_EKtor := Eps_EKtor(QG, EKtor, j, E`h, l);
    PrintList( "eps_EKtor = ", eps_EKtor );

    /* Compute eps( Sha(E/K) ) */
    is_c, eps_Sha := Eps_Sha(QG, E, l);
    is_computable := is_computable and is_c; 
    PrintList( "eps_Sha = ", eps_Sha );


    if l in S or l eq 2 then
        
        /* Compute eps( Tamagawa module at q ). This is only OK if the reduction is split multiplicative.
        */
        eps_Tam := [ ];
        for q in Sl do
            is_c, eps:= Eps_Tam(QG, E`OK, E, E`h, l, q);
            printf "eps_Tam_%o = ", q; PrintList("", eps);
            Append(~eps_Tam, eps);
            is_computable := is_computable and is_c;
        end for;
        Prod_eps_Tam := [* &*[v[i] : v in eps_Tam] : i in [1..#eps_Tam[1]] *];
        PrintList( "Prod_eps_Tam_", Prod_eps_Tam );

        /* For each q in Sl we compute eps( \bar E_{q, ns} (k_v) ). More precisely: Let v be a place
           above q, then consider the non-singular k_v-rational points of the reduction of E mod q as
           a \Zl[D]-module, where D denotes the decomposition group.
        */

        eps_Ens := [ ];
        for q in Sl do
            is_c, eps := Eps_Ens(QG, E`OK, E, E`h, l, q);
	    printf "Eps_Ens_%o = ", q; PrintList("", eps);
            Append(~eps_Ens, eps);
            is_computable := is_computable and is_c; 
        end for;
        Prod_eps_Ens := [* &*[v[i] : v in eps_Ens] : i in [1..#eps_Ens[1]] *];
	PrintList( "Prod_eps_Ens = ", Prod_eps_Ens );

        /* Compute eps( \oplus_{v|l} k_v ) */
        eps_kl := Eps_kl(QG, E`h, E`OK, l);
	PrintList( "eps_kl = ", eps_kl );


        /* Compute the Euler factors for all primes q in Sl. These are elemets of the character fields 
           K_i, i = 1, ..., r. */    
        euler := [ ];
        for q in Sl do
	    eul := DualLocalEulerFactors(QG, E, E`OK, E`h, q);
	    Rat_eul := [* v[1] : v in eul *];
            printf "EulerFactor_%o = ", q; PrintList("", Rat_eul);
            Append(~euler, eul);
        end for;
        Prod_euler := [* &*[v[i][1] : v in euler] : i in [1..#euler[1]] *];
	PrintList( "Prod_Factors = ", Prod_euler );

        for i:= 1 to #E`RatZ do
            xi[i] := eps_kl[i]^-1 * Prod_eps_Ens[i] * Prod_euler[i]^-1 * Prod_eps_Tam[i];
        end for;
    end if;

    for i:= 1 to #E`RatZ do
            xi[i] := xi[i] * eps_EKtor[i]^-1 * eps_Sha[i];
            eta[i] := xi[i]^-1 * E`RatZ[i][1];
    end for; 
    printf "xi_%o = ", l; PrintList("",  xi);
    printf "u * xi_%o^-1 = ", l; PrintList("", eta);

    Krel := RelativeGroup(QG, l);
    log := K0RelLog(Krel, eta);
    // print "******************l = ", l, " completely checked *************************";
    print "log = ",log;

    ETNCIsValid := IdealPartIsTrivial(log) and (log[2] eq Id(Parent(log[2])));
    return is_computable, ETNCIsValid;
end intrinsic;

	
    
/******************************** OldEps_kl ***************************************
 *
 * QG group algebra, h given by AutomorphismGroup or ArtinRepresentations, 
 * OK ring of ints, l rational prime.
 *
 * Computes eps( \prod_{v|l} k_v ) = \ind_{G_v}^G k_v. We use the ses
 *            0 -> < \Delta(D, I), l > -> Z_l[D] -> k_v -> 0
 * where D (resp. I) denotes the decomposition group (resp. inertia group) of the chosen
 * place v above l. 
 *
 * No longer used.
 *
 *******************************************************************************************/

intrinsic OldEps_kl(QG :: AlgGrp, h :: Map, OK :: RngOrd, l :: RngIntElt) -> List
{}
    local D, I, QD, rho, ZD, DGH, BasisDGH, BasisZD, theta, S, nr;

    D := DecompositionGroup(OK, h, l);
    I := InertiaGroup(OK, h, l);

    QD := InitGroupAlgebra(D);
    rho := RegularRep(QD);
    ZD := RegularModule(rho);
    DGH := DeltaGH(QD, I, l, rho);
    BasisDGH := LocalBasis(DGH, l);
    BasisZD := LocalBasis(ZD, l);
    theta := IdentityMatrix(Rationals(), #D);
    S := QGMatrix(QD, theta, DGH, BasisDGH, ZD, BasisZD);
    nr:= NewtonReducedNorm(S);

    E := CommonSplittingField(QG, QD);

    Tau := EmbeddingsToE(QD, E);
    beta := [**];
    for i:=1 to #QG`H do
        chi := Restriction(QG`H[i][2], D);
        v := [];
        for j:=1 to #nr do
            v cat:= [tau(nr[j])^ScalarProduct(D, tau, QD`H[j][2], QD`H[j][3], chi) : tau in Tau[j]];
        end for;
        Append(~beta, &*v);
    end for;
   
    Tau := EmbeddingsToE(QG, E);
    alpha := [**];
    for i:=1 to #beta do
        iota := Tau[i][1]; OKi := QG`H[i][4];
        /* For some groups there are PROBLEMS with coercion and domains of maps in the next line ??? */
        Append(~alpha, (iota^-1)(beta[i]) );
    end for;
    return alpha;
end intrinsic;



/******************************** Eps_kl ***************************************
 *
 * QG group algebra, h given by AutomorphismGroup or ArtinRepresentations,
 * OK ring of ints, l rational prime.
 *
 * Computes eps( \prod_{v|l} k_v ) = \ind_{G_v}^G k_v. We use the ses
 *            0 -> < \Delta(D, I), l > -> Z_l[D] -> k_v -> 0
 * where D (resp. I) denotes the decomposition group (resp. inertia group) of the chosen
 * place v above l.
 * This uses Proposition 5.5 of the paper.
 *
 *******************************************************************************************/

intrinsic Eps_kl(QG :: AlgGrp, h :: Map, OK :: RngOrd, l :: RngIntElt) -> List
{}
    local D, I, QD, rho, ZD, DGH, BasisDGH, BasisZD, theta, S, nr;

    D := DecompositionGroup(OK, h, l);
    I := InertiaGroup(OK, h, l);

    QD := InitGroupAlgebra(D);
    nr := [* H[1]!1 : H in QD`H *];
    for i:=1 to #QD`X do
        if I subset Kernel(QD`X[i,1]) then
            nr[i] := QD`H[i,1] ! l;
	end if;
    end for;
    return Induction(QG, QD, nr);
end intrinsic;


/******************************** Eps_EKtor ***************************************
 *
 * QG group algebra, EKtor is the torsion subgroup of E(K) as abstract abelian group,
 * j: EKtor -> E(K), h as usual, l rational prime.
 *
 * Computes eps( E(K)_{tors} ) * eps( E(K)^*_{tors} )
 *
 *******************************************************************************************/

intrinsic Eps_EKtor(QG :: AlgGrp, EKtor :: GrpAb, j :: Map, h :: Map, l :: RngIntElt) -> BoolElt, List
{}
    local is_regular, G, A, nr; 

    G := Domain(h);
    if IsDivisibleBy(#G, l) and IsDivisibleBy(#EKtor, l) then
        is_regular := false; 
	print "EKtor is not regular for l = ", l;
        return is_regular, [* H[1]!1 : H in QG`H *];
    end if;
    if not IsDivisibleBy(#EKtor, l) then
        is_regular := true; 
        return is_regular, [* H[1]!1 : H in QG`H *];
    else
        A := EKtorRep(EKtor, j, h);
	At := EKtorDualRep(EKtor, j, h);
        nr := ReducedNorm(EKtor, A, QG, l);
	nrt := ReducedNorm(EKtor, At, QG, l);
	nr_times_nrt := [* nr[i]*nrt[i] : i in [1..#nr] *];
	is_regular := true; 
        return is_regular, nr_times_nrt;
    end if;
end intrinsic;


/******************************** Eps_Sha ***************************************
 *
 * QG group algebra, E elliptic curve over Q, Lvals the L-values ordered as before, 
 * OK ring of ints, l rational prime.
 *
 * Computes eps( Sha(E/K) ). This is very weak. Essentially we cannot compute this contribution.
 * We use BSD for E/K to compute #Sha(E/K). If this is divisible by l, then we do not know how
 * to continue, even in the regular case when l does not divide #G.
 * 
 * In some examples the computational results conjecturally tell something about the
 * components of Sha for regular primes l.
 *
 *******************************************************************************************/

intrinsic Eps_Sha(QG :: AlgGrp, E :: CrvEll, l :: RngIntElt) -> BoolElt, List
{}
    local is_regular, ord;

    ord := TateShafarevicGroup(E, QG, E`Lvals, E`Rvals, E`OK);
    if IsDivisibleBy(#Group(QG), l) and IsDivisibleBy(ord, l) then
        print "Sha is not regular for l = ", l;
        is_regular := false; 
        return  is_regular, [* H[1]!1 : H in QG`H *];
    end if;
    if not IsDivisibleBy(ord, l) then
       is_regular := true; 
       return is_regular, [* H[1]!1 : H in QG`H *];
    else
       print "Cannot compute Eps_Sha because l = ", l, " divides #Sha = ", ord;
       return false, [* H[1]!1 : H in QG`H *];
    end if;
end intrinsic;

    
/******************************** Eps_Ens ***************************************
 *
 * QG group algebra, OK ring of ints, EK elliptic curve over K, h given by AutomorphismGroup,
 * l rational prime (we consider the l-part of ETNC), q rational prime (we consider reductions of EK
 * with respect to primes v over q).
 *
 * Computes eps( \oplus_{v|q} \bar E_{ns}(k_v)_l ) = \ind_D^G( eps( \bar E_{ns}(k_v) ) . 
 * If l divides #D and #E_{ns}(k_v)_l then cohomology is not perfect. Even if l does not
 * divide #D, we are not able to compute the Euler characteristic if the reduction is bad because
 * we do not yet consider the D-structure in this case. But this is certainly possible with some
 * effort.
 *
 *******************************************************************************************/

intrinsic Eps_Ens(QG :: AlgGrp, OK :: RngOrd, E :: CrvEll, h :: Map, l :: RngIntElt, q :: RngIntElt) -> BoolElt, List
{}
    local EK, is_regular,is_computable, P, D, Frob, I, QD, EP, Pts, j, lPts, A, nr, alpha, cnt;

    // print "In Eps_Ens with l = ", l, "  q = ", q;

    EK := E`EK;
    P := Factorization(q*OK)[1,1];
    D := DecompositionGroup(OK, h, q);
    Frob := Frobenius(OK, h, q);
    I := InertiaGroup(OK, h, q);
    QD := InitGroupAlgebra(D);

    if Valuation(Discriminant(EK), P) eq 0 then
        EP := Reduction(EK, P);
	Pts, j := TorsionSubgroup(EP);
	lPts := SylowSubgroup(Pts, l);

        if IsDivisibleBy(#D, l) and IsDivisibleBy(#Pts, l) then
	    print "E_ns is not regular for l = ", l, " and q = ", q, "    #E_ns = ", #Pts;
	    is_regular := false; 
            return is_regular, [* H[1]!1 : H in QG`H *];
        end if;

        is_regular := true; 
	A := EnsRep(lPts, j, q);
	
	A := [A^Log(g, Frob, I) : g in D];
        nr := ReducedNorm(lPts, A, QD, l);
        alpha := Induction(QG, QD, nr);
        return is_regular, alpha;
     end if;
     if Valuation(Discriminant(EK), P) gt 0 then
         cnt := NrOfNsPoints(E, OK, P);

         if IsDivisibleBy(#D, l) and IsDivisibleBy(cnt, l) then
	    is_regular := false; 
            print "Cannot compute Eps_Ens for q = ", q, " because l = ", l, " divides #D = ", #D, " and #E_ns = ", cnt;
	    return is_regular, [* H[1]!1 : H in QG`H *];
         end if;

         if not IsDivisibleBy(cnt, l) then
	     is_regular := true; 
	     return is_regular, [* H[1]!1 : H in QG`H *];
         /* In all the following cases l does not divide #D, hence is regular */
         elif ( Valuation(cnt, l) eq 1 and GCD(l-1, #D) eq 1 ) then
             /* D acts trivially on the l-part of \bar E_{ns}(k_v) which is cyclic of order l */
	     is_regular := true; 
             nr := [* H[1]!1 : H in QD`H *];
	     nr[1] := QD`H[1][1] ! l;
	     alpha := Induction(QG, QD, nr);
	     return is_regular, alpha;
	 elif #D eq 1 then
	     is_regular := true; is_computable := true;
	     nr := [* H[1] ! l^Valuation(cnt, l) : H in QD`H *];
             alpha := Induction(QG, QD, nr);
             return is_regular, alpha;
	 else
             assert ReductionType(E, q) eq "Split multiplicative";
             F := ResidueClassField(OK, P);
             is_regular := true; is_computable := true;
             A := Matrix([[q]]);
             U := SylowSubgroup( UnitGroup(F), l );
             A := [A^Log(g, Frob, I) : g in D];
             nr := ReducedNorm(U, A, QD, l);
             alpha := Induction(QG, QD, nr);
             return is_regular, alpha;
	 end if;
    end if;
end intrinsic;

/******************************** Eps_Tam ***************************************
 *
 * QG group algebra, OK ring of ints, E elliptic curve over Q, h given by AutomorphismGroup,
 * l rational prime (we consider the l-part of ETNC), q rational prime (we consider reductions of EK
 * with respect to primes v over q).
 *
 * Computes eps( \oplus_{v|q} E(K_v)/E_0(K_v) ) = ind_D^G( eps(E(K_v)/E_0(K_v)) ). 
 * If l divides #D and #E(K_v)/E_0(K_v) then the cohomology is not perfect. 
 *
 * In the perfect case we are only able to deal with the split multiplicative case.
 * This uses Prop. 5.4.
 *
 *******************************************************************************************/

intrinsic Eps_Tam(QG :: AlgGrp, OK :: RngOrd, E :: CrvEll, h :: Map, l :: RngIntElt, q :: RngIntElt) -> BoolElt, List
{}
    local is_regular, P, D, QD, nr, alpha, tam, type;

    P := Factorization(q*OK)[1,1];
    D := DecompositionGroup(OK, h, q);
    QD := InitGroupAlgebra(D);

    tam := NewTamagawaNumber(E, P);
    if IsDivisibleBy(#D, l) and IsDivisibleBy(tam, l) then
       print "E(K_v) / E_ns(K_v) is not regular for l = ", l, " and  q = ", q, "    # = ", tam;
       is_regular := false; 
       return is_regular, [* H[1]!1 : H in QG`H *];
    end if;

    type := ReductionType(E, q);
    assert type eq "Split multiplicative" or type eq "Good";
    is_regular := true; 
    nr := [* H[1]!1 : H in QD`H *];
    nr[1] := QD`H[1][1] ! l^Valuation(tam, l);
    alpha := Induction(QG, QD, nr);
    return is_regular, alpha;
end intrinsic;


/******************************** OldNrOfNsPoints ***************************************
 *
 * EK elliptic curve over K, OK ring of ints, P prime ideal of OK
 *
 * Computes the number of non-singular points in the reduction of EK mod P. One should have a 
 * function which computes the non-singular points of EK mod P as an abstract abelian group as
 * the magma function TorsionSubgroup. Then one could also consider the D-structure of this group
 * as we do in the case of good reduction.
 *
 * No longer used.
 *
 *******************************************************************************************/

intrinsic OldNrOfNsPoints(EK :: CrvEll, OK :: RngOrd, P :: RngOrdIdl) -> RngIntElt
{}
    local F, t, A, x, y, c, f, C;

    F, t := ResidueClassField(OK, P);
    if #F gt 10^7 then
        print "Cannot compute number of non singular points because #F = ", #F;
        return -1;
    end if;
    A<x,y> := AffineSpace(F, 2);
    c := [t(x) : x in Coefficients(EK)];
    f := y^2 + c[1]*x*y + c[3]*y - x^3 - c[2]*x^2 - c[4]*x - c[5];
    // C := Curve(A,y^2+y-x^3+x^2+10*x+20);
    C := Curve(A,f);
    return #RationalPoints(C, F) + 1 - #SingularPoints(C);
end intrinsic;

/******************************** NrOfNsPoints ***************************************
 *
 * E elliptic curve over Q, OK ring of ints, P prime ideal of OK
 *
 * Computes the number of non-singular points in the reduction of EK mod P. One should have a
 * function which computes the non-singular points of EK mod P as an abstract abelian group as
 * the magma function TorsionSubgroup. Then one could also consider the D-structure of this group
 * as we do in the case of good reduction.
 *
 * One should also cover the case of additive reduction !!!
 *
 *******************************************************************************************/

intrinsic NrOfNsPoints(E :: CrvEll, OK :: RngOrd, P :: RngOrdIdl) -> RngIntElt
{}
    local F, t, p, EK;

    F, t := ResidueClassField(OK, P);
    p := Characteristic(F);
    t := ReductionType(E, p);
    assert t eq "Split multiplicative" or t eq "Nonsplit multiplicative" or t eq "Good" ;
    if t eq "Split multiplicative" then
        return #F - 1;
    elif t eq "Nonsplit multiplicative" then
        return #F + 1;
    elif t eq "Good" then
        EK := BaseChange(E, Order(P));
        return #Reduction(EK, P);
    end if;
end intrinsic;

/******************************** NrOfNsPoints ***************************************
 *
 * E elliptic curve over Q, p rational prime
 *
 * Computes the number of non-singular points in the reduction of E mod p. 
 *
 * One should also cover the case of additive reduction !!!!!
 *
 *******************************************************************************************/

intrinsic NrOfNsPoints(E :: CrvEll, p :: RngIntElt) -> RngIntElt
{}
    local t;

    t := ReductionType(E, p);
    if t eq "Good" then
       return #EllipticCurve([GF(p)!x : x in ElementToSequence(E)]);
    end if;
    assert t eq "Split multiplicative" or t eq "Nonsplit multiplicative";
    if t eq "Split multiplicative" then
        return p - 1;
    else
        return p + 1;
    end if;
end intrinsic;


/******************************** Induction ***************************************
 *
 * QG group algebra, QD group algebra, where D is a subgroup of G, nr an element of the center
 * of QD, i.e.a list of elements in the character fields of QD. 
 *
 * Computes the induction on the level of the centers.
 *
 *******************************************************************************************/

intrinsic Induction(QG :: AlgGrp, QD :: AlgGrp, nr :: List) -> List
{}
    local E, tau, beta, i, chi, v, j, alpha, iota;

    D := Group(QD);
        E := CommonSplittingField(QG, QD);

        Tau := EmbeddingsToE(QD, E);
        beta := [**];
        for i:=1 to #QG`H do
            chi := Restriction(QG`H[i][2], D);
            v := [];
            for j:=1 to #nr do
                v cat:= [tau(nr[j])^ScalarProduct(D, tau, QD`H[j][2], QD`H[j][3], chi) : tau in Tau[j]];
            end for;
            Append(~beta, &*v);
        end for;

        Tau := EmbeddingsToE(QG, E);
        alpha := [**];
        for i:=1 to #beta do
            iota := Tau[i][1]; OKi := QG`H[i][4];
            Append(~alpha, (iota^-1)(beta[i]) );
        end for;
    return alpha;
end intrinsic;


/******************************** ReducedNorm ***************************************
 *
 * C is a finite G-module, A[g] \in M_n(Z) gives the action of g in G on C with respect to
 * the generators of C, l is a rational prime (we always consider l-parts)
 *
 * Computes eps( C ) as an element of the centre of QG.
 *
 *******************************************************************************************/

intrinsic ReducedNorm(C :: GrpAb, A :: SeqEnum, QG :: AlgGrp, l :: RngIntElt) -> List 
{}
    local G, rank, B, i, kern, lambdas, s, lam, rho, R, F, W, Q, BasisF, BasisQ, theta, S, nr;

    G := Group(QG);
    gens := [C.i : i in [1..#Generators(C)]];
    if #gens eq 0 then
        gens := [ Id(C) ];
    end if;
    rank := #gens;
        B := A[1];
        for i:=2 to #A do
            B := HorizontalJoin(B, A[i]);
        end for;
        B := HorizontalJoin(B, DiagonalMatrix([Order(P) : P in gens]));
        kern := [ ElementToSequence(x) : x in Basis(Kernel( Transpose(B))) ];
        lambdas := [];
        for i:=1 to #kern do
            lambda := [];
            for s := 1 to rank do
                lam := [];
                for t:=1 to #G do
                    Append(~lam, kern[i][s+(t-1)*rank]);
                end for;
                Append(~lambda, lam);
            end for;
            Append(~lambdas, lambda);
        end for;
        rho := RegularRep(QG);
        R := RegularModule(rho);
        F := R;
        for i:=2 to rank do
            F := ZGModuleDirectSum(F, R);
        end for;
        W := [&cat lambda : lambda in lambdas];
        Q := ZGSubModule(F, W);

        BasisF := LocalBasis(F, l);
        BasisQ := LocalBasis(Q, l);
        theta := IdentityMatrix(Rationals(), #G*rank);
        S := QGMatrix(QG, theta, Q, BasisQ, F, BasisF);
        nr:= NewtonReducedNorm(S);
    return nr;
end intrinsic;



/******************************** Log ***************************************
 *
 * Return i, if g = Frob^i mod I.
 *
 *******************************************************************************************/

intrinsic Log(g :: GrpPermElt, Frob :: GrpPermElt, I :: GrpPerm) -> RngIntElt
{}
    local i;

    for i:=0 to Order(Frob) do
        if g^-1 * Frob^i in I then
	    return i;
	end if;
    end for;
end intrinsic;

/******************************** EnsRep ***************************************
 *
 * lPts are generators of the l-part of the non-singular points, j is a map which gives actual
 * points on the reduced elliptic curve, q is a rational prime.
 *
 * Computes the representation matrix of the Frobenius.
 *
 * Presently we are only able to use this function in the case of good reduction !!!
 * If we could compute the non-singular points also in the case of bad reduction, then we we can
 * use this function also in this case.
 *
 *******************************************************************************************/

intrinsic EnsRep(lPts :: GrpAb, j :: Map, q :: RngIntElt) -> Mtrx
{}
    local A, G, gens, r, g, M;

    
    // gens := Generators(lPts);
    gens := [lPts.i : i in [1..#Generators(lPts)]];
    if #gens eq 0 then
         return Matrix([[1]]);
    end if;
    M := [ EnsLog(FrobAction(j(P), q), lPts, j) : P in gens ];
    return Transpose( Matrix(M) );
end intrinsic;
						

/******************************** FrobAction ***************************************
 *
 * P is a point on E / F_q, q a rational prime. Returns P^Frob.
 *
 *******************************************************************************************/

intrinsic FrobAction(P :: PtEll, q :: RngIntElt) -> PtEll
{}
    return Parent(P) ! [x^q : x in ElementToSequence(P)];
end intrinsic;
 
/******************************** GAction ***************************************
 *
 * P is a point on E / K, g in Gal(K/Q), h the map from Automorphismgroup or ArtinRepresentations.
 * Returns P^g.
 *
 *******************************************************************************************/

intrinsic GAction(P :: PtEll, g :: GrpPermElt, h :: Map) -> PtEll
{}
    return Parent(P) ! [h(g)(x) : x in ElementToSequence(P)];
end intrinsic;


/******************************** EnsLog ***************************************
 *
 * P is a point in the lPart of the non-singular points, lPts the group of the non-singular
 * points (as abstract abelian group), 
 * j the map which gives the coordinates of points on the reduced curve.
 * Solves the discrete log problem (in the most naive way).
 *
 *******************************************************************************************/

intrinsic EnsLog(P :: PtEll, lPts :: GrpAb, j :: Map) -> SeqEnum
{}
    local q;

    for q in lPts do
        if j(q) eq P then
            return ElementToSequence(q);
        end if;
    end for;
end intrinsic;


/******************************** HeightPairing ***************************************
 *
 * P, Q are points in E(K), K/Q a number field. Computes the height pairing. 
 * Attention: Note the normalization !!!
 * This is not what we need in BSD for E/K.
 *
 *******************************************************************************************/

intrinsic HeightPairing(P :: PtEll, Q :: PtEll) -> FldReElt
{}
    return CanonicalHeight(P+Q) - CanonicalHeight(P) - CanonicalHeight(Q);
end intrinsic;

/******************************** HeightPairing ***************************************
 *
 * P, Q are points in E(K), K/Q a number field. Computes the height pairing.
 * Attention: Note the normalization !!!
 * That the form which is needed for BSD over number fields or the ETNC for E/K. 
 *
 *******************************************************************************************/

intrinsic HeightPairing(P :: PtEll, Q :: PtEll, K :: FldNum) -> FldReElt
{}
    local prec;

    prec := Precision(RealField());
    return ( CanonicalHeight(P+Q:Precision:=prec) - 
             CanonicalHeight(P:Precision:=prec) - 
             CanonicalHeight(Q:Precision:=prec) ) * Degree(K) / 2;
end intrinsic;

/******************************** HeightPairing1 ***************************************
 *
 * P, Q are points in E(K), K/Q a number field. Computes the height pairing.
 * Attention: Note the normalization !!!
 *
 * This is NOT the normalization which is needed for ETNC.
 *
 *******************************************************************************************/

intrinsic HeightPairing1(P :: PtEll, Q :: PtEll, K :: FldNum) -> FldReElt
{}
    local prec;

    prec := Precision(RealField());
    return ( CanonicalHeight(P+Q:Precision:=prec) -
             CanonicalHeight(P:Precision:=prec) -
             CanonicalHeight(Q:Precision:=prec) ) / 2;
end intrinsic;



/******************************** MyRegulator ***************************************
 *
 * Computes the regulator of a sequence of points of an elliptic curve defined over
 * a number field. Uses the height pairing with respect to the number field K, when the points are
 * given as points in E(K).
 *
 *******************************************************************************************/

intrinsic MyRegulator(pts :: SeqEnum) -> FldReElt
{}

    if #pts eq 0 then return 1; end if;
    K := Parent( Coordinate(pts[1], 1) );
    R := [ [HeightPairing(P,Q,K) : P in pts] :  Q in pts];
    R := Matrix(R);
    return( Determinant(R) );

end intrinsic;


/******************************** IdealPartIsTrivial ***************************************
 *
 * log as returned by K0RelLog. Checks whether log is trivial modulo torsion.
 *
 *******************************************************************************************/

intrinsic IdealPartIsTrivial(log :: List) -> BoolElt
{}

    return &and[ P eq 1*Order(P) : P in log[1] ];
end intrinsic;



/********************************  ReadReal ***************************************
 *
 * Very basic input routine for reading reals.
 *
 *******************************************************************************************/

intrinsic ReadReal(str :: MonStgElt) -> FldReElt
{}
     local s, a, b;

     read s, str;
     a:=Substring(s, 1,Index(s, ".")-1);
     b:=Substring(s, Index(s, ".")+1, #s);
     return( RealField() ! ( StringToInteger(a) + StringToInteger(b) / 10^#b ) );
end intrinsic;


/********************************  PrintList ***************************************
 *
 * Very basic output routine for lists.
 *
 *******************************************************************************************/

intrinsic PrintList(str :: MonStgElt, L :: Any) 
{}
    local i, len;

    printf(str);
    len := #L;
    printf "(";
    for i:=1 to len-1 do
        printf " %o,", L[i];
    end for;
    if len gt 0 then
        printf " %o )\n", L[len];
    else
        printf " )\n";
    end if;
end intrinsic;




/********************************  PrintList ***************************************
 *
 * pols is a list of polynomials which define tame Dp-extension K/Q, Es is a list of
 * elliptic curves with split multiplicative reduction. Looks for examples where we can
 * perform the numerical verification of ETNC at l = p. In particular, we need the analytic
 * rank of E/K to be trivial and l nmid I_2.
 *
 *******************************************************************************************/

intrinsic FindGoodDpExample(pols :: SeqEnum, Es :: SeqEnum) -> SeqEnum
{}
    local GoodExamples, pol, L, A, K, G, QG, act, h, E, Lvec;

    GoodExamples := [];

    for pol in pols do
        L := NumberField(pol);
        A := ArtinRepresentations(L);
        K :=L`artinrepdata`K;
	OK := MaximalOrder(K);
        G := L`artinrepdata`G;
        QG := InitGroupAlgebra(G);
        act := L`artinrepdata`act;
        h := map<G->Parent(act(G.1)) | g:->act(g^-1)>;
        p := Integers()!(#G / 2);
        prec := 5;
        for E in Es do
            S := ComputeS(E, OK);
            good := true;
            for q in S do
                Iq := InertiaGroup(OK, h, q);
                good := good and not IsDivisibleBy(#Iq, p);
            end for;
            c := &*TamagawaNumbers(E);
	    print CremonaReference(E);
            good := good and (GCD(Conductor(E), Discriminant(OK)) eq 1 and
                              not IsDivisibleBy(c, p) and 
                              not p in ComputeS(E, OK));
	    if good then
                   Lvec := InitLSeries(QG, E, K, prec);
                   if AnalyticRankIsZero(Lvec) then
	               print "Good example found";
	               Append(~GoodExamples, [* pol, E *]);
	           end if;
	    end if;
	end for;
    end for;
    return GoodExamples;
end intrinsic;


	    


