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

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

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

E := EllipticCurves(DB, "73A")[1];
pol := x^10 - 2*x^9 - 20*x^8 + 2*x^7 + 69*x^6 - x^5 - 69*x^4 + 2*x^3 + 20*x^2 - 2*x - 1;
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);

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 */
Z := ComputeZ(QG, E, Lvals, R, Ovals, Rvals);
print "Z = ", Z;
readi d, "Input a denominator";
2
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);
2304
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(3*OK)[1,1];
#Reduction(EK, P);
P := Factorization(5*OK)[1,1];
#Reduction(EK, P);
P := Factorization(73*OK)[1,1];
NrOfNsPoints(E, OK, P);
P := Factorization(401*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 = 5.
   From the above computation we have 
           eta := u * xi_5^-1 = ( 1/2, 18, -12*$.1 + 8 ) = RatZ.
*/
OF := MaximalOrder( Parent(RatZ[3,1]) );
P := Factorization(5*OF)[1,1];
eta := RatZ;
Valuation(OF!(eta[1,1]*eta[2,1]) - eta[3,1], P);


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