r"""

This file contains various stpes to verify numerically the eTNC_p conjecture
for an elliptic curve defined over Q and a dihedral extension F/Q of degree
2p with p an odd prime. The curve must have rank 1 over the quadratic
extension K/Q within F.

This should work for real and imaginary case and also for the curve having
rank 0 or 1 over Q, although the cases differ quite a bit in how they are
computed.

The script often calls magma to compute twisted L-series, where the very
good implementation by Tim Dokchitser is used. The rest is in sage and,
underlying, pari-gp.

Reference:
David Burns, Daniel Macias Castillo, Christian Wuthrich, 
'On Mordell-Weil groups and congruences between derivatives of 
Hasse-Weil L-functions', preprint 2012.

AUTHOR:

- Chris Wuthrich 05/12

How to use it:

This file contains one class called ``etnc_example`` and it can be initiated
with a elliptic curve (or just its Cremona label) and a polynomial f. The
Galois splitting field of the degree p polynomial f should be a dihedral
extension. The file also contains at the end a function dihedral_fields(p)
which returns a list of a few such polynomials of small discriminant found
in the tables by Klueners and Malle:
http://www.math.uni-duesseldorf.de/~klueners/minimum/

Once the class is initiated, step by step, one can try to compute information
about this example. Or one can use the method do_all() to attempt to evaluate
all of them in one go.

The intermediate steps (in this order) are

* check_conditions()
* over_Q()
* over_K()
* over_L()
* find_Q()
* over_F()
* eq_regs()
* eq_L_values()
* truncated_L_values()
* check_congruence()

Any of these will print some information and set a few attributes to the class.

The procedures over_? just determine the Mordell-Weil group over this field and
verify numerically the BSD conjecture, i.e. it computes the analytic order of
the Tate Shafarevich group.

In some cases the function over_L() will fail as it needs to find a point in E(L)
which is not defined over Q. If such a point is known from another source, then
one can feed this to the algorithm (see second example below)

WARNING : These are no algorithms that work in all cases. Sometimes a precision 
loss will cause a problem, sometime there could be an implicit assumption that
is violated etc. When the output is that the congruence does not hold, it is 
usually because of an intermediate error that can be corrected.

FIRST EXAMPLE ::

    sage: R.<X> = QQ[]
    sage: d = etnc_example("43a1",X^3-4*X+2,prec=16)
    sage: d.do_all()
    CONDITIONS HOLD
    -- start over_Q 
    bsd over Q asserts 1.00000000000000 = 1 * Sha 
    -- start over_K 
    saturate at 2
    bsd over K=Q(sqrt(37)) asserts 3.99999999999999 = 1 * Sha 
    -- start find_points_over_L 
    magma finds the point  (-2*t^2 - t + 7 : 5*t^2 + 3*t - 19 : 1)
    saturate at 2
    (2 : -4 : 1) is divisible by 2
    saturate at 2
    (-1 : 0 : 1) is divisible by 2
    saturate at 2
    saturate at 3
    saturate at 5
    saturate at 7
    this may not yet be a basis of E(L). We assume it however for the rest !
    -- start eq_L_values 
    the L-values twisted by 2-dimensional representations are  [13.0341187850503]
    -- start over_L 
    bsd over L of discriminant 2^2 * 37 asserts 0.999999999999999 = 1 * Sha 
    -- start find_Q 
    norm of P is 4 * R
    the point Q is 1 * P + -1 * R = (-30/2861*s^5 + 441/11444*s^4 + 1400/2861*s^3 - 3945/2861*s^2 - 28773/5722*s + 18877/2861 : 183/11444*s^5 - 601/11444*s^4 - 2135/2861*s^3 + 11317/5722*s^2 + 42949/5722*s - 34724/2861 : 1)
    -- start over_F 
    saturate at 2
    saturate at 3
    bsd over F of discriminant 2^4 * 37^3 asserts 3.99999999999998 = 1 * Sha 
    -- start eq_regs 
    the equivariant regulators are H_eins = 0.12563301417497529853141758293393737263798407454502063305996, H_eps = 1 and H_rho = [5.3020765772838866517814520217761932242668258787941996311024]
    -- start alg_L_values 
    permutation found. Setting the trace to 1 with error 1.44328993201270e-15
    the equivariant algebraic L-values are Qtilde_one = 1/2, Qtilde_eps = 8, Qtilde_rho has minimal polynomial  x - 1
    -- start truncated_L_values 
    the truncated equivariant L-values are Q_one = 95/74, Q_eps = 4, sum Q_rho = -38/37
    etnc_3 holds numerically for this example as 190/37 == 76/37 (mod 3)

 
SECOND EXAMPLE ::

    sage: R.<X> = QQ[]
    sage: f = X**5-X**4-5*X**3+4*X**2+3*X-1 
    sage: d = etnc_example("37a1",f, prec=16)
    
First we check if all the conditions in the paper hold::

    sage: d.check_conditions(verbose=True)
    conditions for section 1 hold
    all conditions hold, except maybe (h) and (i)
    True
  
Then we check BSD over Q::
    
    sage: d.over_Q()
    -- start over_Q 
    bsd over Q asserts 2.00000000000000 = 2 * Sha 

Then over K, the quadratic field in F/Q::

    sage: d.over_K()
     -- start over_K 
    saturate at 2
    bsd over K=Q(sqrt(401)) asserts 4.00000000000000 = 4 * Sha 

Now, the news points are quite difficult to find. However
the Stark-Heegner points for this examples were computed by
Darmon and his collaborators. So we feed one to the algorithm.
L here is a degree p=5 extension within F ::

    sage: EL = d.EL
    sage: RL.<x> = d.L[]
    sage: f_x = x^5 - 73/9*x^4 + 1195/81*x^3 - 173/81*x^2 - 976/81*x + 527/81
    sage: x_P = f_x.roots()[0][0]
    sage: P = EL.lift_x(x_P)
    sage: P
    (-7/9*t^4 + t^3 + 3*t^2 - 13/3*t + 16/9 : 16/9*t^4 - 82/27*t^3 - 187/27*t^2 + 338/27*t - 88/27 : 1)
    sage: d.find_points_over_L(known_point = P)
    -- start find_points_over_L 
    saturate at 2
    saturate at 3
    saturate at 5
    saturate at 7
    this may not yet be a basis of E(L). We assume it however for the rest !
    [(-7/9*t^4 + t^3 + 3*t^2 - 13/3*t + 16/9 : 16/9*t^4 - 82/27*t^3 - 187/27*t^2 + 338/27*t - 88/27 : 1), (1203457/3207681*t^4 - 1050725/3207681*t^3 - 1104262/1069227*t^2 + 1745119/3207681*t + 2637008/3207681 : -814604161/5744956671*t^4 + 4605928994/5744956671*t^3 - 390085805/1914985557*t^2 - 1938248974/5744956671*t - 3802270868/5744956671 : 1), (-3*t^4 + 5*t^3 + 12*t^2 - 19*t + 4 : -13*t^4 + 22*t^3 + 51*t^2 - 86*t + 18 : 1)]

The last warning just says that we have not used the Siksek bound
to verify that we have saturated the Mordell-Weil group enough to
be certain the have the full Mordell Weil group. For the entc it
only matters that we have p-saturated.

Now we can check BSD over L::

    sage: d.over_L()   # long time
    -- start eq_L_values 
    the L-values twisted by 2-dimensional representations are  [58.5757519880864, 7.99398779371418]
    -- start over_L 
    bsd over L of discriminant 401^2 asserts 31.9999999999997 = 32 * Sha 

Now, the first part of our paper says that there exists a point Q in E(F)
such that E(F) is Z[G]Q, up to torsion and maybe some index 2::

    sage: d.find_Q()
    -- start find_Q 
    norm of P is 1 * R
    the point Q is 1 * P + 0 * R = (99196036141/331373902316634*s^9 + 285048480883/331373902316634*s^8 - 19179474261305/994121706949902*s^7 - 16304483654039/331373902316634*s^6 + 12967315524831/36819322479626*s^5 + 705748434839215/994121706949902*s^4 - 2450104724313527/994121706949902*s^3 - 1065793132491143/331373902316634*s^2 + 324947801308067/55228983719439*s + 3656691390456553/994121706949902 : -2395659744361/2982365120849706*s^9 - 3308484787691/1491182560424853*s^8 + 76427108355079/1491182560424853*s^7 + 187153901511740/1491182560424853*s^6 - 2718347796562069/2982365120849706*s^5 - 294397070273263/165686951158317*s^4 + 9001686574880344/1491182560424853*s^3 + 24348622303169119/2982365120849706*s^2 - 2124663860560766/165686951158317*s - 29998483307659769/2982365120849706 : 1)

Next we check BSD over F::

    sage: d.over_F()  # long time
    -- start over_F 
    saturate at 2
    saturate at 3
    saturate at 5
    bsd over F of discriminant 401^5 asserts 1023.99999999998 = 1024 * Sha 

The equivariant regulators H_rho::

    sage: d.eq_regs()
    -- start eq_regs 
    the equivariant regulators are H_eins = 0.10222281647993768047177219951388404321907640456170592849849, H_eps = 1 and H_rho = [4.4661072698670985345601764762836972799865563344802408782635, 32.725292875432142475197512055286549091997859324196023957098]

    sage: prod(d.h_rhos)*d.h_eins/d.reg_L
    1.9999999999999999999999999999999999999999999999999999999935

The last line check that the 4*regulator of E(L) is H_eins * prod(H_rho for dim(rho)=2 ).
The twisted L-values will be computed. Note that we set the precision to 16 in the beginning::

    sage: d.alg_L_values()
    permutation found. Setting the trace to 8 with error 3.55271367880050e-14
    the equivariant algebraic L-values are Qtilde_one = 1, Qtilde_eps = 4, Qtilde_rho has minimal polynomial  x - 4

The lat information are the algebraic values \tilde{Q}. We print a likely
minimal polynomial for the ones of dimension 2.
Note that at this step there is some guessing involved. Magma gives us the
twisted L-values and our previous computation the regulators. These two
list however need not be ordered in the same way. So check all combinations
until we found one with trace very close to a rational number of small denominator.
The second to last line tell us about the error.
In this particular example Q_tildes were already rational numbers.

Finally, we need to truncate the L-series and Gauss sums at ramified places::

    sage: d.truncated_L_values()
    -- start truncated_L_values 
    the truncated equivariant L-values are Q_one = -384/401, Q_eps = 4, sum Q_rho = -3072/401

and then check if the predicted congruences hold ::

    sage: d.check_congruence()
    etnc_5 holds numerically for this example as -1536/401 == 6144/401 (mod 5)

THIRD EXAMPLE::

    sage: R.<X> = QQ[]
    sage: d= etnc_example("50b1",X^3 - X^2 - 5*X + 4,prec=16)
    sage: d.do_all() # long time
    CONDITIONS HOLD
    -- start over_Q 
    bsd over Q asserts 0.200000000000000 = 1/5 * Sha 
    -- start over_K 
    found point on twist (-447371/1681 : 1403243866/68921 : 1)
    gives point (-110738/112627 : 56900929/5568954642*s0 + 175802378/2784477321 : 1) over K=QQ(sqrt(469))
    saturate at 2
    bsd over K=Q(sqrt(469)) asserts 0.199999999999999 = 1/5 * Sha 
    -- start find_points_over_L 
    magma finds the point  (2*t^2 - 6*t + 1 : 8*t^2 - 24*t + 14 : 1)
    saturate at 2
    saturate at 3
    saturate at 5
    saturate at 7
    this may not yet be a basis of E(L). We assume it however for the rest !
    -- start eq_L_values 
    the L-values twisted by 2-dimensional representations are  [16.5156716492158]
    -- start over_L 
    bsd over L of discriminant 7 * 67 asserts 1.00000000000000 = 1 * Sha 
    -- start find_Q 
    the point Q was obtained from S = (2) and is equal to (15491/898537*s^5 - 159559/3594148*s^4 - 1034283/898537*s^3 + 8924971/3594148*s^2 + 28301039/1797074*s - 45735689/3594148 : 3542155/60201979*s^5 - 37240335/481615832*s^4 - 662158945/180605937*s^3 + 7056897217/1444847496*s^2 + 30645994063/722423748*s - 58833455135/1444847496 : 1) 
    -- start over_F 
    saturate at 2
    (29379072145/20966043639758*s^5 - 33066111595/10483021819879*s^4 - 974425433385/10483021819879*s^3 + 586889666111/2995149091394*s^2 + 1899468804082/1497574545697*s - 11798535328161/10483021819879 : -23168047852649705/12694499136957034082*s^5 + 41294512992631735/12694499136957034082*s^4 + 5055905094999026677/38083497410871102246*s^3 - 3475236082820115220/19041748705435551123*s^2 - 39757812344226769172/19041748705435551123*s + 586255880394452510/2720249815062221589 : 1) is divisible by 2
    saturate at 2
    (8143002327/20966043639758*s^5 - 1854038100/1497574545697*s^4 - 275166095955/10483021819879*s^3 + 1057590977769/20966043639758*s^2 + 543147921012/1497574545697*s + 722717031283/1497574545697 : 540014836544103/6347249568478517041*s^5 + 1994392699541259/12694499136957034082*s^4 + 5950405794621667/777214232874920454*s^3 + 557805912249165001/38083497410871102246*s^2 - 1029524953836511598/2720249815062221589*s - 3694246096679909206/2720249815062221589 : 1) is divisible by 2
    saturate at 2
    saturate at 3
    bsd over F of discriminant 7^3 * 67^3 asserts 5.00000000000002 = 5 * Sha 
    -- start eq_regs 
    the equivariant regulators are H_eins = 1, H_eps = 24.918817743633364873940975569685984917863261491634259562394 and H_rho = [12.502017857243808799928029364651848928000017580974916618997]
    -- start alg_L_values 
    permutation found. Setting the trace to 5/4 with error 4.44089209850063e-15
    the equivariant algebraic L-values are Qtilde_one = 1/5, Qtilde_eps = 1, Qtilde_rho has minimal polynomial  4*x - 5
    -- start truncated_L_values 
    the truncated equivariant L-values are Q_one = 110/469, Q_eps = 1, sum Q_rho = 1375/938
    etnc_3 holds numerically for this example as 110/469 == -1375/469 (mod 3)

FOURTH EXAMPLE (ramified extension of degree 5 over Q(i) for a curve of rank 0 over Q) ::

    sage: R.<x> = QQ[]
    sage: f = x^5 - 2*x^4 - 6*x^3 + 10*x^2 + 17*x - 12
    sage: d = etnc_example("21a1",f,prec = 32)
    sage: d.do_all()  #very long time (several hours)
    CONDITIONS HOLD
    -- start over_Q 
    bsd over Q asserts 0.2500000000000000000000000000000 = 1/4 * Sha 
    -- start over_K 
    found point on twist (-6 : 14 : 1)
    gives point (3/2 : 7/1368*s0 - 116/171 : 1) over K=QQ(sqrt(-4))
    saturate at 2
    bsd over K=Q(sqrt(-4)) asserts 0.2500000000000000000000000000000 = 1/4 *
    Sha 
    -- start find_points_over_L 
    magma finds the point  (-1/8*t^4 - 1/8*t^2 - 3/2*t - 1 : 39/32*t^4 -
    7/32*t^3 - 213/32*t^2 - 127/32*t + 49/8 : 1)
    saturate at 2
    saturate at 3
    saturate at 5
    saturate at 7
    this may not yet be a basis of E(L). We assume it however for the rest !
    -- start eq_L_values 
    the L-values twisted by 2-dimensional representations are 
    [8.608827844782438006881130813540, 7.532156342610583068208326962226]
    -- start over_L 
    bsd over L of discriminant 2^4 * 19^4 asserts
    64.00000000000000000000000000027 = 16 * Sha 
    -- start find_Q 
    the point Q was obtained from S = (4, 3) and is equal to
    (-800383197/277206133477504*s^9 - 860396345/39600876211072*s^8 +
    9153474137/39600876211072*s^7 + 421083260413/277206133477504*s^6 -
    219481985381/25200557588864*s^5 - 11335743186093/277206133477504*s^4 +
    40612740828241/277206133477504*s^3 + 139478602959811/277206133477504*s^2
    - 65226789478787/69301533369376*s - 340192091570241/69301533369376 :
    -1167290653/1663236800865024*s^9 + 1644685747/79201752422144*s^8 +
    13975484095/79201752422144*s^7 - 730407950067/554412266955008*s^6 -
    1774054650205/151203345533184*s^5 + 43834828000181/1663236800865024*s^4
    + 303780729547741/1663236800865024*s^3 +
    34764157743505/1663236800865024*s^2 - 570321570098807/415809200216256*s
    - 739362983558827/415809200216256 : 1) 
    -- start over_F 
    saturate at 2
    saturate at 3
    saturate at 5
    bsd over F of discriminant -1 * 2^10 * 19^8 asserts
    16384.00000000000000000000000014 = 512 * Sha 
    -- start eq_regs 
    the equivariant regulators are H_eins = 1, H_eps =
    1.4925984173890274178188752072631063537841265398716775762473 and H_rho =
    [15.523044529000670898117290401029177191873730624085000217421,
    1.9815345320563939359206623685543995364232610232645726932082]
    -- start alg_L_values 
    permutation found. Setting the trace to 48 with error
    5.048709793414475554635062817810e-29
    permutation found. With larger error
    2.382576922837997948024977792504e-6. We ignore it, its trace would be
    10485/197
    Warning: found more than one value. The chosen one might be wrong!!!
    the equivariant algebraic L-values are Qtilde_one = 1/4, Qtilde_eps = 1,
    Qtilde_rho has minimal polynomial  549255893*x^2 - 25968106226*x -
    15712796421
    -- start truncated_L_values 
    the truncated equivariant L-values are Q_one = 8/19, Q_eps = 24/19, sum
    Q_rho = -96
    etnc_5 holds numerically for this example as 192/361 == 192 (mod 5)

FIFTH EXAMPLE (Stark-Heegner point over a degree seven extension of a real quadratic field for the curve 37a1 of rank1 over Q) ::

    fi = dihedral_fields(7, "real");
    f = fi[0]; f
    d = etnc_example("37a1",f, prec= 16)
    d.over_Q()
    -- start over_Q 
    bsd over Q asserts 2.00000000000000 = 2 * Sha 

    d.over_K()
    -- start over_K 
    saturate at 2
    bsd over K=Q(sqrt(577)) asserts 4.00000000000001 = 4 * Sha 
    
    R.<x> = d.L[]
    fx = x^7-29*x^6+245*x^5-633*x^4+515*x^3-15*x^2-18*x-1
    xP = fx.roots()[0][0]
    P = d.EL.lift_x(xP)
    d.find_points_over_L(known_point=P)
    -- start find_points_over_L 
    saturate at 2
    saturate at 3
    saturate at 5
    saturate at 7
    this may not yet be a basis of E(L). We assume it however for the rest !
    [(2*t^6 - 4*t^5 - 14*t^4 + 17*t^3 + 30*t^2 - 7*t - 3 : 2*t^6 - 19*t^4 - 2*t^3 + 32*t^2 + t - 5 : 1), (2*t^6 - 2*t^5 - 14*t^4 + 8*t^3 + 26*t^2 - 3*t - 3 : -3*t^6 - 7*t^5 + 8*t^4 + 29*t^3 + 15*t^2 - 5*t - 3 : 1), (-331438/380689*t^6 + 472039/380689*t^5 + 2470111/380689*t^4 - 1397324/380689*t^3 - 5088084/380689*t^2 - 1579987/380689*t + 1611853/380689 : 109098724/234885113*t^6 - 36139678/234885113*t^5 - 1510274811/234885113*t^4 + 1229703263/234885113*t^3 + 3616665689/234885113*t^2 - 3385061968/234885113*t + 585716643/234885113 : 1), (105*t^6 - 242*t^5 - 660*t^4 + 1249*t^3 + 975*t^2 - 1337*t + 320 : 2916*t^6 - 6767*t^5 - 18245*t^4 + 35015*t^3 + 26701*t^2 - 37744*t + 9148 : 1)]

    sage: d.over_F(dont_saturate=True)
    -- start over_F 
    bsd over F of discriminant 577^7 asserts 16384.0000000021 = 16384 * Sha 

    sage: d.do_all()
    Already done.
    Already done.
    Points already found.
    Already done.
    Already done.
    Already done.
    Already done.
    -- start eq_regs 
    the equivariant regulators are H_eins =
    0.10222281647993768047177219951388404321907640456170592849849, H_eps = 1
    and H_rho =
    [3.0836851237573553962975584583954456723973726035554223935599,
    33.520071682824013363299006243826299571948977069703228460044,
    9.3011344262851043506212274587876816351574576345806444172552]
    -- start alg_L_values 
    permutation found. Setting the trace to 12 with error
    2.50466314355435e-13
    the equivariant algebraic L-values are Qtilde_one = 1, Qtilde_eps = 4,
    Qtilde_rho has minimal polynomial  264*x^3 - 2808*x^2 - 14013*x - 5321
    -- start truncated_L_values 
    the truncated equivariant L-values are Q_one = -578/577, Q_eps = 4, sum
    Q_rho = -6936/577
    etnc_7 holds numerically for this example as -2312/577 == 13872/577 (mod
    7)
"""
#****************************************************
#     Copyright (C) 2012 Chris Wuthrich
#       <christian.wuthrich@gmail.com>
#       GNU General Public License (GPL)
#****************************************************

# First a few functions that will be useful.

def generalised_Nv(E,v):
    r"""
    Given an elliptic curve over a number field and v a place,
    this returns the number of non-singular points on the reduction. 
    
    INPUT:
    
        - ``E`` - an elliptic curve over a number field
        - ``v`` - a place in this field
        
    OUTPUT:
        an natural number
        
    EXAMPLES ::
    
        sage: E = EllipticCurve("37a1")
        sage: K.<t> = QuadraticField(101)
        sage: EK = E.base_extend(K)
        sage: v = K.primes_above(7)[0]
        sage: generalised_Nv(EK,v)
        63

        sage: K.<t> = QuadraticField(103)
        sage: v = K.primes_above(37)[0]
        sage: EK = E.base_extend(K)
        sage: generalised_Nv(EK,v)
        1368

    """
    dav = E.local_data(v)
    qv = v.smallest_integer() ** v.residue_class_degree()
    if dav.has_good_reduction():
        Nv = E.reduction(v).cardinality()
    elif dav.has_split_multiplicative_reduction():
        Nv = qv - 1
    elif dav.has_nonsplit_multiplicative_reduction():
        Nv = qv + 1
    else: #additive
        Nv = qv
    return Nv


def nt_height(P):
    r"""
    Returns the canonical height of the point P on an elliptic curve over a number field.
    
    Rather than sage's usual normalisation, we make it dependent on the field here, so it
    is simply the degree times the height computed by sage
    
    EXAMPLES::
    
        sage: E = EllipticCurve("37a1")
        sage: K.<t> = QuadraticField(103)
        sage: EK = E.base_extend(K)
        sage: P = EK([0,0])
        sage: nt_height(P)
        0.10222281647993768047177219951388404321907640456170592849849                                                                                         
        sage: nt_height(E([0,0]))                                                                                                                             
        0.051111408239968840235886099756942021609538202280852964249243 
        
    """
    K = P.curve().base_field()
    return K.degree() * P.height(200)

def reg(pts):
    r"""
    Returns the regulator of a list of points on an elliptic curve over a number field
    
    Rather than sage's usual normalisation, we make it dependent on the field here, so it
    is simply the height is the degree times the height computed by sage.
    
    EXAMPLES::
    
        sage: E = EllipticCurve("37a1")
        sage: K.<t> = QuadraticField(-5)                                                                                                                   
        sage: EK = E.base_extend(K)
        sage: P = EK([0,0])                                                                                                                                             
        sage: Et= E.quadratic_twist(-5)                                                                                                                      
        sage: i = Et.base_extend(K).isomorphism_to(EK)
        sage: Q = i( Et.base_extend(K) ( Et.gens()[0] ) )
        sage: Q
        (3/4 : 1/8*t - 1/2 : 1)
        sage: reg([P,Q])
        0.28870082693041765657406710843325924521647963315306127987116

    """
    E = pts[0].curve()
    K = E.base_field()
    r = len(pts)
    return K.degree()**r * E.height_pairing_matrix(pts, precision=200).determinant()

def hh(a,b):    
    r"""
    Returns the canonical height pairing of two points on an elliptic curve over a number field.
    
    Rather than sage's usual normalisation, we make it dependent on the field here, so it
    is simply the degree times the pairing computed by sage
 
    EXAMPLES::
    
        sage: E = EllipticCurve("37a1")
        sage: K.<t> = QuadraticField(-5)                                                                                                                   
        sage: EK = E.base_extend(K)
        sage: P = EK([0,0])                                                                                                                                             
        sage: Et= E.quadratic_twist(-5)                                                                                                                      
        sage: i = Et.base_extend(K).isomorphism_to(EK)
        sage: Q = i( Et.base_extend(K) ( Et.gens()[0] ) )
        sage: hh(P,Q)
        -2.4892061111444566828576256215120496962361008674884668532404e-60
        sage: hh(Q,Q)
        2.8242308016144152840625208880105533709255350297682727455320
        sage: hh(P+Q,Q)
        2.8242308016144152840625208880105533709255350297682727455320

    """
    E = a.curve()
    K = E.base_field()
    return K.degree() * E.height_pairing_matrix([a,b], precision=200)[0][1]


def saturate(pts, ell):
    r""" 
    Given a list of points on an elliptic curve over a number field and a prime ell
    this returns a set of linearly independent point which generate the 
    ell-stauration of the group generated by the given points.
    
    EXAMPLES::
    
        sage: E = EllipticCurve("37a1")
        sage: K.<t> = QuadraticField(-5)                                                                                                                   
        sage: EK = E.base_extend(K)
        sage: P = EK([0,0])                                                                                                                                             
        sage: Et= E.quadratic_twist(-5)                                                                                                                      
        sage: i = Et.base_extend(K).isomorphism_to(EK)
        sage: Q = i( Et.base_extend(K) ( Et.gens()[0] ) )
        sage: saturate([P,Q],2)
        saturate at 2
        [(0 : 0 : 1), (3/4 : 1/8*t - 1/2 : 1)]
        sage: saturate([P,Q],7)
        saturate at 7
        [(0 : 0 : 1), (3/4 : 1/8*t - 1/2 : 1)]

        sage: R1 = 3*Q
        sage: R2 = 3*Q + P
        sage: saturate([R1,R2],3)
        saturate at 3
        (131363/362404 : 24846383/218167208*t - 1/2 : 1) is divisible by 3
        saturate at 3
        [(-14957522566/17256237769*t + 18061852956/17256237769 : -3455800971171206/2266831162049147*t - 2080701023629480/2266831162049147 : 1), (3/4 : 1/8*t - 1/2 : 1)]

    """
    print "saturate at %s"%(ell)
    E = pts[0].curve()
    # this is a basis of E_tors/ell
    Ts = [a[0] for a in E._p_primary_torsion_basis(ell)]
    Ba =  pts + Ts
    
    Q = Ba[0]
    if Q.is_divisible_by(ell):
        P1 = Q.division_points(ell)[0]
        print "%s is divisible by %s"%(Q,ell)
        tmp = pts.pop(pts.index(Q))
        pts += [P1]
        return saturate(pts, ell) 
   
    k = GF(ell)
    
    for i in range(1,len(Ba)):
        P = Ba[i]
        V = VectorSpace(k, i)
        for v in V:
            Q = P + sum( lift(v[j])*Ba[j] for j in range(i) ) 
            if Q.is_divisible_by(ell):
                P1 = Q.division_points(ell)[0]
                print "%s is divisible by %s"%(Q,ell)
                if P in pts:
                    tmp = pts.pop(pts.index(P))
                else:
                    j = 0
                    while v[j] == 0:
                        j += 1
                    tmp = pts.pop(pts.index(Ba[j]))
                pts += [P1]
                return saturate(pts,ell)
                    
    return pts

# =====================================
# now to the main class

class etnc_example(SageObject):
    r"""
    Returns a class that holds all information about the example on our way to compute them.
    """
    
    def __init__(self,E,f, prec = 20):
        r""" 
        Initialises the class 
        
        INPUT:
        
            - ``E`` - an elliptic curve over `\QQ` or just a Cremona label for it
            - ``f`` - a polynomial over `\QQ`
            - ``prec`` (optional) - the precision with which the L-values are computed, default 20
            
            
        The polynomial f must be of odd prime degree and its splitting field F must be a dihedral extension.
        
        This class has the following methods that can be evaluated in this order only:
        
        * check_conditions()
        * over_Q()
        * over_K()
        * eq_L_values()
        * over_L()
        * find_Q()
        * over_F()
        * eq_regs()
        * alg_L_values()
        * truncated_L_values()
        * check_congruence()

        For examples how to use them see the individual doc-string. The can also be evaluated at once by do_all()
        
        """
        if type(E) == type("11a1"):
            E = EllipticCurve(E)
            
        self.E = E
        self.prec = prec
        
        self.L = NumberField(f, "t")
        self.F = self.L.galois_closure("s")
        self.K = self.F.subfields(2)[0][0]
        self.EF = E.base_extend(self.F)
        self.EL = E.base_extend(self.L)
        self.EK = E.base_extend(self.K)
        
        self.i_LF = self.L.embeddings(self.F)[0]
        self.i_KF = self.K.embeddings(self.F)[0]
   
        self.p = self.F.degree() // 2
        assert self.p.is_prime(), "the extension must be prime dihedral"
        self.periods()
 
    def do_all(self, bound=30, known_point=None):
        r"""
        Performs all methods on this example to try to determine if the eTNC_p holds.
        They are in order check_conditions(), over_Q(), over_K(), eq_L_values(), over_L(), find_Q(),
        over_F(), eq_regs(), alg_L_values(), truncated_L_values(), check_congruence().

        See the individual methods for the docstrings. The optional arguments are passed on to find_points_over_L
        
        EXAMPLE::
        
            sage: R.<X> = QQ[]
            sage: d = etnc_example("17a1", X^3 - X^2 - 3*X + 1, prec=8)
            sage: d.do_all()
            CONDITIONS HOLD
            -- start over_Q 
            bsd over Q asserts 0.2500000 = 1/4 * Sha 
            -- start over_K 
            found point on twist (150455134/974169 : 1539419885296/961504803 : 1)
            gives point (4303315/974169 : 1178721325/17307086454*s0 + 1553922079/2472440922 : 1) over K=QQ(sqrt(37))
            saturate at 2
            bsd over K=Q(sqrt(37)) asserts 0.2500001 = 1/4 * Sha 
            -- start find_points_over_L 
            magma finds the point  (-1/2*t + 9/2 : 2*t - 21/2 : 1)
            saturate at 2
            saturate at 3
            saturate at 5
            saturate at 7
            this may not yet be a basis of E(L). We assume it however for the rest !
            -- start eq_L_values 
            the L-values twisted by 2-dimensional representations are  [5.209018]
            -- start over_L 
            bsd over L of discriminant 2^2 * 37 asserts 0.9999996 = 1 * Sha 
            -- start find_Q 
            the point Q was obtained from S = (1) and is equal to (290/11073*s^5 + 546/3691*s^4 - 9416/11073*s^3 - 16241/3691*s^2 + 23004/3691*s + 363788/11073 : -1853/14764*s^5 - 18773/44292*s^4 + 62866/11073*s^3 + 447467/22146*s^2 - 2104897/44292*s - 2495615/14764 : 1) 
            -- start over_F 
            saturate at 2
            saturate at 3
            bsd over F of discriminant 2^4 * 37^3 asserts 3.999998 = 4 * Sha 
            -- start eq_regs 
            the equivariant regulators are H_eins = 1, H_eps = 31.110050571768522086089555234021941844238064927019940608904 and H_rho = [6.6191394875629486406090365973490256838043869035639617464676]
            -- start alg_L_values 
            permutation found. Setting the trace to 4 with error 1.788139e-6
            the equivariant algebraic L-values are Qtilde_one = 1/4, Qtilde_eps = 1, Qtilde_rho has minimal polynomial  x - 4
            -- start truncated_L_values 
            the truncated equivariant L-values are Q_one = 20/37, Q_eps = 1, sum Q_rho = -160/37
            etnc_3 holds numerically for this example as 20/37 == 320/37 (mod 3)
                
        """       
        if self.check_conditions(verbose=False):
            print "CONDITIONS HOLD"
        self.over_Q()
        self.over_K()
        self.find_points_over_L(bound=bound,known_point=known_point)
        self.eq_L_values()
        self.over_L()
        self.find_Q()
        self.over_F()
        self.eq_regs()
        self.alg_L_values()
        self.truncated_L_values()
        self.check_congruence()

    def check_conditions(self, verbose=True):
        r"""
        Checks if the (rather restrictive) conditions in our paper hold for this example.
        The computations can still be performed, but it is not usual that the congruence
        will hold without the conditions here verified.
        
        This test all the conditions, except the finiteness of Sha and the condition that
        the analytic and algebraic rank is 1 over the quadratic field. The latter is done later in computations anyway.
        
        It optional argument `verbose` (True by default) allows to see which was the 
        first condition that did not hold.
        
        OUTPUT: a boolean.
        
        EXAMPLES::
        
            sage: R.<X> = QQ[]
            sage: d = etnc_example("17a1", X^3 - X^2 - 3*X + 1, prec=8)
            sage:  d.check_conditions()
            conditions for section 1 hold
            all conditions hold, except maybe (h) and (i)
            True
 
            sage: d = etnc_example("11a1", X^3 - X^2 - 3*X + 1, prec=8)
            sage: d.check_conditions()
            conditions for section 1 hold
            all conditions hold, except maybe (h) and (i)
            True

            sage: d = etnc_example("37a1", X^3 - X^2 - 3*X + 1, prec=8)
            sage: d.check_conditions()
            (e) There are bad places that ramifies 
            False

        """
        E = self.E
        EK = self.EK
        K = self.K
        F = self.F  
        p = self.p

        # conditions 
        # a
        if EK.torsion_order() % p == 0:
            if verbose: print "(a) E has %s-torsion over K=Q(sqrt(%s))"%(p,K.disc()) 
            return False

        A = F.disc().support()
        self.SrK = []
        for ell in A:
            vs = F.primes_above(ell)
            v = vs[0]
            if v.ramification_index() == p:
                self.SrK += K.primes_above(ell)
        da = EK.local_data()
        self.SbK = [dav.prime() for dav in da if not dav.has_good_reduction()]
        self.SpK = K.primes_above(p)
        # if verbose: print "S_b^K = %s, S_r^K = %s, S_p^K = %s"%(self.SbK,self.SrK,self.SpK)

        # b
        for v in self.SbK : 
            dav = EK.local_data(v)
            if dav.tamagawa_number() % p == 0:
                if verbose: print "(b) %s divides the Tamagawa number %s at %s"%(p,dav.tamagawa_number(),v)
                return False
 
        # c
        if E.discriminant() % p == 0:
            if verbose: print "(c) E does not have good reduction at %s"%(p)
            return False

        # d
        for v in self.SpK:
            if v in self.SrK:
                if not E.is_ordinary(p):
                    if verbose: print "(d) E does not have ordinary reduction at %s at which it ramifies in F/K"%(p)
                    return False
                if generalised_Nv(EK, v) % p == 0:
                    if verbose: print "(d) E does not have non-anomalous reduction at %s at which it ramifies in F/K"%(v)
                    return False

        # e
        if gcd( F.disc(), E.conductor() ) > 1:
            if verbose: print "(e) There are bad places that ramifies "%(gcd( F.disc(), E.conductor() ).support()) 
            return False

        if verbose: print "conditions for section 1 hold"

        # f
        if p == 2:
            if verbose: print "(f) p = 2"
            return False  

        # g
        for ell in F.disc().support():
            for v in K.primes_above(ell):
                if generalised_Nv(EK,v) % p == 0:
                    if verbose: print "(g) E has %s-torsion in the reduction at %s"%(p,v)
                    return False

        # h 
        # we assume shas are finite

        # i
        # L-functions extend for dihedral extensions.
        # we assume the order of vanishing will work out well later

        # j
        if F.disc() % p == 0:
            if verbose: print "(j) p=%s is ramified in F"%p
            return False

        # addtional condition in Thm 4.1
        if len(K.primes_above(p)) < 2:
            if verbose: "(Thm 4.1) %s should split in K"%p
            return False

        if verbose: print "all conditions hold, except maybe (h) and (i)"

        return True

    def periods(self):
        r"""
        Set the periods of the elliptic curve. The normalisation here is that the plus-period is the 
        integral around the generator gamma+ for the Neron differential given by the globale 
        minimal model. (So the number of connected components is in the value Q_rho).
        
        This function sets Om, Om_F, Om_K, Om_L and Om_eps, Om_rho (for any dim(rho)=2).
        This method is evaluated at __init__
        
        EXAMPLES::
        
            sage: R.<X> = QQ[]
            sage: d = etnc_example("11a1", X^3 - X^2 - 3*X + 1, prec=8)
            sage: d.Om
            1.2692093042795534216887946167545473052194922418306086679671
            sage: d = etnc_example("11a1", X^3 - 7, prec=8)
            sage: d.Om_F
            50.779899055235133494420766435144855246912416178394265718078

        """
        p = self.p
        E = self.E
        self.cinf = E.real_components()
        Ba = E.period_lattice().basis(200)
        self.Om = Ba[0]
        self.Om_minus = Ba[1].imag()*2/self.cinf
        
        self.Om_RR = self.Om
        self.Om_CC =  self.Om * self.Om_minus
        if self.K.disc() < 0:
            self.Om_K = self.Om_CC
            self.Om_eps = self.Om_CC/self.Om_RR
            self.Om_L = self.Om_RR * (self.Om_CC) ** ( (p-1)//2 )
            self.Om_rho = (self.Om_CC) ** ( (p-1)//2)
            self.Om_F = (self.Om_CC ) ** p
        else:
            self.Om_K = self.Om_RR ** 2
            self.Om_eps = self.Om_RR
            self.Om_L = self.Om_RR ** p
            self.Om_rho =  self.Om_RR ** ( p - 1)
            self.Om_F = self.Om_RR ** (2* p)
             
    def over_Q(self):
        r"""
        This method checks if the BSD holds over `\QQ`.
        It sets the values for reg_Q, basis_Q and l_Q,
        which is the leading term of the L-series.
        
        EXAMPLES::

            sage R.<X> = QQ[]
            sage: d = etnc_example("11a1", X^3 - 7, prec=8)
            sage: d.over_Q()
            -- start over_Q 
            bsd over Q asserts 0.2000 = 1/5 * Sha 
            
            sage: d = etnc_example("37a1", X^3 +2*X - 1, prec=32)
            sage: d.over_Q()
            -- start over_Q 
            bsd over Q asserts 2.00000000000000000 = 2 * Sha 
        
        """
        if hasattr(self,"l_Q"):
            print "Already done."
        else:
            print "-- start over_Q "
            E = self.E
            assert E.rank() < 2
            if E.rank() == 1:
                self.rkQ = 1
            else:
                self.rkQ = 0
 
            if self.rkQ == 1:
                self.R_Q = E.gens()[0]
                self.basis_Q = [self.R_Q]
                self.reg_Q = nt_height(self.R_Q)
            else:
                self.basis_Q = []
                self.reg_Q = 1

            Em = magma(E)
            LEm = Em.LSeries(Precision = self.prec)
            lEQ = LEm.Evaluate(1, Derivative = self.rkQ)
            self.l_Q = RealField(3.32*self.prec)( lEQ )

            # bsd
            rhs = self.cinf * E.tamagawa_product()/E.torsion_order()**2 
            lhs = self.l_Q/self.Om_RR/self.reg_Q
            print "bsd over Q asserts %s = %s * Sha "%(lhs, rhs)

    def over_K(self):
        r"""
        This method checks if bsd holds over K. The method sets
        reg_K, basis_K and l_K. An error occurs if the rank is 
        not 1.
        
        EXAMPLE::
        
            sage: R.<X> = QQ[]
            sage: d = etnc_example("37a1", X^3 +2*X - 5, prec=8)
            sage: d.over_K()
            -- start over_Q 
            bsd over Q asserts 2.000000 = 2 * Sha 
            -- start over_K 
            saturate at 2
            bsd over K=Q(sqrt(-707)) asserts 50.00001 = 2 * Sha 

        """
        if hasattr(self,"l_K"):
            print "Already done."
        else:
            if not hasattr(self,"l_Q"):
                self.over_Q()
            print "-- start over_K "
            EK = self.EK
            K = self.K

            if self.rkQ == 1:
                R_K = EK( self.R_Q )
                
            else:
                Et = self.E.quadratic_twist(K.disc())
                Et.two_descent(verbose=False,second_limit=13)
                assert Et.rank() == 1, "the rank over K has to be 1"
                Pt = Et.gens()[0]
                print "found point on twist %s"%Pt     
                Et_K = Et.base_extend(K)
                Pt_K = Et_K(Pt)
                iso_K = Et_K.isomorphism_to(EK)
                R_K = iso_K(Pt_K)
                print "gives point %s over K=QQ(sqrt(%s))"%(R_K,K.disc())
        
            pts = saturate([R_K],2)
            R1_K = pts[0]
            # now R1_K is a basis of the free part of E(K)
            self.basis_K = [R1_K]       
            self.reg_K = nt_height(R1_K)     

            EKm = magma(EK)
            LEKm = EKm.LSeries(Precision=self.prec)
            lEK = LEKm.Evaluate(1, Derivative = 1)
            self.l_K = RealField(3.32*self.prec)( lEK )

            assert self.l_K > 0.0001, "the rank must be 1 over K"

            tam = EK.tamagawa_product_bsd()
            tor = EK.torsion_order()
         
            # bsd
            rhs = self.cinf * tam / tor**2 
            if K.disc() > 0 :
                rhs *= self.cinf 
            sqK =  sqrt(abs(RealField(5*self.prec)(K.disc())))
            lhs = self.l_K/self.Om_K/self.reg_K * sqK
            print "bsd over K=Q(sqrt(%s)) asserts %s = %s * Sha "%(K.disc(),lhs, rhs)

    def galois_group(self):
        r"""
        This method identifies the Galois group. It sets
        sigma and tau, an element of order p and 1 repsectively.
        tau fixes L, sigma fixes K.
        
        EXAMPLES::
        
            sage: R.<X> = QQ[]
            sage: d = etnc_example("37a1", X^3 +2*X - 14, prec=8)
            sage: d.galois_group()
            sage: d.sigma
            Ring endomorphism of Number Field in s with defining polynomial X^6 + 28*X^4 + 280*X^3 + 196*X^2 + 3920*X + 67516
            Defn: s |--> -35/454589*s^5 + 7203/1818356*s^4 - 2450/1363767*s^3 + 64435/909178*s^2 + 127903/2727534*s + 196196/1363767
            sage: d.tau
            Ring endomorphism of Number Field in s with defining polynomial X^6 + 28*X^4 + 280*X^3 + 196*X^2 + 3920*X + 67516
            Defn: s |--> 40/454589*s^5 - 2058/454589*s^4 + 2800/1363767*s^3 - 36820/454589*s^2 - 657559/1363767*s - 224224/1363767
            
        """
        if hasattr(self,"tau"):
            print "Already done"
        else:
            # print "-- start galois_group "
            F = self.F            
            s = F.gens()[0]
            i_K = self.i_KF
            i_L = self.i_LF
            t = i_L(self.L.gens()[0])
            G = F.automorphisms()
            p = self.p
            i = 0
            found = False
            while not found and i < len(G):
                si = G[i]
                if si(s) != s and (si**p)(s) == s:
                    found = True
                i += 1
            assert i < len(G)
            self.sigma = si
            i = 0
            found = False
            while not found and i < len(G):
                tau = G[i]
                if tau(s) != s and tau(t) == t:
                    found = True
                i += 1
            self.tau = tau

    def sigma_action(self,pt):
        r"""
        Given a point P, this returns sigma(P)
        where sigma is the chosen generator of
        the cyclic group of order p.
        
        needs: galois_group()
        
        EXAMPLES::
    
            sage: R.<X> = QQ[]
            sage: f = X**5-X**4-5*X**3+4*X**2+3*X-1 
            sage: d = etnc_example("37a1",f, prec=16)
            sage: RF.<x> = d.F[]
            sage: f_x = x^5 - 73/9*x^4 + 1195/81*x^3 - 173/81*x^2 - 976/81*x + 527/81
            sage: x_P = f_x.roots()[0][0]
            sage: P = d.EF.lift_x(x_P)
            sage: d.galois_group()
            sage: d.sigma_action(P)
            (23232417241/398771984143746*s^9 - 73869901/312028156607*s^8 - 4111213480669/1196315952431238*s^7 + 3559820814049/199385992071873*s^6 + 18009507951367/398771984143746*s^5 - 456684363436381/1196315952431238*s^4 - 35568895654099/1196315952431238*s^3 + 954007839815915/398771984143746*s^2 - 578463490542373/398771984143746*s - 1102689294141253/1196315952431238 : 965696443978/1794473928646857*s^9 + 6416085887/25274280685167*s^8 - 121425690631819/3588947857293714*s^7 - 4837613638381/1794473928646857*s^6 + 2016283273691699/3588947857293714*s^5 - 58916173878437/199385992071873*s^4 - 5086849671834962/1794473928646857*s^3 + 5795968505566841/1794473928646857*s^2 + 220137140496889/199385992071873*s - 9132050569326803/3588947857293714 : 1)

        """
        return self.EF( [self.sigma(pt[0]), self.sigma(pt[1])])
        
    def sigma_inverse_action(self,pt):
        r"""
        Given a point P, this returns sigma^(-1)(P)
        where sigma is the chosen generator of
        the cyclic group of order p.
        
        needs: galois_group()
        
        EXAMPLES::
    
            sage: R.<X> = QQ[]
            sage: f = X**5-X**4-5*X**3+4*X**2+3*X-1 
            sage: d = etnc_example("37a1",f, prec=16)
            sage: RF.<x> = d.F[]
            sage: f_x = x^5 - 73/9*x^4 + 1195/81*x^3 - 173/81*x^2 - 976/81*x + 527/81
            sage: x_P = f_x.roots()[0][0]
            sage: P = d.EF.lift_x(x_P)
            sage: d.galois_group()
            sage: d.sigma_inverse_action(P)
            (-3960449405225/70582641193443042*s^9 + 18501908533/36819322479626*s^8 + 45199356076355/7842515688160338*s^7 - 2025303761026435/70582641193443042*s^6 - 4223361986063161/23527547064481014*s^5 + 26648495578651417/70582641193443042*s^4 + 41216083599311711/23527547064481014*s^3 - 107454130029402661/70582641193443042*s^2 - 150809937566777951/35291320596721521*s + 272505635371808035/70582641193443042 : -31987064989079/211747923580329126*s^9 + 2004460504226/1491182560424853*s^8 + 1684571288272697/105873961790164563*s^7 - 8152458293095291/105873961790164563*s^6 - 106908375981280783/211747923580329126*s^5 + 108901951811472802/105873961790164563*s^4 + 523383512849808038/105873961790164563*s^3 - 912810860458624721/211747923580329126*s^2 - 1252946408086014796/105873961790164563*s + 551526472907812747/70582641193443042 : 1)

        """
        return self.EF( [((self.sigma)**(-1))(pt[0]), ((self.sigma)**(-1))(pt[1])])
    
    def tau_action(self,pt):
        r"""
        Given a point P, this returns tau(P)
        where tau is an element of order 2
        which fixed the chosen field L.
        
        needs: galois_group()
        
        EXAMPLES::
    
            sage: R.<X> = QQ[]
            sage: f = X**5-X**4-5*X**3+4*X**2+3*X-1 
            sage: d = etnc_example("37a1",f, prec=16)
            sage: RF.<x> = d.F[]
            sage: f_x = x^5 - 73/9*x^4 + 1195/81*x^3 - 173/81*x^2 - 976/81*x + 527/81
            sage: x_P = f_x.roots()[0][0]
            sage: P = d.EF.lift_x(x_P)
            sage: d.galois_group()
            sage: d.tau_action(P) == P
            True
            sage: siP = d.sigma_action(P)
            sage: d.tau_action(siP)
            (-3960449405225/70582641193443042*s^9 + 18501908533/36819322479626*s^8 + 45199356076355/7842515688160338*s^7 - 2025303761026435/70582641193443042*s^6 - 4223361986063161/23527547064481014*s^5 + 26648495578651417/70582641193443042*s^4 + 41216083599311711/23527547064481014*s^3 - 107454130029402661/70582641193443042*s^2 - 150809937566777951/35291320596721521*s + 272505635371808035/70582641193443042 : -31987064989079/211747923580329126*s^9 + 2004460504226/1491182560424853*s^8 + 1684571288272697/105873961790164563*s^7 - 8152458293095291/105873961790164563*s^6 - 106908375981280783/211747923580329126*s^5 + 108901951811472802/105873961790164563*s^4 + 523383512849808038/105873961790164563*s^3 - 912810860458624721/211747923580329126*s^2 - 1252946408086014796/105873961790164563*s + 551526472907812747/70582641193443042 : 1)
            
        """
        return self.EF( [self.tau(pt[0]), self.tau(pt[1])])

    def find_points_over_L(self, bound = 100, known_point = None, dont_saturate = False):
        r"""
        This method tries to find a point in E(L) that is not defined over Q. If found
        the Mordell-Weil group is saturated at a few primes.
        
        INPUT (all optional):
        
            - `bound` - Bound to give to the magma function Points. Default=30.
            - `known_point` - Can be a point in E(L) that is known from another source. Default=None
            - `dont_saturate` - By default (False) the group is saturated in the end.
                 
        OUTPUT: a list of points in E(L)
        
        A warning is issued that we do not guarantee that the full saturation of E(L) is found.
        
        EXAMPLES::
        
            sage: d = etnc_example("43a1",X^3-4*X+2,prec=16)
            sage: d.find_points_over_L()
            - start over_Q 
            bsd over Q asserts 1.00000000000000 = 1 * Sha 
            -- start find_points_over_L 
            magma finds the point  (-5*t^2 - 3*t + 18 : 19*t^2 + 10*t - 71 : 1)
            saturate at 2
            (11/49 : -363/343 : 1) is divisible by 2
            saturate at 2
            (2 : -4 : 1) is divisible by 2
            saturate at 2
            (-1 : 0 : 1) is divisible by 2
            saturate at 2
            saturate at 3
            saturate at 5
            saturate at 7
            this may not yet be a basis of E(L). We assume it however for the rest !
            [(-5*t^2 - 3*t + 18 : 19*t^2 + 10*t - 71 : 1), (0 : -1 : 1)]

            sage: E = EllipticCurve([-4,5])
            sage: d = etnc_example(E,X^3-4*X+1,prec=16)
            sage: t = d.L.gens()[0]
            sage: P = d.EL.lift_x(t)
            sage: P.order()
            +Infinity
            sage: d.find_points_over_L(known_point=P)
            -- start over_Q 
            bsd over Q asserts 1.00000000000000 = 1 * Sha 
            -- start find_points_over_L 
            saturate at 2
            saturate at 3
            saturate at 5
            saturate at 7
            this may not yet be a basis of E(L). We assume it however for the rest !
            [(t : 2 : 1)]

        """
        if hasattr(self, "_points_over_L"):
            print "Points already found."
            return self._points_over_L
        else:
            if not hasattr(self,"l_Q"):
                self.over_Q()
            print "-- start find_points_over_L "
            EL = self.EL
            L = self.L
            if self.rkQ == 1 :
                R_L = EL( self.R_Q )

            # find a single new point   
            def is_new(P):
                if self.rkQ == 1:
                    return  P.order() == +Infinity and abs( hh(R_L, P) ) > 0.00001
                else:
                    return P.order() == +Infinity
                    
            found = False    
            if known_point != None:
                P = known_point
                if P in EL and is_new(P):
                    found = True
            
            if not found:
                # use magma's Points
                ELm = magma(EL)
                pts = ELm.Points(Bound = bound)
                i = 1
                while not found and i <= len(pts):
                    P = pts[i]
                    xP = P[1]
                    xP = xP.sage()
                    if xP in L and xP not in QQ:
                        P = EL.lift_x(xP)
                        if is_new(P):
                            print "magma finds the point ", P       
                            found = True
                    i += 1

            #use simon BUGGY
            #if not found:
                #si = EL.simon_two_descent()
                #pts = si[2]
                #i = 0
                #while not found and i < len(pts):
                    #P = pts[i]
                    #xP = P[0]
                    #if xP in L and xP not in QQ:
                        #if is_new(P):
                            #print "simon's 2-descent finds the point ", P       
                            #found = True
                    #i += 1
    
            if found:
                # now get a full basis of E(L)
                self.galois_group()
                EF = self.EF
                pts = [P]
                i_L = self.i_LF
                P_F = EF([ i_L(P[0]), i_L(P[1]) ])
                siP = self.sigma_action(P_F)
                siiP = self.sigma_inverse_action(P_F)
                for i in range(1, (self.p-3)//2 + self.rkQ + 1):
                    R = siP + siiP # should be in EL
                    RL = EL( [i_L.preimage(R[0]),i_L.preimage(R[1]) ] )
                    pts.append( RL )
                    siP = self.sigma_action(siP)
                    siiP = self.sigma_inverse_action(siiP)
                
                if not dont_saturate:
                    for ell in union([2,3,5,7,self.p]):
                        pts = saturate(pts,ell)
                    # now we should use the Siksek bound to make sure.
                
                print "this may not yet be a basis of E(L). We assume it however for the rest !"
                self._points_over_L = pts
                return pts
            else:
                print "no new point found, you may try with a larger bound."
                return []
          
    def eq_L_values(self, prec = None):
        r"""
        This method evaluates the twisted L-values, using Tim's magma implementation.
         
        The optional argument ``prec`` allows to choose a different precision for these computations.
        
        EXAMPLE::
            
            sage: R.<X> = QQ[]
            sage: d = etnc_example("40a1",X^3-4*X+1,prec= 8)
            sage: d.eq_L_values()
            -- start eq_L_values 
            the L-values twisted by 2-dimensional representations are  [8.266613]
    
        """
        if hasattr(self,"l_rhos"):
            print "Already done."
        else:
            print "-- start eq_L_values "
            if prec == None:
                prec = self.prec
            F = self.F
            E = self.E
            Em = magma(E)
            Fm = magma(F)
            Am = Fm.ArtinRepresentations()
            Am = [a for a in Am][2:]
            ls = [Em.LSeries(a, Precision = prec) for a in Am]
            # magma prec is in number of decimal digits, so we want at least 3.32*prec
            vals = [RealField(3.32 * prec)( l.Evaluate(1, Derivative = 1) ) for l in ls]
            print "the L-values twisted by 2-dimensional representations are ", vals
            self.l_rhos = vals
            
    def over_L(self, dont_saturate = False):
        r"""
        Tests BSD over L, a non-Galois extension of degree p in F.
        It sets the values of basis_L, reg_L and l_L
          
        EXAMPLES::
        
            sage: R.<X> = QQ[]
            sage: d = etnc_example("17a1",X^3-X^2-3*X+1,prec=8)
            sage: d.over_L()
            -- start eq_L_values 
            the L-values twisted by 2-dimensional representations are  [5.209018]
            -- start over_Q 
            bsd over Q asserts 0.2500000 = 1/4 * Sha 
            -- start find_points_over_L 
            magma finds the point  (-1/2*t + 9/2 : 2*t - 21/2 : 1)
            saturate at 2
            saturate at 3
            saturate at 5
            saturate at 7
            this may not yet be a basis of E(L). We assume it however for the rest !
            -- start over_L 
            bsd over L of discriminant 2^2 * 37 asserts 0.9999996 = 1 * Sha 

        """
        if hasattr(self,"l_L"):
            print "Already done."
        else:
            if not hasattr(self,"l_rhos"):
                self.eq_L_values()
            if not hasattr(self,"l_Q"):
                self.over_Q()
            # find a basis of E(L)
            if not hasattr(self,"_points_over_L"):
                pts = self.find_points_over_L(dont_saturate=dont_saturate)
            else:
                pts = self._points_over_L
            print "-- start over_L "
            
            EL = self.EL
            L = self.L
            p = self.p
            
            rkL = (p-1)//2 + self.rkQ
            assert len(pts) == rkL, "not enough points over L, yet"
            self.basis_L = pts
            self.reg_L = reg(pts)

            #ELm = magma(EL)
            #LELm = ELm.LSeries(Precision = self.prec)
            #lEL = LELm.Evaluate(1, Derivative = rkL )
            self.l_L = prod(self.l_rhos)*self.l_Q

            assert self.l_L > 0.0001

            tam = EL.tamagawa_product_bsd()
            tor = EL.torsion_order()

            # bsd
            sqL =  sqrt(abs(RealField(5*self.prec)(L.disc())))
            r1 = len(L.real_places())
            rhs = self.cinf**r1 * tam / tor**2 
            lhs = self.l_L/self.Om_L/self.reg_L * sqL
            print "bsd over L of discriminant %s asserts %s = %s * Sha "%(L.disc().factor(),lhs, rhs)

    def find_Q(self):
        r"""
        Under our assumptions, there exists a point Q in E(F) such that Z[G] Q generates a finite
        prime-to-p index subgroup of E(F) and tau(Q) = +- Q depending on the rank of E over Q.
        
        This method finds such a point Q. The algorithm is very different for the rank = 1 and 
        rank = 0 cases. This is not fully tested but worked well so far.
        
        It sets Q
                
        EXAMPLES::
        
            sage: R.<X> = QQ[]
            sage: d = etnc_example("17a1",X^3-X^2-3*X+1,prec=8)
            sage: d.find_Q()
            -- start eq_L_values 
            the L-values twisted by 2-dimensional representations are  [5.209018]
            -- start over_Q 
            bsd over Q asserts 0.2500000 = 1/4 * Sha 
            -- start find_points_over_L 
            magma finds the point  (-1/2*t + 9/2 : 2*t - 21/2 : 1)
            saturate at 2
            saturate at 3
            saturate at 5
            saturate at 7
            this may not yet be a basis of E(L). We assume it however for the rest !
            -- start over_L 
            bsd over L of discriminant 2^2 * 37 asserts 0.9999996 = 1 * Sha 
            -- start over_K 
            found point on twist (150455134/974169 : 1539419885296/961504803 : 1)
            gives point (4303315/974169 : 1178721325/17307086454*s0 + 1553922079/2472440922 : 1) over K=QQ(sqrt(37))
            saturate at 2
            bsd over K=Q(sqrt(37)) asserts 0.2500001 = 1/4 * Sha 
            -- start find_Q 
            the point Q was obtained from S = (1) and is equal to (290/11073*s^5 + 546/3691*s^4 - 9416/11073*s^3 - 16241/3691*s^2 + 23004/3691*s + 363788/11073 : -1853/14764*s^5 - 18773/44292*s^4 + 62866/11073*s^3 + 447467/22146*s^2 - 2104897/44292*s - 2495615/14764 : 1)
            
        """
        if hasattr(self,"Q"):
            print "Already done."
        else:  
            if not hasattr(self,"l_L"):
                self.over_L()
            if not hasattr(self,"l_K"):
                self.over_K()
            print "-- start find_Q "
            
            EF = self.EF
            p = self.p
            si = self.sigma_action
            i_K = self.i_KF
            i_L = self.i_LF
            
            P = self.basis_L[0]
            P_F = EF([ i_L(P[0]), i_L(P[1]) ])
            R_K = self.basis_K[0]
            R_F = EF([ i_K(R_K[0]), i_K(R_K[1]) ])
            self.R_F = R_F
            
            if self.rkQ == 1:
                # computes the norm of P
                NP = P_F
                sP = P_F
                for i in range(1,p):
                    sP = si(sP)
                    NP += sP
            
                # find Q now.
                n = hh(NP,R_F)/hh(R_F,R_F)
                n_r = round(n)
                assert abs(n-n_r) < 0.0001, "failed to determine the norm of P correctly"
                n = ZZ(n_r)
                print "norm of P is %s * R"%n
                g, a,b = xgcd(n, p)
                assert g == 1    # otherwise we did not saturate well
                self.Q = a * P_F + b * R_F
                print "the point Q is %s * P + %s * R = %s"%(a,b, self.Q)
            
            else : # rkQ = 0
            
                # this uses the formula that N + (si-1) * sum( i*(si^i-si^-i) for i in [1.. (p-1)/2] ) = p * si ^( (p+1)/2 )
                # the point T below is the generator of the kernel of A/p -> E(F)/p where A is the group generated by E(K), E(L) and si(E(L)) such that its norm to K is R. We then know that T = p * si^( - (p-1)/2 )(Q)
                # S will be the sum( i* (si^i - si^(-i))(Q) ) above
            
                # building a basis of E(L) of P, (si+si^(-1)) P , ...
                baL = self._points_over_L
                baL = [ EF([ i_L(sP[0]), i_L(sP[1]) ]) for sP in baL ]
                    
                V = VectorSpace(GF(p),(p-1)//2)
                for v in V:
                    if v != V(0):
                        S = sum( ZZ(v[i])*baL[i] for i in range((p-1)//2) )
                        T = R_F - S + si(S)
                        if T.is_divisible_by(p):
                            break
                
                # now T is divisible by p
                Q = T.division_points(p)[0]
                for i in range( (p-1)//2 ):
                    Q = si(Q)
                self.Q = Q
                
                print "the point Q was obtained from S = %s and is equal to %s "%(v, self.Q)
  
    def over_F(self, dont_saturate = False):
        r"""
        tests BSD over F. It sets reg_F, basis_F and l_F.
                
        EXAMPLE::          
        
            sage: R.<X> = QQ[]
            sage: d = etnc_example("17a1",X^3-X^2-3*X+1,prec=8)
            sage: d.over_F()
            -- start eq_L_values 
            the L-values twisted by 2-dimensional representations are  [5.209018]
            -- start over_Q 
            bsd over Q asserts 0.2500000 = 1/4 * Sha 
            -- start find_points_over_L 
            magma finds the point  (-1/2*t + 9/2 : 2*t - 21/2 : 1)
            saturate at 2
            saturate at 3
            saturate at 5
            saturate at 7
            this may not yet be a basis of E(L). We assume it however for the rest !
            -- start over_L 
            bsd over L of discriminant 2^2 * 37 asserts 0.9999996 = 1 * Sha 
            -- start over_K 
            found point on twist (150455134/974169 : 1539419885296/961504803 : 1)
            gives point (4303315/974169 : 1178721325/17307086454*s0 + 1553922079/2472440922 : 1) over K=QQ(sqrt(37))
            saturate at 2
            bsd over K=Q(sqrt(37)) asserts 0.2500001 = 1/4 * Sha 
            -- start find_Q 
            the point Q was obtained from S = (1) and is equal to (290/11073*s^5 + 546/3691*s^4 - 9416/11073*s^3 - 16241/3691*s^2 + 23004/3691*s + 363788/11073 : -1853/14764*s^5 - 18773/44292*s^4 + 62866/11073*s^3 + 447467/22146*s^2 - 2104897/44292*s - 2495615/14764 : 1) 
            -- start over_F 
            saturate at 2
            saturate at 3
            bsd over F of discriminant 2^4 * 37^3 asserts 3.999998 = 4 * Sha 
            
        """
        if hasattr(self,"l_F"):
            print "Already done."
        else:
            if not hasattr(self,"Q"):
                self.find_Q()
            print "-- start over_F "

            EF = self.EF
            F = self.F
            Q = self.Q
            p = self.p
            si = self.sigma_action
            
            ba = [Q]
            sQ = Q
            for i in range(1,p):
                sQ = si(sQ)
                ba.append(sQ)
            
            assert reg(ba) > 0.0001, "Q can not be a generator !!"

            if not dont_saturate:
                for ell in union([2,3,self.p]):
                    pts = saturate(ba, ell)
                if pts != ba:
                    "points had to be saturated. Careful!"
            else:
                pts = ba

            self.basis_F = pts
            self.reg_F = reg(pts)

            self.l_F = self.l_L**2  * self.l_K / self.l_Q**2

            tam = EF.tamagawa_product_bsd()
            tor = EF.torsion_order()
            
            # bsd
            sqF =  sqrt(abs(RealField(5*self.prec)(F.disc())))
            rhs = tam / tor**2 
            r1 = len(F.real_places())
            # r2 = (F.degree() - r1 ) // 2
            rhs *= self.cinf**r1
            lhs = self.l_F/self.Om_F /self.reg_F * sqF
            print "bsd over F of discriminant %s asserts %s = %s * Sha "%(F.disc().factor(),lhs, rhs)
  
    def eq_regs(self):
        r"""
        This method evaluates the equivariant regulators. These are the numbers H_rho appearing in the
        denominator of Q. (They are not normalised to satisfy the Artin formalism here.) 
        
        It sets h_eins, h_eps and h_rhos
         
        EXAMPLE::
        
            sage: R.<X> = QQ[]
            sage: d = etnc_example("55a1",X^3-X^2-3*X+1,prec=8)
            sage: d.eq_regs()
            -- start eq_L_values 
            the L-values twisted by 2-dimensional representations are  [10.36345]
            -- start over_Q 
            bsd over Q asserts 0.5000000 = 1/2 * Sha 
            -- start find_points_over_L 
            magma finds the point  (1/4*t^2 + 1/2*t + 1/4 : 1/4*t^2 + 1/2*t - 3/2 : 1)
            saturate at 2
            saturate at 3
            saturate at 5
            saturate at 7
            this may not yet be a basis of E(L). We assume it however for the rest !
            -- start over_L 
            bsd over L of discriminant 2^2 * 37 asserts 4.000003 = 4 * Sha 
            -- start over_K 
            found point on twist (16 : 107 : 1)
            gives point (25/37 : 77/16428*s0 - 1777/16428 : 1) over K=QQ(sqrt(37))
            saturate at 2
            (19/9 : -7/972*s0 - 1369/972 : 1) is divisible by 2
            saturate at 2
            bsd over K=Q(sqrt(37)) asserts 2.000000 = 2 * Sha 
            -- start find_Q 
            the point Q was obtained from S = (2) and is equal to (589/276825*s^5 + 3007/184550*s^4 + 51619/1107300*s^3 - 52213/1107300*s^2 - 1285117/1107300*s + 188003/1107300 : -4477/276825*s^5 - 387913/2768250*s^4 + 35916/461375*s^3 + 3298582/1384125*s^2 + 1358534/1384125*s - 21816331/2768250 : 1) 
            -- start eq_regs 
            the equivariant regulators are H_eins = 1, H_eps = 1.5271486849022197351774638047526182962134576410193697226404 and H_rho = [14.893361251872325072715266360904169325938941176678794029968]

        """       
        if hasattr(self,"h_eps"):
            print "Already done."
        else:            
            if not hasattr(self,"Q"):
                self.find_Q()
            print "-- start eq_regs "

            # compute equivariant regulators
            Q = self.Q
            EF = self.EF
            p = self.p
            si = self.sigma_action
            assert Q in EF
            
            if self.rkQ == 0:
                self.h_eins = 1
                self.h_eps = self.reg_K
            else:
                self.h_eins = 2 * self.reg_Q
                self.h_eps = 1
            
            self.h_rhos = []
            pir = RealField(200)(pi)
            for j in range(1, (p-1)//2+1):
                hrhoj = 0
                sQ = Q
                for k in range(0,p):
                    hrhoj += cos(2*pir*j*k/p)*hh(Q,sQ)
                    sQ = si(sQ)    
                self.h_rhos.append(hrhoj)
            
            print "the equivariant regulators are H_eins = %s, H_eps = %s and H_rho = %s"%(self.h_eins, self.h_eps, self.h_rhos)
 
    def alg_L_values(self, bound = None):   
        r"""
        Compute the algebraic \tilde{Q}_rho for all rho.
        
        There is a bit of guessing in this as the order of the L-values and the order
        of the equivariant regulators may differ. We test all possibility and take one
        where the trace of the Qtildes is sufficiently close to a rational number of
        small denominator.
        
        It sets l_rhos, Qtilde_eins, Qtilde_eps, Qtilde_rhos, Qtilde_trace_rhos, Qtilde_norm_rhos
              
        The optional argument bound can be set as the bound to the denominator in the algebraic values.
        
        EXAMPLE::
        
            sage: R.<X> = QQ[]
            sage: d = etnc_example("17a1",X^3-X^2-3*X+1,prec=8)
            sage: d.alg_L_values()
            -- start eq_L_values 
            the L-values twisted by 2-dimensional representations are  [5.209018]
            -- start over_Q 
            bsd over Q asserts 0.2500000 = 1/4 * Sha 
            -- start find_points_over_L 
            magma finds the point  (-1/2*t + 9/2 : 2*t - 21/2 : 1)
            saturate at 2
            saturate at 3
            saturate at 5
            saturate at 7
            this may not yet be a basis of E(L). We assume it however for the rest !
            -- start over_L 
            bsd over L of discriminant 2^2 * 37 asserts 0.9999996 = 1 * Sha 
            -- start over_K 
            found point on twist (150455134/974169 : 1539419885296/961504803 : 1)
            gives point (4303315/974169 : 1178721325/17307086454*s0 + 1553922079/2472440922 : 1) over K=QQ(sqrt(37))
            saturate at 2
            bsd over K=Q(sqrt(37)) asserts 0.2500001 = 1/4 * Sha 
            -- start find_Q 
            the point Q was obtained from S = (1) and is equal to (290/11073*s^5 + 546/3691*s^4 - 9416/11073*s^3 - 16241/3691*s^2 + 23004/3691*s + 363788/11073 : -1853/14764*s^5 - 18773/44292*s^4 + 62866/11073*s^3 + 447467/22146*s^2 - 2104897/44292*s - 2495615/14764 : 1) 
            -- start eq_regs 
            the equivariant regulators are H_eins = 1, H_eps = 31.110050571768522086089555234021941844238064927019940608904 and H_rho = [6.6191394875629486406090365973490256838043869035639617464676]
            -- start alg_L_values 
            permutation found. Setting the trace to 4 with error 1.788139e-6
            the equivariant algebraic L-values are Qtilde_one = 1/4, Qtilde_eps = 1, Qtilde_rho has minimal polynomial  x - 4

            """           
        if hasattr(self,"Qtilde_eps"):
            print "Already done."
        else:            
            if not hasattr(self,"h_eps"):
                self.eq_regs()
            print "-- start alg_L_values "

            K = self.K
            F = self.F
            p = self.p
            vals = self.l_rhos
            
            d_eps = sqrt(abs(RealField(5*self.prec)(K.disc())))
            sqNfphi = F.disc() / K.disc()**p
            sqNfphi = sqNfphi.nth_root(p-1)
            d_rhos = sqrt(abs(RealField(5*self.prec)(sqNfphi))) * d_eps
                
            l_eins = self.l_Q/self.Om/self.h_eins
            l_eps = self.l_K/ self.Om_eps/self.l_Q * d_eps / self.h_eps
            tor = self.E.torsion_order()
            if bound == None:
                bo = tor**2 * 10
            else:
                bo = bound
            self.Qtilde_eins = l_eins.nearby_rational(max_denominator=bo)
            self.Qtilde_eps = l_eps.nearby_rational(max_denominator=bo)
            assert abs(l_eins-self.Qtilde_eins) < 0.0001 and abs(l_eps-self.Qtilde_eps) < 0.0001, " precision loss "
            
            # we have to try to pair up these (p-1)/2 values against the (p-1)/2 heights until we find one whose trace is sufficiently close to a rational number          
            hs = self.h_rhos
            n = (p-1)//2
            assert len(hs) == n and len(vals) == n
            Sn = SymmetricGroup(n)
            co = 0
            mini = 1000
            for g in Sn:
                Qt = []
                for i in range(1,n+1):
                    Qrho = vals[i-1] * d_rhos/ hs[g(i)-1]/self.Om/self.Om_eps
                    Qt.append(Qrho)
                trQt = sum(Qt)
                trQt_approx = trQt.nearby_rational(max_denominator=bo)
                err = abs(trQt_approx-trQt)
                if err < 0.0001:
                    print "permutation found.",
                    if err < mini:
                        print "Setting the trace to %s with error %s"%(trQt_approx,err)
                        self.Qtilde_rhos = Qt
                        self.Qtilde_trace_rhos = trQt_approx
                        self.Qtilde_norm_rhos = prod(Qt).nearby_rational(max_denominator=bo)
                        mini = err
                    else:
                        print "With larger error %s. We ignore it, its trace would be %s"%(err, trQt_approx)
                    co += 1
                      
            if co > 1:
                print "Warning: found more than one value. The chosen one might be wrong!!!"
                
            self.Qtilde_rho_poly = Qt[0].algdep( (p-1)//2 )
            
            print "the equivariant algebraic L-values are Qtilde_one = %s, Qtilde_eps = %s, Qtilde_rho has minimal polynomial  %s"%(self.Qtilde_eins,self.Qtilde_eps,self.Qtilde_rho_poly)
 
 
    def truncated_L_values(self):
        r"""
        This method truncates the previously computed algebraic L-values.
        It sets the valeus Q_eins, Q_eps, Q_rhos, Q_trace_rhos and Q_norm_rhos 
        
        needs over_Q, over_K, over_L, find_Q, eq_regs and eq_L_values
        
        EXAMPLE::
            
            sage: R.<X> = QQ[]
            sage: d = etnc_example("50b1", X^3 - X^2 - 5*X + 4, prec=16)
            sage: d.truncated_L_values()
            -- start eq_L_values 
            the L-values twisted by 2-dimensional representations are  [16.5156716492158]
            -- start over_Q 
            bsd over Q asserts 0.200000000000000 = 1/5 * Sha 
            -- start find_points_over_L 
            magma finds the point  (2*t^2 - 6*t + 1 : 8*t^2 - 24*t + 14 : 1)
            saturate at 2
            saturate at 3
            saturate at 5
            saturate at 7
            this may not yet be a basis of E(L). We assume it however for the rest !
            -- start over_L 
            bsd over L of discriminant 7 * 67 asserts 1.00000000000000 = 1 * Sha 
            -- start over_K 
            found point on twist (-447371/1681 : 1403243866/68921 : 1)
            gives point (-110738/112627 : 56900929/5568954642*s0 + 175802378/2784477321 : 1) over K=QQ(sqrt(469))
            saturate at 2
            bsd over K=Q(sqrt(469)) asserts 0.199999999999999 = 1/5 * Sha 
            -- start find_Q 
            the point Q was obtained from S = (2) and is equal to (15491/898537*s^5 - 159559/3594148*s^4 - 1034283/898537*s^3 + 8924971/3594148*s^2 + 28301039/1797074*s - 45735689/3594148 : 3542155/60201979*s^5 - 37240335/481615832*s^4 - 662158945/180605937*s^3 + 7056897217/1444847496*s^2 + 30645994063/722423748*s - 58833455135/1444847496 : 1) 
            -- start eq_regs 
            the equivariant regulators are H_eins = 1, H_eps = 24.918817743633364873940975569685984917863261491634259562394 and H_rho = [12.502017857243808799928029364651848928000017580974916618997]
            -- start alg_L_values 
            permutation found. Setting the trace to 5/4 with error 4.44089209850063e-15
            the equivariant algebraic L-values are Qtilde_one = 1/5, Qtilde_eps = 1, Qtilde_rho has minimal polynomial  4*x - 5
            -- start truncated_L_values 
            the truncated equivariant L-values are Q_one = 110/469, Q_eps = 1, sum Q_rho = 1375/938

            sage: d.Q_rhos
            [1.46588486140725]

 
        """
        if hasattr(self,"Q_rhos"):
            print "Already done."
        else:            
            if not hasattr(self,"Qtilde_eps"):
                self.alg_L_values()
            print "-- start truncated_L_values "

            F = self.F
            p = self.p
            Sr = F.disc().support()
            u_eins = (-1)**len(Sr)
            u_eps = +1
            u_rho = +1
            for v in Sr:
                w = F.primes_above(v)[0]
                fw = w.residue_class_degree() 
                ew = w.ramification_index()
                if fw == 1 and ew == p:
                    u_eps *= -1
                if ew == 2:
                    assert fw == 1
                    u_rho *= -1

            E = self.E
            t_eins = 1; t_eps = 1; t_rho = 1
            for v in Sr:
                Nv = E.Np(v)
                w = F.primes_above(v)[0]
                fw = w.residue_class_degree() 
                ew = w.ramification_index()
                t_eins *= Nv/v
                if fw == 1 and ew == p:
                    t_eps *= Nv/v
                elif ew == 2:
                    t_rho *= Nv/v
                else: #ew = p fw = 2
                    w = self.K.primes_above(v)[0]
                    qw = v ** 2 
                    Nw = generalised_Nv(self.EK,w)
                    t_eps *= Nw/qw * v/Nv

            self.Q_eins  = u_eins * t_eins * self.Qtilde_eins
            self.Q_eps   = u_eps  * t_eps  * self.Qtilde_eps
            self.Q_rhos  = [u_rho * t_rho  * qq for qq in self.Qtilde_rhos]
            self.Q_trace_rhos = u_rho * t_rho * self.Qtilde_trace_rhos
            self.Q_norm_rhos =  u_rho * t_rho * self.Qtilde_norm_rhos
            print "the truncated equivariant L-values are Q_one = %s, Q_eps = %s, sum Q_rho = %s"%(self.Q_eins,self.Q_eps,self.Q_trace_rhos)

    def check_congruence(self): 
        r"""
        This method finally checks the congruences predicted by etnc_p. It needs the prior
        evaluation of all the other methods.
       
        EXAMPLE::
        
            sage: R.<X> = QQ[]
            sage: d = etnc_example("43a1", X^3 - X^2 - 5*X + 4, prec=16)
            sage: d.check_congruence()
            -- start eq_L_values 
            the L-values twisted by 2-dimensional representations are  [32.0764704468899]
            -- start over_Q 
            bsd over Q asserts 1.00000000000000 = 1 * Sha 
            -- start find_points_over_L 
            magma finds the point  (17/4*t^2 + 5*t - 6 : 181/8*t^2 + 267/8*t - 81/2 : 1)
            saturate at 2
            saturate at 3
            saturate at 5
            saturate at 7
            this may not yet be a basis of E(L). We assume it however for the rest !
            -- start over_L 
            bsd over L of discriminant 7 * 67 asserts 0.999999999999990 = 1 * Sha 
            -- start over_K 
            saturate at 2
            bsd over K=Q(sqrt(469)) asserts 1.00000000000000 = 1 * Sha 
            -- start find_Q 
            norm of P is 1 * R
            the point Q is 1 * P + 0 * R = (29/13411*s^5 + 1341/26822*s^4 - 2828/40233*s^3 - 420827/160932*s^2 + 81445/80466*s + 1258951/40233 : 5251/80466*s^5 + 60637/321864*s^4 - 104117/26822*s^3 - 3016765/321864*s^2 + 8598967/160932*s + 10656817/80466 : 1)
            -- start eq_regs 
            the equivariant regulators are H_eins = 0.12563301417497529853141758293393737263798407454502063305996, H_eps = 1 and H_rho = [23.227713433565741183086444372106462739426731778729526964157]
            -- start alg_L_values 
            permutation found. Setting the trace to 1 with error 1.05471187339390e-14
            the equivariant algebraic L-values are Qtilde_one = 1/2, Qtilde_eps = 2, Qtilde_rho has minimal polynomial  x - 1
            -- start truncated_L_values 
            the truncated equivariant L-values are Q_one = 284/469, Q_eps = 2, sum Q_rho = 568/469
            etnc_3 holds numerically for this example as 568/469 == -1136/469 (mod 3)

        """
        if not hasattr(self,"Q_eps"):
            self.truncated_L_values()
            
        lhs = self.Q_eins * self.Q_eps
        rhs = -2* self.Q_trace_rhos
        p = self.p
        if valuation(lhs,p) == 0 and valuation(rhs,p) == 0 and lhs % p == rhs % p:
            print "etnc_%s holds numerically for this example as %s == %s (mod %s)"%(p,lhs,rhs,p)
        else:
            print "etnc_%s DOES NOT hold numerically for this example. lhs = %s, rhs =%s"%(p,lhs,rhs)
            
# ========================
# function for generating interesting fields
            
def dihedral_fields(p,at_infinity=None):
    r"""
    Returns a list of polynomials that give dihedral extensions of small disciminant.
    
    These polynomials were copied from the tables by Klueners and Malle:
    http://www.math.uni-duesseldorf.de/~klueners/minimum/
    
    INPUT ::
    
        - ``p`` - an odd small prime, currently p<= 11.
        - ``at_infinity`` (optional) - either "real", "imaginary" or "all" (default)
        
    EXAMPLES::
    
        sage: dihedral_fields(5,"real")
        [X^5 - X^4 - 5*X^3 + 4*X^2 + 3*X - 1, X^5 - X^4 - 6*X^3 + 5*X^2 + 3*X - 1]
        
    """
    if not p.is_prime() or p == 2:
        raise ValueError
    
    R = PolynomialRing(QQ,"X")
    X = R.gens()[0] 

    if p == 3 :
        re = [X**3-X**2-3*X+1,X**3-4*X-1,X**3-X**2-4*X+3,X**3-X**2-4*X+2,X**3-X**2-4*X+1,X**3-X**2-5*X-1,X**3-X**2-5*X+4,X**3-5*X-1,X**3-X**2-5*X+3,X**3-X**2-6*X-2,X**3-6*X-3,X**3-7*X-5,X**3-X**2-7*X+8,X**3-6*X-2,X**3-X**2-6*X-1,X**3-X**2-6*X+5,X**3-X**2-7*X-3,X**3-6*X-1,X**3-X**2-8*X+10,X**3-7*X-4,X**3-X**2-6*X+1,X**3-X**2-6*X+3]
        im = [X**3-X**2+1,X**3-X**2-1,X**3-X**2+X+1,X**3-X**2-X+2,X**3-X**2+3*X-1,X**3-X**2+X-2,X**3-X**2+2*X+1,X**3-X-2,X**3-X**2+3*X-2,X**3-2,X**3-X**2-2,X**3+3*X-1,X**3-X**2+X+2,X**3+2*X-2,X**3-X**2-2*X-2,X**3-X**2-X+3,X**3-X**2+2*X-3, X**3-X**2+4*X-1,X**3-X**2+2*X+2,X**3-X**2+X-3,X**3-2*X-3,X**3-X**2+4*X-2,X**3+3*X-2,X**3-X**2+3,X**3-X-3, X**3-3,X**3-X**2-4*X+6,X**3+X-3, X**3-X**2-3,X**3-X**2-3*X+5,X**3+4*X-1,X**3-X**2-3*X-3, X**3-X**2+3*X+2, X**3-3*X-4, X**3-X**2-2*X-3,X**3-X**2+3*X-4,X**3-X**2+4*X+1,X**3-X**2-X+4,X**3+3*X-3,X**3-X**2+X+7,X**3+4*X-2,X**3-X**2+2*X+3,X**3-X**2+X-4]
    elif p == 5 :
        re = [X**5-X**4-5*X**3+4*X**2+3*X-1,X**5-X**4-6*X**3+5*X**2+3*X-1]
        im = [X**5-2*X**4+2*X**3-X**2+1,X**5-2*X**4+3*X**2-2*X+1,X**5-2*X**4+3*X**3-3*X**2+X+1,X**5-X**4-X**2+3*X-1,X**5-X**4-2*X**3+X**2+3*X-1,X**5-X**4+2*X**3-X**2+X+2, X**5-X**4-X**3+3*X-1, X**5+X**3-3*X**2+X-3,X**5-X**4+3*X**2-X+2,X**5-X**3-2*X**2+3*X+4,X**5-2*X**4+4*X**3-5*X**2+2*X+1,X**5-X**4+3*X**3-3*X**2+2*X+2,X**5+3*X**3-X**2+3*X+3,X**5-X**4-X**3-3*X**2+6*X+1,X**5-X**4-X**3+5*X**2-2*X+2,X**5-X**4+4*X**3-X**2+5*X-4,X**5-5*X**2-3,X**5-X**4-9*X**2+13*X-1,X**5-3*X**3-4*X**2+8*X+1,X**5-X**4-X**3+5*X**2+3*X-5,X**5-2*X**4+3*X**3-9*X**2+5*X-5,X**5-2*X**4+5*X**3-7*X**2+3*X+1,X**5-2*X**4+3*X**3+10*X+8,X**5-5*X**3+10*X-4,X**5-X**4+3*X**3+5*X**2-X+3,X**5-X**4-2*X**3+7*X**2-3*X+2,X**5-X**4+3*X**3-9*X**2+6*X-5,X**5+3*X**3-6*X**2-X-4,X**5-X**4+7*X**3+2*X**2-5*X+5,X**5-X**4+4*X**3-5*X**2+3*X+2,X**5-X**4+3*X**3-7*X**2-X-5,X**5-X**4-4*X**2-X-27,X**5-2*X**4-X**3+6*X**2+2*X-13,X**5-2*X**4+5*X**3+6*X**2-7*X+2,X**5-2*X**4+3*X**3-8*X**2-X-6, X**5-2*X**4-X**3+8*X**2+3*X-6]
    elif p == 7 :
        re = [X**7-2*X**6-7*X**5+10*X**4+13*X**3-10*X**2-X+1, X**7-X**6-9*X**5+2*X**4+21*X**3+X**2-13*X-1, X**7-3*X**6-13*X**5+13*X**4+31*X**3-17*X**2-18*X+9]
        im = [X**7-2*X**6+2*X**5+X**3-3*X**2+X-1, X**7-X**6+X**5+3*X**3-X**2+3*X+1 , X**7-3*X**6+3*X**5-5*X**4+6*X**3-X**2+5*X-1, X**7+X**5-4*X**4-X**3+5*X+1,X**7+4*X**5-2*X**4+2*X**3-3*X+2,X**7-2*X**5-X**4-2*X**3+6*X+5,X**7+7*X**3-7*X**2+7*X+1,X**7-2*X**6+X**5-X**4-5*X**2-6*X-4,X**7-2*X**6+2*X**5-X**3-6*X**2-6*X-4,X**7-2*X**6-2*X**4+X**3-4*X**2-6*X-4,X**7-2*X**6+3*X**5+X**4-2*X**3-7*X**2-6*X-4,X**7-2*X**6-X**5-3*X**4+2*X**3-3*X**2-6*X-4,X**7-2*X**6-2*X**5-4*X**4+3*X**3-2*X**2-6*X-4,X**7-2*X**6+4*X**5+2*X**4-3*X**3-8*X**2-6*X-4,X**7-2*X**6-3*X**5-5*X**4+4*X**3-X**2-6*X-4,X**7-2*X**6-4*X**5-6*X**4+5*X**3-6*X-4,X**7+7*X**5-21*X**4+70*X**3+273*X**2+210*X+900,X**7+14*X**5-70*X**4+133*X**3-238*X**2+112*X-96,X**7-21*X**5+14*X**4+294*X**3-812*X**2+784*X-160]
    elif p == 11 :
        re = [ X**11-5*X**10-4*X**9+54*X**8-53*X**7-127*X**6+208*X**5+69*X**4-222*X**3+29*X**2+56*X-5,X**11+5*X**10-18*X**9-114*X**8+37*X**7+856*X**6+746*X**5-2117*X**4-3648*X**3-135*X**2+2677*X+1301,X**11+5*X**10-36*X**9-72*X**8+333*X**7+132*X**6-604*X**5-314*X**4+306*X**3+303*X**2+90*X+9,X**11-8*X**10-347*X**9+5306*X**8-21275*X**7-19744*X**6+259087*X**5-298510*X**4-363642*X**3+787828*X**2-440424*X+79056, X**11-859*X**9+7731*X**8+180390*X**7-2931767*X**6-1205177*X**5+270949216*X**4-1779001885*X**3+1542227125*X**2+20251397450*X-52943155025]
        im = [X**11-X**10+5*X**8+8*X**5+6*X**4-X**3+X**2+3*X+1, X**11-79*X**10-57*X**9-33*X**8-68*X**7-73*X**6-59*X**5-67*X**4-57*X**3-30*X**2-8*X-1]
    else :
        print "None listed yet for this p=%s"%p
        return []
    
    if at_infinity == None or at_infinity == "all":
        return re + im
    elif at_infinity == "real":
        return re
    elif at_infinity == "imaginary":
        return im
    else:
        raise ValueError
        
# real for p^2 = 9 is  [X**9+X**8-16*X**7+6*X**6+64*X**5-65*X**4-48*X**3+78*X**2-19*X-1, X**9-95*X**7+430*X**6-285*X**5-1026*X**4+1416*X**3-152*X**2-304*X+64,X**9+4*X**8-111*X**7-37*X**6+3389*X**5-8421*X**4-796*X**3+10067*X**2+5176*X+644, X**9-X**8-112*X**7-43*X**6+2580*X**5-73*X**4-8921*X**3-9620*X**2-3764*X-508,X**9-3*X**8-174*X**7+254*X**6+7527*X**5-7137*X**4-54306*X**3+59022*X**2+102816*X-124848,X**9+3*X**8-229*X**7-670*X**6+10150*X**5+51226*X**4+72093*X**3+39191*X**2+8955*X+720,X**9+3*X**8-265*X**7+1109*X**6+9297*X**5-61461*X**4+15125*X**3+400127*X**2-467422*X-193218]

# and im is  [ X**9-3*X**8+4*X**7-5*X**6+6*X**5-X**4-5*X**3+4*X**2-2, X**9-X**8-3*X**6+3*X**3+3*X**2+5*X+1,X**9+2*X**8+3*X**7+5*X**6-2*X**5+X**4+2*X**3-2*X**2+20*X-7,X**9-X**8+7*X**7-10*X**6-20*X**5+18*X**4+36*X**3-19*X**2-10*X-1,X**9+4*X**8+4*X**7+X**6+5*X**5-40*X**4-17*X**3+82*X**2-38*X+23,X**9+4*X**8+12*X**7+26*X**6+32*X**5+59*X**4-3*X**3+176*X**2-89*X+125,X**9+9*X**7-6*X**6+27*X**5-36*X**4+27*X**3-54*X**2-32,X**9+6*X**8+5*X**7+36*X**6-37*X**5+54*X**4-16*X**3+33*X**2-2*X+9,X**9-9*X**6+27*X**3-3,X**9+9*X**7-51*X**6+72*X**5+27*X**4-150*X**3+180*X**2-144*X+64,X**9-2*X**8+25*X**7+9*X**6+114*X**5+247*X**4+158*X**3+572*X**2+302*X+49,X**9-2*X**8-8*X**7+41*X**6+44*X**5-156*X**4-13*X**3+382*X**2+236*X+59]