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

ZX<x> := PolynomialRing( Integers() );
load "D5Mipos.m";

DB := CremonaDatabase();
prec := 20;
SetDefaultRealFieldPrecision(prec);

/************************************************************************************
*                      D_5-Beispiel                                            *
************************************************************************************/

E := EllipticCurves(DB, "11A")[1];
pol := D5_pols[1,1];
K := SplittingField(NumberField(pol));
/* Compute an integral normal basis element. */
X := 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 outside 2 */

/* Compute general information */
CharacterTable(G);
EK := BaseChange(E, K);

Factorization(401*OK);

P := Factorization(2*OK)[1,1];
#Reduction(EK, P);
P := Factorization(5*OK)[1,1];
#Reduction(EK, P);
P := Factorization(401*OK)[1,1];
#Reduction(EK, P);
P := Factorization(11*OK)[1,1];
NrOfNsPoints(E, OK, 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;

/* Compute the equivariant regulator using pre-computed points */
F, iF := FixedField(K, [h(g) : g in Kernel(QG`X[2,1])]);
d := Discriminant( MaximalOrder( F ) );
d;
Et := QuadraticTwist(E, d);
EF := BaseChange(E, F);
EK := BaseChange(E, K);

FX := PolynomialRing(F);
roots := Roots(FX ! x^2 - 401);
a := roots[1,1];

P1 := EK ! [iF(74/9), iF(53/54*a - 1/2)];
P2 := EK ! [iF(6), iF(1/2*a - 1/2)];

-P1 eq EK ! [iF(74/9), iF(53/54*(-a) - 1/2)];
-P2 eq EK ! [iF(6), iF(1/2*(-a) - 1/2)];

RegMat := Matrix(RealField(), 2, 2,
                 [HeightPairing(P1,P1,K), HeightPairing(P1,P2,K),
                  HeightPairing(P2,P1,K), HeightPairing(P2,P2,K)]);
Rvals := [ [1], [Determinant(RegMat)], [1, 1]];


/* 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";
5
RatZ := MakeRational(Z, QG, d, Round(prec/2));

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 */
TateShafarevicGroup(E, QG, Lvals, Rvals, OK);

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;

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

for l in HP do
    CheckHardPrime(QG, E, S, l);
    print "*********************************************************************";
end for;