Attach("RelAlgKTheory.m");
Attach("INB.m");
Attach("etnc.m");

ZX<x> := PolynomialRing( Integers() );
DB := CremonaDatabase();
prec := 30;
SetDefaultRealFieldPrecision(prec);

/************************************************************************************
*                      The D_7 Example of the paper EllETNC I
************************************************************************************/

E := EllipticCurves(DB, "11A")[1];
pol := x^14 - 2*x^13 - 25*x^12 + 69*x^11 + 161*x^10 - 632*x^9 - 147*x^8 + 2146*x^7 - 1171*x^6 -
    2669*x^5 + 2682*x^4 + 667*x^3 - 1466*x^2 + 336*x + 49;
L := NumberField(pol);
A := ArtinRepresentations(L);
K :=L`artinrepdata`K;

/* Compute an integral normal basis element. */
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);   /* delta generates an integral normal basis */

/* Compute general information */
EK := BaseChange(E, K);
D := Discriminant(OK);
N := Conductor(E);
Factorization( D );
Factorization( N );
Factorization(N*OK);
Factorization(577*OK);

P := Factorization(N*OK)[1,1];
NrOfNsPoints(E, OK, P);

P := Factorization(D*OK)[1,1];
#Reduction(EK, P);

TorsionSubgroup(E);
TorsionSubgroup(EK);

TamagawaNumbers(E, OK);


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 conjecturally is ", orders;
print "Lvals = ", Lvals;

/* The Mordell-Weil group is torsion, so the regulator is trivial */
Rvals := [ [ComplexField()!1 : i in [1..#z]] : z in Lvals];

/* Compute the ratios of L-values, periods and resolvents and check rationality */
/* Note that the values for Z depend on the choice of the integral normal basis element delta.
   More explicitly, Z may change by the reduced norm of an element lambda in ZG^\times.
   Note that such reduced norms always satisfy the explicit congruences !

   Depending on the choice of delta there may be problems with rounding.
*/
Z := ComputeZ(QG, E, Lvals, R, Ovals, Rvals);
print "Z = ", Z;
readi d, "Input a denominator";
5
RatZ, err := MakeRational(Z, QG, d, Round(4));
print "u = ", RatZ;
print "Maximal rounding error: ", err; 
g := MinimalPolynomial(RatZ[3,1]);
g;
Roots( PolynomialRing(Parent(RatZ[3,1])) ! g);


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

/* Compute a conjectural order for the Tate-Shafarevic group */
E`ts := TateShafarevicGroup(E, QG, Lvals, Rvals, OK);
244140625
print "Conjectural order of Sha(E/K): ", E`ts;

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 "S = ", S;
print "HP = ", HP;

P := Factorization(2*OK)[1,1];
#Reduction(EK, P);
P := Factorization(5*OK)[1,1];
#Reduction(EK, P);
P := Factorization(7*OK)[1,1];
#Reduction(EK, P);
P := Factorization(11*OK)[1,1];
NrOfNsPoints(E, OK, P);
P := Factorization(577*OK)[1,1];
#Reduction(EK, P);



/* Check primes which are not in HP */
CheckEasyPrimes(RatZ, HP);

/* Check the difficult primes */
for l in HP do
    CheckHardPrime(QG, E, S, l);
    print "*********************************************************************";
end for;

/* Finally check the explicit congruence in the case l = 7.
   From the above computation we have 
           eta := u * xi_7^-1 = ( -1/5, -5, 25*$.1 ) = RatZ.
*/
l := 7;
OF := MaximalOrder( Parent(RatZ[3,1]) );
P := Factorization(l*OF)[1,1];
eta := RatZ;
Valuation(OF!(eta[1,1]*eta[2,1]) - eta[3,1], P);


/* Or very explicit */
C<z> := CyclotomicField(l);
K3 := sub<C | z+z^-1>;
K3<beta> := sub<C | z+z^-1>;
OK3 := MaximalOrder(K3);
P := Factorization(l*OK3)[1,1];
Valuation(1 - 25*beta, P);