Ticket #10973: Trac10973.7.patch

File Trac10973.7.patch, 50.7 KB (added by John Cremona, 11 years ago)

Replaces all previous: reviewer's patch

  • sage/schemes/elliptic_curves/all.py

    # HG changeset patch
    # User Author Justin C. Walker <justin@mac.com>
    # Date 1301015058 25200
    # Node ID 3402c51950f4f73ea4dd57ef5abb0d6fba4570d4
    # Parent  9e29a3d84c48c399daaf3920bcb8b17273a0e876
    Trac 10973: Corrected report2.txt issues
    
    diff --git a/sage/schemes/elliptic_curves/all.py b/sage/schemes/elliptic_curves/all.py
    a b  
    4040from ell_curve_isogeny import EllipticCurveIsogeny, isogeny_codomain_from_kernel
    4141
    4242from heegner import heegner_points, heegner_point
     43
     44
  • sage/schemes/elliptic_curves/ell_egros.py

    diff --git a/sage/schemes/elliptic_curves/ell_egros.py b/sage/schemes/elliptic_curves/ell_egros.py
    a b  
    418418        # will preserve S-integrality of points.
    419419        E2 = E.local_minimal_model(2) if 2 in S else E
    420420        E23 = E2.local_minimal_model(3) if 3 in S else E2
    421         urst = E23.isomorphism_to(E)    
     421        urst = E23.isomorphism_to(E)
    422422       
    423423        try:
    424424            pts = E23.S_integral_points(S,proof=proof)
  • new file sage/schemes/elliptic_curves/ell_int_points.py

    diff --git a/sage/schemes/elliptic_curves/ell_int_points.py b/sage/schemes/elliptic_curves/ell_int_points.py
    new file mode 100644
    - +  
     1r"""
     2Integral Points of Elliptic Curves over Number Fields
     3
     4The following examples are from trac ticket #10152.  Previously Magma was finding more
     5integral points on these curves than Sage.  We resolved this issue.
     6
     7EXAMPLES::
     8
     9        sage: E = EllipticCurve('2082a1')
     10        sage: E.integral_points(algorithm = "old")
     11        [(-11 : 29 : 1), (-2 : 29 : 1), (4 : 11 : 1), (13 : 29 : 1)]
     12        sage: E.integral_points()
     13        [(-11 : -19 : 1), (-2 : 29 : 1), (4 : -16 : 1), (13 : -43 : 1), (507525709 : 11433453531221 : 1)]
     14   
     15    ::
     16   
     17        sage: E = EllipticCurve('6104b1')
     18        sage: E.integral_points(algorithm = "old")
     19        [(-41 : 266 : 1), (-27 : 294 : 1), (-13 : 266 : 1), (19 : 74 : 1), (22 : 49 : 1), (27 : 6 : 1),
     20        (29 : 14 : 1), (43 : 154 : 1), (53 : 266 : 1), (71 : 490 : 1), (414 : 8379 : 1), (1639 : 66346
     21        : 1), (1667 : 68054 : 1)]
     22        sage: E.integral_points()
     23        [(-41 : -266 : 1), (-27 : 294 : 1), (-13 : 266 : 1), (19 : 74 : 1), (22
     24        : -49 : 1), (27 : -6 : 1), (29 : -14 : 1), (43 : -154 : 1), (53 : 266 :
     25        1), (71 : 490 : 1), (414 : -8379 : 1), (1639 : 66346 : 1), (1667 :
     26        -68054 : 1), (23036693 : 110568192854 : 1)]
     27   
     28    ::
     29   
     30        sage: E = EllipticCurve('8470g1')
     31        sage: E.integral_points(algorithm = "old")
     32        [(-1 : 0 : 1), (0 : 6 : 1), (3 : 12 : 1), (10 : 33 : 1), (15 : 56 : 1), (24 : 110 : 1), (43 : 264
     33        : 1), (98 : 924 : 1), (395 : 7656 : 1)]
     34        sage: E.integral_points()
     35        [(-1 : 0 : 1), (0 : -7 : 1), (3 : 12 : 1), (10 : -44 : 1), (15 : 56 :
     36        1), (24 : 110 : 1), (43 : 264 : 1), (98 : -1023 : 1), (395 : -8052 : 1),
     37        (1681690 : 2179974753 : 1)]
     38   
     39    ::
     40   
     41        sage: E = EllipticCurve('8470g2')
     42        sage: E.integral_points(algorithm = "old")
     43        [(-14 : 17 : 1), (-12 : 33 : 1), (-7 : 38 : 1), (-1 : 22 : 1), (0 : 17 : 1), (13 : 8 : 1), (14
     44        : 17 : 1), (21 : 66 : 1), (33 : 158 : 1), (44 : 257 : 1), (63 : 458 : 1), (98 : 913 : 1), (175
     45        : 2222 : 1), (1533 : 59258 : 1), (1561 : 60896 : 1)]
     46        sage: E.integral_points()
     47        [(-14 : 17 : 1), (-12 : -22 : 1), (-7 : 38 : 1), (-1 : 22 : 1), (0 : -18
     48        : 1), (13 : -22 : 1), (14 : -32 : 1), (21 : 66 : 1), (33 : -192 : 1),
     49        (44 : 257 : 1), (63 : 458 : 1), (98 : -1012 : 1), (175 : -2398 : 1),
     50        (1533 : -60792 : 1), (1561 : 60896 : 1), (18888968 : -82103620687 : 1)]
     51       
     52    ::
     53   
     54        sage: from sage.schemes.elliptic_curves.ell_int_points import integral_points
     55        sage: len(EllipticCurve('15a1').integral_points(algorithm = "old")) == len(EllipticCurve('15a1').integral_points())
     56        True
     57
     58
     59
     60
     61AUTHORS:
     62
     63- this file is based on work of Thotsaphon "Nook" Thongjunthug that was initially
     64  ported to Sage by Radoslav Kirov and Jackie Anderson.
     65
     66- Radoslav Kirov, Jackie Anderson (2010): first version
     67
     68- Justin Walker, Gagan Sekhon (2011): clean up
     69
     70- Alyson Deines, Jen Balakrishnan (2011): Documentation, clean up, fixing bugs so this works for all examples in the paper [SMART]_
     71
     72REFERENCE:
     73
     74.. [SMART] Smart, N. P. and Stephens, N. M., Integral Points on Elliptic Curves over Number Fields.
     75Math. Proc. Camb. Phil. Soc. (1997), 122, 9
     76
     77"""
     78
     79
     80##############################################################################
     81#       Copyright (C) 2011 William Stein <wstein@gmail.com>
     82#
     83#  Distributed under the terms of the GNU General Public License (GPL)
     84#
     85#    This code is distributed in the hope that it will be useful,
     86#    but WITHOUT ANY WARRANTY; without even the implied warranty of
     87#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
     88#    General Public License for more details.
     89#
     90#  The full text of the GPL is available at:
     91#
     92#                  http://www.gnu.org/licenses/
     93##############################################################################
     94
     95
     96from copy import copy
     97from sage.functions.all import sqrt
     98from sage.matrix.all import zero_matrix, MatrixSpace, matrix
     99from sage.modules.free_module_element import vector
     100from sage.misc.all import verbose, n, prod
     101from sage.rings.all import is_RealField, polygen, ZZ, RealField, ComplexField
     102from sage.rings.all import QQ, integer_ceil, integer_floor, PolynomialRing
     103from sage.schemes.elliptic_curves.ell_number_field import EllipticCurve
     104
     105
     106def abs_log_height(X, gcd_one = True, precision = None):
     107    r"""
     108    Computes the height of a point in a projective space over field `K`.
     109    It assumes the coordinates have gcd equal to 1.
     110    If not use gcd_one flag.
     111   
     112    INPUT:
     113   
     114    - ``X`` -- Point in projective space over a number field `K`
     115   
     116    - ``gcd_one`` -- (default: True) if false this throws in the
     117      places at the numerators in addition to the places at the denominator.
     118   
     119    - ``precision`` -- (default: None) bits of precision of the real and complex fields
     120   
     121    OUTPUT:
     122   
     123    The absolute logarithmic height of the point.
     124   
     125   
     126    EXAMPLES::
     127   
     128        sage: import sage.schemes.elliptic_curves.ell_int_points as ellpts
     129        sage: E = EllipticCurve('5077a')
     130        sage: P1,P2,P3 = E.gens()
     131        sage: ellpts.abs_log_height(list(P2))
     132        1.09861228866811
     133       
     134    ::
     135   
     136        sage: import sage.schemes.elliptic_curves.ell_int_points as ellpts
     137        sage: K.<a> = NumberField(x^2-x-1)
     138        sage: F = EllipticCurve(K,[0,-a,1,-a-1,2*a+1])
     139        sage: P1,P2 = F.gens()
     140        sage: P2
     141        (-3/4*a + 1/4 : -5/4*a - 5/8 : 1)
     142        sage: ellpts.abs_log_height(list(P2))
     143        2.56625746464637
     144
     145   
     146    """
     147    assert isinstance(X,list)
     148    K = X[0].parent()
     149    if precision is None:
     150        RR = RealField()
     151        CC = ComplexField()
     152    else:
     153        RR = RealField(precision)
     154        CC = ComplexField(precision)
     155    places = set([])
     156    if K == QQ:
     157        embs = K.embeddings(RR)
     158        Xideal = X
     159    else:
     160        embs = K.places(precision)
     161        # skipping zero as it currently over K breaks Sage
     162        Xideal = [K.ideal(x) for x in X if not x.is_zero()]
     163    for x in Xideal:
     164        places = places.union(x.denominator().prime_factors())
     165        if not gcd_one:
     166            places = places.union(x.numerator().prime_factors())
     167    if K == QQ:
     168        non_arch = sum([max([RR(P)**(-x.valuation(P)) for x in X]).log() for P in places])
     169    else:
     170        non_arch = sum([P.residue_class_degree() *
     171                        P.absolute_ramification_index() *
     172                        max([x.abs_non_arch(P, precision) for x in X]).log() for P in places])
     173    arch = 0
     174    r,s = K.signature()
     175    for i,f in enumerate(embs):
     176        if i<r:
     177            arch += max([f(x).abs() for x in X]).log()
     178        else:
     179            arch += max([f(x).abs()**2 for x in X]).log()
     180    return (arch+non_arch)/K.degree()
     181   
     182def _compute_c6(E,emb):
     183    r"""
     184    Computes the `c_6` invariant from [SMART]_
     185   
     186    EXAMPLES::
     187   
     188        sage: import sage.schemes.elliptic_curves.ell_int_points as ellpts
     189        sage: K.<a> = NumberField(x^2+5)
     190        sage: E = EllipticCurve(K,[-112,400])
     191        sage: embs = K.embeddings(CC)
     192        sage: ellpts._compute_c6(E,embs[0])
     193        24.0994429807464
     194
     195    """
     196    F = emb.codomain()
     197    if is_RealField(F):
     198        x = polygen(ComplexField(F.prec()))
     199    else:
     200        x = polygen(F)
     201    #f = x**3-27*emb(E.c4())*x-54*emb(E.c6())
     202    f = x**3-emb(E.c4()/48)*x-emb(E.c6()/864)
     203    R = f.roots(multiplicities = False)
     204    m = max([r.abs() for r in R])
     205    return 2*m
     206
     207def _compute_c8(L):
     208    r"""
     209    Computes the `c_8` invariant from [SMART]_
     210   
     211    Original code transformed the lattice generators.
     212    Here we assume they are of standard form.
     213   
     214    INPUT:
     215   
     216    - `L` -- a basis for the period lattice of an elliptic curve over a number field
     217      via a given embedding
     218   
     219    EXAMPLES::
     220   
     221        sage: import sage.schemes.elliptic_curves.ell_int_points as ellpts
     222        sage: K.<a> = NumberField(x^2+5)
     223        sage: E = EllipticCurve(K,[-112,400])
     224        sage: embs = K.embeddings(CC)
     225        sage: Periods = E.period_lattice(embs[0]).normalised_basis();
     226        sage: ellpts._compute_c8(Periods)
     227        1.27486778253862
     228   
     229    """
     230    CC=ComplexField()
     231    w1, w2 = L
     232    m = max(CC(w1).abs(), CC(w2).abs(), CC(w1+w2).abs())
     233    return m
     234
     235def _c3(E,L):
     236    r"""
     237    Computes the `c_3` invariant of [SMART]_ which is the least eigenvalue of the
     238    height pairing matrix of the Mordell-Weil basis.
     239   
     240    INPUT:
     241   
     242    - `E` -- an elliptic curve over a number field
     243   
     244    - `L` -- a list of points of `E` which generate the Mordell-Weil group
     245   
     246    EXAMPLES::
     247   
     248        sage: import sage.schemes.elliptic_curves.ell_int_points as ellpts
     249        sage: K.<a> = NumberField(x^2+5)
     250        sage: E = EllipticCurve(K,[-112,400])
     251        sage: ellpts._c3(E,[E((4,-4)),E((0,-20)),E((-4,-28)),E((-89/5,-637*a/25))])
     252        0.366347735724767
     253
     254   
     255    """
     256    import warnings
     257    with warnings.catch_warnings():
     258       warnings.simplefilter("ignore")
     259       return min(map(abs,E.height_pairing_matrix(L).eigenvalues()))
     260
     261def _h_E(E):
     262    r"""
     263    Computes the height of the ellpitc curve `E`.
     264       
     265    INPUT:
     266   
     267    -`E` -- an elliptic curve
     268   
     269    EXAMPLES::
     270   
     271        sage: import sage.schemes.elliptic_curves.ell_int_points as ellpts
     272        sage: K.<a> = NumberField(x^2+5)
     273        sage: E = EllipticCurve(K,[-112,400])
     274        sage: ellpts._h_E(E)
     275        17.4513334798896
     276   
     277    """
     278    K = E.base_field()
     279    j = E.j_invariant()
     280    c4, c6 = E.c_invariants()
     281    g2, g3 = c4/12, c6/216
     282    return max(abs_log_height([K(1), g2, g3]), abs_log_height([K(1), j]))
     283
     284def _h_m(E, P, ElogEmbedP, D7):
     285    r"""
     286    Computes the list of modified height of a point `P` on `E(K)` given a list of embeddings
     287   
     288    INPUT:
     289   
     290    - `E` -- an elliptic curve
     291   
     292    - `P` -- a point on `E`
     293   
     294    - `ElogEmbedP` -- elliptic logarithm of `P` in some embedding
     295   
     296    - `D7` -- constant `d_7`
     297   
     298    EXAMPLES::
     299   
     300        sage: import sage.schemes.elliptic_curves.ell_int_points as ellpts
     301        sage: K.<a> = NumberField(x^2+5)
     302        sage: embs = K.embeddings(CC)
     303        sage: E = EllipticCurve(K,[-112,400])
     304        sage: L = [E((4,-4)),E((0,-20)),E((-4,-28)),E((-89/5,-637*a/25))]
     305        sage: Elog = [ P.elliptic_logarithm(embs[0], precision = 300) for P in L]
     306        sage: ellpts._h_m(E, L[0], Elog[0], 3.28)
     307        17.4513334798896
     308   
     309 
     310    """
     311    K = E.base_field()
     312    M1 =  max([P.height(), _h_E(E), D7/K.degree()*abs((ElogEmbedP).log())**2])
     313    M2 =  max([P.height(), _h_E(E), D7/K.degree()*abs(ElogEmbedP)**2])
     314    return M2
     315   
     316def _Extra_h_m(E, Periods, D7):
     317    r"""
     318    Computes two extra `h_m`'s given periods.
     319   
     320    INPUT:
     321   
     322    - `E` -- elliptic curve
     323   
     324    - `Periods` -- list of two periods of the fundamental parallelogram
     325   
     326    - `D7` -- constant `d_7`
     327   
     328    EXAMPLES::
     329   
     330        sage: import sage.schemes.elliptic_curves.ell_int_points as ellpts
     331        sage: K.<a> = NumberField(x^2+5)
     332        sage: embs = K.embeddings(CC)
     333        sage: E = EllipticCurve(K,[-112,400])
     334        sage: Periods = E.period_lattice(embs[0]).normalised_basis()
     335        sage: ellpts._Extra_h_m(E, Periods, 90.3)
     336        [48.6392854292010, 24.7424615832146]
     337       
     338    """
     339
     340    D = E.base_field().degree()
     341    h = _h_E(E)
     342    return map(lambda x: max([0, h, D7/D*abs(x)**2]), Periods)
     343
     344def _d8(E, L, Elog, Periods, D7):
     345    r"""
     346    The `c_8` constant (which we have renamed to `d_8` to avoid confusion with the other set of
     347    `c` constants) arrising from David's lower bound for linear forms in elliptic logarithms. 
     348    Notation as in Appendix A of Smart's book, "The Algorithmic Resolution of Diophantine Equations".
     349   
     350    INPUT:
     351   
     352    - `E` -- elliptic curve
     353   
     354    - `L` -- a list of points in `E(K)` (e.g. a basis for the Mordell-Weil group)
     355   
     356    - `Elog` -- a sequence of precomputed elliptic logs of each point in `L`
     357   
     358    - `D7` -- constant `d_7`
     359   
     360    EXAMPLES::
     361     
     362        sage: import sage.schemes.elliptic_curves.ell_int_points as ellpts
     363        sage: K.<a> = NumberField(x^2+5)
     364        sage: embs = K.embeddings(CC)
     365        sage: E = EllipticCurve(K,[-112,400])
     366        sage: Periods = E.period_lattice(embs[0]).normalised_basis()
     367        sage: L = [E((4,-4)),E((0,-20)),E((-4,-28)),E((-89/5,-637*a/25))]
     368        sage: Elog = [ P.elliptic_logarithm(embs[0], precision = 300) for P in L]
     369        sage: ellpts._d8(E, L, Elog, Periods, 90.3)
     370        47.4376426807629
     371
     372    """
     373    RR=RealField()
     374    C = [RR(1).exp() * _h_E(E)]
     375    D = E.base_field().degree()
     376    for i in range(len(L)):
     377        C.append(_h_m(E, L[i], Elog[i], D7) / D)
     378    C += [t / D for t in _Extra_h_m(E, Periods, D7)]
     379    return max(C);
     380
     381def _d9(E, L, Elog, Periods, D7):
     382    r"""
     383    The `c_9` constant (which we have renamed to `d_9` to avoid confusion with the other set of
     384    `c` constants) arrising from David's lower bound for linear forms in elliptic logarithms. 
     385    Notation as in Appendix A of Smart's book, "The Algorithmic Resolution of Diophantine Equations".
     386   
     387    INPUT:
     388   
     389    - `E` -- elliptic curve
     390   
     391    - `L` -- a list of points in `E(K)` (e.g. a basis for the Mordell-Weil group)
     392   
     393    - `Elog` -- a sequence of precomputed elliptic logs of each point in `L`
     394   
     395    - `D7` -- constant `d_7`
     396   
     397    EXAMPLES::
     398     
     399        sage: import sage.schemes.elliptic_curves.ell_int_points as ellpts
     400        sage: K.<a> = NumberField(x^2+5)
     401        sage: embs = K.embeddings(CC)
     402        sage: E = EllipticCurve(K,[-112,400])
     403        sage: Periods = E.period_lattice(embs[0]).normalised_basis()
     404        sage: L = [E((4,-4)),E((0,-20)),E((-4,-28)),E((-89/5,-637*a/25))]
     405        sage: Elog = [ P.elliptic_logarithm(embs[0], precision = 300) for P in L]
     406        sage: ellpts._d9(E, L, Elog, Periods, 90.3)
     407        2.71828182845904
     408
     409    """
     410    RR=RealField()
     411    D = E.base_field().degree()
     412    C = []
     413    for i in range(len(L)):
     414        tmp = RR(1.0).exp() * sqrt(D * _h_m(E, L[i], Elog[i], D7) / D7)
     415        C.append( tmp / abs(Elog[i]))
     416    Ehm = _Extra_h_m(E, Periods, D7)
     417    C += [RR(1).exp() * sqrt(D*Ehm[i]/D7) / abs(Periods[i]) for i in [0,1]]
     418    verbose("d9: C = %s"%(C,))
     419    return min(C);
     420
     421def _d10(E, L, Elog, Periods, D7):
     422    r"""
     423    The `c_10` constant (which we have renamed to `d_10` to avoid confusion with the other set of
     424    `c` constants) arrising from David's lower bound for linear forms in elliptic logarithms. 
     425    Notation as in Appendix A of Smart's book, "The Algorithmic Resolution of Diophantine Equations".
     426   
     427    INPUT:
     428   
     429    - `E` -- elliptic curve
     430   
     431    - `L` -- a list of points in `E(K)` (e.g. a basis for the Mordell-Weil group)
     432   
     433    - `Elog` -- a sequence of precomputed elliptic logs of each point in `L`
     434   
     435    - `D7` -- constant `d_7`
     436   
     437    EXAMPLES::
     438     
     439        sage: import sage.schemes.elliptic_curves.ell_int_points as ellpts
     440        sage: K.<a> = NumberField(x^2+5)
     441        sage: embs = K.embeddings(CC)
     442        sage: E = EllipticCurve(K,[-112,400])
     443        sage: Periods = E.period_lattice(embs[0]).normalised_basis()
     444        sage: L = [E((4,-4)),E((0,-20)),E((-4,-28)),E((-89/5,-637*a/25))]
     445        sage: Elog = [ P.elliptic_logarithm(embs[0], precision = 300) for P in L]
     446        sage: ellpts._d10(E, L, Elog, Periods, 90.3)
     447        1.19716657228519e226
     448   
     449    """
     450    RR=RealField()
     451    D = E.base_field().degree()
     452    n = len(L)+2
     453    scalar = 2 * 10**(8 + 7*n) * (2/RR(1).exp())**(2*n**2)
     454    scalar *= (n+1)**(4*n**2 + 10*n) * D**(2*n + 2)
     455    scalar *= (RR(_d9(E, L, Elog, Periods, D7)).log())**(-2*n-1)
     456    for i in range(len(L)):
     457        scalar *= _h_m(E, L[i], Elog[i], D7)
     458    scalar *= prod(_Extra_h_m(E, Periods, D7))   
     459    return scalar
     460
     461def _RHS(D, r, C9, C10, D9, D10, h, Q, expTors):
     462    r"""
     463    Computes the right hand side of the principal inequality.
     464   
     465    INPUT:
     466   
     467    - `D` -- degree of field extension
     468   
     469    - `r` -- Mordell-Weil rank
     470   
     471    - `C9` -- `c_9` invariant from [SMART]_
     472   
     473    - `C19` -- `c_{10}` invariant from [SMART]_
     474   
     475    - `D9` -- `D_9` invariant
     476   
     477    - `D10` -- `d_{10}` invariant
     478   
     479    - `h` -- `h_E(E)` where `E` is an elliptic curve
     480   
     481    - `Q` -- initial bound for the coefficients of `P_i`'s, where
     482      `P_i`'s are in the Mordell-Weil basis
     483       
     484    - `expTors` -- the exponent of the torsion subgroup of `E`
     485   
     486    EXAMPLES: 
     487   
     488        sage: import sage.schemes.elliptic_curves.ell_int_points as ellpts
     489        sage: K.<a> = NumberField(x^2+5)
     490        sage: E = EllipticCurve(K,[-112,400])
     491        sage: L = [E((4,-4)),E((0,-20)),E((-4,-28)),E((-89/5,-637*a/25))]
     492        sage: C9 =  3425.58431073226
     493        sage: C10 = 0.183173867862383
     494        sage: Q0 =  5.12358476153997
     495        sage: D8 =  47.4376426807629
     496        sage: D9 =  6.19404412950310
     497        sage: D10 =  2.39513899842427e221
     498        sage: h =  17.4513334798896
     499        sage: expTors = 1
     500        sage: ellpts._RHS(K.degree(), len(L), C9, C10, D9, D10, h, 1000000000000000000000, expTors)
     501        3.02137037238302e233
     502   
     503    """
     504    RR=RealField()
     505    bound = (RR(RR(Q*r*expTors).log()).log() + h + RR(D*D9).log())**(r+3)
     506    bound *= D10*(RR(Q*r*expTors).log() + RR(D*D9).log())
     507    bound += RR(C9*expTors).log()
     508    bound /= C10
     509    return bound
     510
     511def _InitialQ(D, r, Q0, C9, C10, D8, D9, D10, h, expTors):
     512    r"""
     513    Computes the initial bound on `Q = max_{1 \leq i \leq r}(|q_i|)`.
     514   
     515    INPUT:
     516   
     517    - `D` -- degree of field extension
     518   
     519    - `r` -- Mordell-Weil rank
     520   
     521    - `Q0` -- the `K_0` constant in [SMART]_
     522   
     523    - `C9` -- `c_9` invariant from [SMART]_
     524   
     525    - `C10` -- `c_{10}` invariant from [SMART]_
     526   
     527    - `D9` -- `D_9` invariant
     528   
     529    - `D10` -- `d_{10}` invariant
     530   
     531    - `h` -- `h_E(E)` where `E` is an elliptic curve
     532       
     533    - `expTors` -- the exponent of the torsion subgroup of `E`
     534   
     535    EXAMPLES::
     536   
     537        sage: import sage.schemes.elliptic_curves.ell_int_points as ellpts
     538        sage: K.<a> = NumberField(x^2+5)
     539        sage: E = EllipticCurve(K,[-112,400])
     540        sage: L = [E((4,-4)),E((0,-20)),E((-4,-28)),E((-89/5,-637*a/25))]
     541        sage: C9 =  3425.58431073226
     542        sage: C10 = 0.183173867862383
     543        sage: Q0 =  5.12358476153997
     544        sage: D8 =  47.4376426807629
     545        sage: D9 =  6.19404412950310
     546        sage: D10 =  2.39513899842427e221
     547        sage: h =  17.4513334798896
     548        sage: expTors = 1
     549        sage: ellpts._InitialQ(K.degree(), len(L), Q0, C9, C10, D8, D9, D10, h, expTors)
     550        118
     551
     552    """
     553
     554    RR=RealField()
     555    minQ = max(Q0, RR(D8).exp())
     556    Q = minQ + 1
     557    x = integer_ceil(Q.log(10)) # x = log_10(Q)
     558    verbose("x: %s, %s"%(type(Q), type(x)))
     559    exp_OK = 0   # the exponent that satisfies P.I.
     560    exp_fail = 0 # the exponent that fails P.I.
     561    verbose("Looping on %s < RHS(%s, %s, %s, %s %s, %s, %s, %s, %s"%( 10**(2*x), D, r, C9, C10, D9, D10, h, 10**x, expTors))
     562    while 10**(2*x) < _RHS(D, r, C9, C10, D9, D10, h, 10**x, expTors):
     563        exp_OK = x # Principal Inequality satisfied
     564        x *= 2     # double x, and retry
     565    exp_fail = x # x that fails the Principal Inequality
     566   
     567   
     568    # So now x = log_10(Q) must lie between exp_OK and exp_fail
     569    # Refine x further using "binary search"
     570    while True:
     571        x = integer_floor((exp_OK + exp_fail)/2)
     572        if 10**(2*x) >= _RHS(D, r, C9, C10, D9, D10, h, 10**x, expTors):
     573            exp_fail = x
     574        else:
     575            exp_OK = x
     576        if (exp_fail - exp_OK) <= 1:
     577            break
     578    return exp_fail # over-estimate
     579
     580
     581def _ReducedQ(E, f, LGen, Elog, C9, C10, Periods, expTors, initQ):
     582    r"""
     583    Internal reduction function which reduces the bound on `Q`.
     584   
     585    INPUT:
     586   
     587    - `E` -- elliptic curve over `K`
     588   
     589    - `f` -- embedding of `K` into `\mathbb{C}`
     590   
     591    - `LGen` -- a list of points in `E(K)` (e.g. a basis for the Mordell-Weil group)
     592   
     593    - `Elog` -- a sequence of precomputed elliptic logs of each point in `L`
     594
     595    - `C9` -- `c_9` invariant from [SMART]_
     596   
     597    - `C10` -- `c_{10}` invariant from [SMART]_
     598   
     599    - `Periods` -- list of two periods of the fundamental parallelogram
     600   
     601    - `expTors` -- the exponent of the torsion subgroup of `E`
     602     
     603    - `initQ` -- inital guess for `Q` which is then reduced by LLL
     604     
     605    EXAMPLES::
     606     
     607        sage: import sage.schemes.elliptic_curves.ell_int_points as ellpts
     608        sage: K.<a> = NumberField(x^2+5)
     609        sage: E = EllipticCurve(K,[-112,400])
     610        sage: L = [E((4,-4)),E((0,-20)),E((-4,-28)),E((-89/5,-637*a/25))]
     611        sage: embs = K.embeddings(CC)
     612        sage: Elog = [ P.elliptic_logarithm(embs[0], precision = 300) for P in L]
     613        sage: Periods = Periods = E.period_lattice(embs[0]).normalised_basis()
     614        sage: C9 =  3425.58431073226
     615        sage: C10 = 0.183173867862383
     616        sage: expTors = 1   
     617        sage: initQ = 10000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000
     618        sage: ellpts._ReducedQ(E, embs[0], L, Elog, C9, C10, Periods, expTors, initQ)
     619        13
     620   
     621    """
     622
     623    RR=RealField()
     624    w1, w2 = Periods
     625    r = len(LGen)
     626    newQ = initQ
     627    # Repeat LLL reduction until no further reduction is possible
     628    while True:
     629        Q = newQ
     630        S = r*(Q**2)*(expTors**2)
     631        T = 3*r*Q*expTors/sqrt(2.0)
     632        # Create the basis matrix
     633        C = 1
     634        while True:
     635            C *= Q**integer_ceil((r+2)/2)
     636            precbits = int(RR(C).log(2)+10)*4
     637            L = copy(zero_matrix(ZZ, r+2))
     638            # Elliptic logarithm should have precision "suitable to" C
     639            # e.g. If C = 10**100, then Re(Elog[i]) should be computed
     640            # correctly to at least 100 decimal places
     641            if precbits > Elog[0].prec():
     642                verbose( "Need higher precision, recompute elliptic logarithm ...")
     643                # Re-compute elliptic logarithm to the right precision
     644                verbose( "precision in bits %i" % precbits)
     645                Elog = [ P.elliptic_logarithm(f, precision = precbits) for P in LGen]
     646                verbose( "Elliptic logarithm recomputed")
     647                w1, w2 = E.period_lattice(f).normalised_basis( prec = precbits)
     648            # Assign all non-zero entries
     649            for i in range(r):
     650                L[i, i] = 1
     651                L[r, i]   = (C*Elog[i].real_part()).trunc()
     652                L[r+1, i] = (C*Elog[i].imag_part()).trunc()
     653            L[r , r ]   = (C*w1.real_part()).trunc()
     654            L[r , r+1 ] = (C*w2.real_part()).trunc()
     655            L[r+1, r]   = (C*w1.imag_part()).trunc()
     656            L[r+1, r+1] = (C*w2.imag_part()).trunc()
     657            # LLL reduction and constants
     658            L = L.transpose()
     659            L = L.LLL()
     660            b1 = L[0] # 1st row of reduced basis
     661            #Norm(b1) = square of Euclidean norm
     662            normb1 = sum([i**2 for i in b1])
     663            lhs = RR(2)**(-r-1)*normb1 - S
     664            if (lhs >= 0) and (sqrt(lhs) > T):
     665                break
     666        newQ = ( RR(C*C9*expTors).log() - RR(sqrt(lhs) - T).log() ) / C10
     667        newQ = integer_floor(sqrt(newQ))
     668        verbose( "After reduction, Q <= %f" % newQ)
     669        if ((Q - newQ) <= 1) or (newQ <= 1):
     670            break
     671    return newQ
     672
     673
     674def _cyc_iter(id, gens, mult, both_signs = False):
     675    r"""
     676    Cyclic group iterator
     677   
     678    INPUT:
     679   
     680    - `id` -- identity element
     681   
     682    - `gens` -- generators of the group
     683   
     684    - `mult` -- orders of the generators
     685   
     686    - ``both_signs`` --(Default: False) whether to return all elements or one per each {element, inverse} pair
     687   
     688    OUTPUT:
     689   
     690    Returns all elements of the cyclic group, remembering all intermediate steps
     691   
     692    EXAMPLE::
     693   
     694        sage: import sage.schemes.elliptic_curves.ell_int_points as ellpts
     695        sage: E = EllipticCurve('57b')
     696        sage: ord = [P.order() for P in E.torsion_points()]
     697        sage: P_list = [i for i in ellpts._cyc_iter(E(0),E.torsion_points(),ord)]
     698        sage: P_list
     699        [((0 : 1 : 0), [0, 0, 0, 0]), ((7/4 : -11/8 : 1), [0, 0, 0, 1]), ((1 :
     700        -1 : 1), [0, 0, 1, 0]), ((-3 : 1 : 1), [0, 0, 1, 1]), ((-3 : 1 : 1), [1,
     701        0, 0, 0]), ((1 : -1 : 1), [1, 0, 0, 1]), ((7/4 : -11/8 : 1), [1, 0, 1,
     702        0]), ((0 : 1 : 0), [1, 0, 1, 1])]
     703   
     704   
     705    """
     706    if len(gens) == 0:
     707        yield id,[]
     708    else:
     709        P = gens[0]
     710        cur = id
     711        if both_signs:
     712            ran = xrange(mult[0])
     713        else:
     714            ran = xrange(mult[0]/2 + 1)
     715        for k in ran:
     716            for rest, coefs in _cyc_iter(id, gens[1:], mult[1:], both_signs or k != 0):
     717                yield cur + rest, [k] + coefs
     718            cur += P
     719
     720#generates elements of form a_1G_1 + ... + a_nG_n
     721#where |a_i| <= bound and the first non-zero coefficient is positive
     722def _L_points_iter(id, gens, bound, all_zero = True):
     723    r"""
     724    Generates elements of form `a_1G_1 + ... + a_nG_n`
     725    where `|a_i| \leq` bound and the first non-zero coefficient is positive
     726   
     727    INPUT:
     728   
     729    - ``id`` -- identity element
     730   
     731    - ``gens`` -- a list of generators, `G_1,...,G_n`
     732   
     733    - ``bound`` -- a bound on the absolute value of the `a_i`'s, the coefficients of the elements in terms
     734      of the generators.
     735   
     736    - ``all_zero`` --(default: True)
     737   
     738    OUTPUT:
     739   
     740    An iterator that generates elements of the form `a_1G_1 + ... + a_nG_n`
     741   
     742    EXAMPLES::
     743   
     744        sage: import sage.schemes.elliptic_curves.ell_int_points as ellpts
     745        sage: K.<a> = NumberField(x^2+2)
     746        sage: E = EllipticCurve(K,[-16,16])
     747        sage: P_list = [[P,Pcoeff] for P,Pcoeff in ellpts._L_points_iter(E(0),[E((0,4)),E((2,-2*a))],24)]
     748        sage: P_list[0:5]
     749        [[(0 : 1 : 0), [0, 0]], [(2 : -2*a : 1), [0, 1]], [(-9/2 : -5/4*a : 1),
     750        [0, 2]], [(418/169 : 4514/2197*a : 1), [0, 3]], [(-30241/200 :
     751        5257039/4000*a : 1), [0, 4]]]
     752       
     753   
     754    """
     755    if len(gens) == 0:
     756        yield id, []
     757    else:
     758        P = gens[0]
     759        cur = id
     760        for k in xrange(bound+1):
     761            for rest, coefs in _L_points_iter(id, gens[1:], bound, all_zero = (all_zero and k == 0)):
     762                yield cur + rest, [k] + coefs 
     763                if k!=0 and not all_zero:
     764                   yield -cur + rest, [-k] + coefs
     765            cur += P
     766
     767def _integral_points_with_Q(E, L, Q):
     768    r"""
     769    Given an elliptic curve `E` over a number field, its Mordell-Weil basis L,and a positive integer Q,
     770    return the sequence of all integral points modulo `[-1]` of the form `P = q_1*L[1] + ... + q_r*L[r] + T`
     771    with some torsion point `T` and `|q_i| \leq Q`.
     772       
     773    INPUT:
     774   
     775    - ``E`` -- an elliptic curve
     776   
     777    - ``L`` -- a basis for the Mordel-Weil group of `E`
     778   
     779    - ``Q`` -- positive integer, maximum for the absolute bound on all coefficients in the linear combination of points in `L`
     780   
     781    OUTPUT:
     782   
     783    Returns the sequence of all integral points modulo `[-1]`
     784   
     785    EXAMPLE::
     786   
     787        sage: import sage.schemes.elliptic_curves.ell_int_points as ellpts
     788        sage: K.<a> = NumberField(x^2-x-1)
     789        sage: E = EllipticCurve(K,[2-a,-a-1,1-a,a,1])
     790        sage: L = E.gens()
     791        sage: ellpts._integral_points_with_Q(E,L,25)
     792        [(a : 1 : 1)]
     793       
     794    """
     795    CC=ComplexField()
     796    if L == []:
     797        try:
     798            L = E.gens()
     799        except ValueError:
     800            raise ValueError, "Generators of the Mordel-Weil group must be supplied."
     801    Tors = E.torsion_subgroup()
     802    expTors = Tors.exponent()
     803    Tpoints = E.torsion_points()
     804    id = E([0,1,0])
     805    verbose( "Generators = %s" % L)
     806    PtsList = []
     807     
     808    verbose( "Exact arithmetic")
     809    for P, Pcoeff in _L_points_iter(id, L, Q):
     810            for T in Tpoints:
     811                R = P+T
     812                if R[0].is_integral() and R[1].is_integral() and R != id:
     813                    if R not in PtsList and -R not in PtsList:
     814                        PtsList.append(R)
     815    verbose( "*"*45)
     816    return PtsList
     817   
     818
     819def _calculate_Q(E,L):
     820    r"""
     821    This calculates `Q`, one of the bounds in [SMART]_.  `Q` is a positive integer,
     822    maximum for the absolute bound on all coefficients in the linear combination of points in `L`
     823   
     824    INPUT:
     825   
     826    - ``E`` -- Elliptic Curve
     827   
     828    - ``L`` -- a basis for the Mordel-Weil group of `E`
     829   
     830    OUTPUT:
     831   
     832    A positive integer, maximum for the absolute bound on all coefficients in the linear combination of points in `L`
     833   
     834    EXAMPLE::
     835   
     836        sage: import sage.schemes.elliptic_curves.ell_int_points as ellpts
     837        sage: K.<a> = NumberField(x^2+5)
     838        sage: E = EllipticCurve(K,[-112,400])
     839        sage: L = [E((4,-4)),E((0,-20)),E((-4,-28)),E((-89/5,-637*a/25))]
     840        sage: ellpts._calculate_Q(E,L)
     841        13
     842    """
     843
     844    from math import pi as Rpi
     845    CC=ComplexField()
     846   
     847    if len(L) == 0:
     848        return 0
     849    K = E.base_ring()
     850    expTors = E.torsion_subgroup().exponent()
     851    r, s = K.signature()
     852    RR = RealField()
     853    pi = RR(Rpi)
     854    b2 = E.b_invariants()[0]
     855    # Global constants (independent to the embedding of E)
     856    gl_con = {}
     857    #C2 = E.silverman_height_bounds[1] matches with [SMART]
     858    #C2 = -E.silverman_height_bound[0] matches with Nook's implementation in Magma
     859    #both return correct sets of integral points . . .
     860    gl_con['c2'] = C2 = -E.silverman_height_bounds()[0]
     861    gl_con['c3'] = C3 = _c3(E,L)
     862    gl_con['h_E'] = h = _h_E(E)
     863    verbose( "Global constants")
     864    for name, val in gl_con.iteritems():
     865        verbose( '%s = %s'%(name, val))
     866    verbose( "-"*45)
     867    Q = []
     868    # Find the most reduced bound on each embedding of E
     869    # Sage bug, QQ.places() is not implemented
     870    # and K.embeddings gives each complex places twice
     871    if K is QQ:
     872        infplaces = K.embeddings(CC) 
     873    else:
     874        infplaces = K.places()
     875    for i, f in enumerate(infplaces):
     876        # Bug in P.complex_logarithm(QQ.embeddings(CC)[0])
     877        if K is QQ:
     878            emb = None
     879        else:
     880            emb = f
     881        if i < r:
     882            nv = 1
     883            verbose( "Real embedding #%i" % i )
     884        else:
     885            nv = 2
     886            verbose( "Complex embedding #%i" % (i-r))
     887        # Create complex embedding of E
     888        # Embedding of all points in Mordell-Weil basis
     889        # Find complex elliptic logarithm of each embedded point
     890        precbits = integer_floor(100*RR(10).log(2))
     891        Elog = [P.elliptic_logarithm(emb, precision = precbits) for P in L]
     892        Periods = E.period_lattice(emb).normalised_basis();
     893        # Local constants (depending on embedding)
     894        # C9, C10 are used for the upper bound on the linear form in logarithm
     895        loc_con = {}
     896        loc_con['c4'] = C4 = C3 * K.degree() / RR(nv*(r+s))
     897        loc_con['c5'] = C5 = C2 * K.degree() / RR(nv*(r+s))
     898        loc_con['c6'] = C6 = _compute_c6(E,f)
     899        loc_con['delta'] = delta = 1 + (nv-1)*pi
     900        loc_con['c8'] = C8 = _compute_c8(Periods)
     901        loc_con['c7'] = C7 = RR(8)*(delta**2) + (C8**2)*abs(f(b2))/RR(12)
     902        #C9 = sqrt(C7)*RR(C5/2).exp() is a tigher bound, but [SMART] only uses
     903        #C9 = C7*RR(C5/2).exp() . . .
     904        #additionally, we could never get our C9 values to match with [SMART]
     905        loc_con['c9'] = C9  = C7*RR(C5/2).exp()
     906        #C10 was often correct to 3 or so digits
     907        loc_con['c10'] = C10 = C4/2
     908        #Q0 (K0 in [SMART]) is often correct to 1 or 2 digits
     909        #These descrepancies don't seem to matter . . .
     910        loc_con['Q0'] = Q0 = sqrt( ( RR(C6+abs(f(b2))/12).log() + C5) / C4 )
     911        # Constants for David's lower bound on the linear form in logarithm
     912        # Magma and Sage output periods in different order, need to swap w1 and w2
     913        w2, w1 = Periods # N.B. periods are already in "standard" form
     914        loc_con['d7'] = D7 = 3*pi / ((abs(w2)**2) * (w2/w1).imag_part())
     915        loc_con['d8'] = D8 = _d8(E, L, Elog, Periods, D7)
     916        loc_con['d9'] = D9 = _d9(E, L, Elog, Periods, D7)
     917        loc_con['d10'] = D10 = _d10(E, L, Elog, Periods, D7)
     918        for name, val in loc_con.iteritems():
     919            verbose( "{0} = {1}".format(name, val))
     920        # Find the reduced bound for the coefficients in the linear logarithmic form
     921        loginitQ = _InitialQ(K.degree(), len(L), Q0, C9, C10, D8, D9, D10, h, expTors)
     922        verbose( "Initial Q <= 10^%f" % loginitQ)
     923        initQ = 10**loginitQ
     924        Q.append(_ReducedQ(E, emb, L, Elog, C9, C10, Periods, expTors, initQ))
     925        verbose( "-"*45)
     926    Q = max(Q)
     927    verbose( "Maximum absolute bound on coefficients = %i"% Q)
     928    return Q
     929
     930def integral_points(self, L=None, both_signs = False, algorithm = "new"):
     931    r"""
     932    Given an elliptic curve over a number field and its Mordell-Weil basis, return the sequence of all integral points modulo [-1]. 
     933    both_signs = True gives the set of all integral points.
     934   
     935    INPUT:
     936       
     937    - ``L`` -- (default: `L` = self.gens()) a basis for the Mordel-Weil group of self
     938   
     939    - ``both_signs`` -- (default: False) This gives the entire set of points, `P` and `-P` for each integral point.  If set to True,
     940      this just gives the `P` and not `-P`.
     941   
     942    - ``algorithm`` -- (default: "new") When "new", this runs the generic code in ell_int_points.py, when "old" this runs the old (and
     943      sometimes incorrect) code in ell_rational_field.py.
     944                 
     945    OUTPUT:
     946   
     947    Returns the sequence of all integral points modulo `[-1]`
     948   
     949    .. note::
     950
     951        The complexity increases exponentially in the rank of curve
     952        E. The computation time (but not the output!) depends on
     953        the Mordell-Weil basis. If mw_base is given but is not a
     954        basis for the Mordell-Weil group (modulo torsion), integral
     955        points which are not in the subgroup generated by the given
     956        points will almost certainly not be listed.
     957       
     958        All examples from [SMART]_ have been tested and return the correct integral points. 
     959
     960   
     961    EXAMPLES::
     962   
     963        sage: import sage.schemes.elliptic_curves.ell_int_points as ellpts
     964        sage: E = EllipticCurve('5077a')
     965        sage: ellpts.integral_points(E,E.gens(),both_signs = True)
     966        [(-3 : -1 : 1), (-3 : 0 : 1), (-2 : -4 : 1), (-2 : 3 : 1), (-1 : -4 : 1), (-1 : 3 : 1), (0 : -3 : 1),
     967        (0 : 2 : 1), (1 : -1 : 1), (1 : 0 : 1), (2 : -1 : 1), (2 : 0 : 1), (3 : -4 : 1), (3 : 3 : 1), (4 :
     968        -7 : 1), (4 : 6 : 1), (8 : -22 : 1), (8 : 21 : 1), (11 : -36 : 1), (11 : 35 : 1), (14 : -52 : 1),
     969        (14 : 51 : 1), (21 : -96 : 1), (21 : 95 : 1), (37 : -225 : 1), (37 : 224 : 1), (52 : -375 : 1), (52
     970        : 374 : 1), (93 : -897 : 1), (93 : 896 : 1), (342 : -6325 : 1), (342 : 6324 : 1), (406 : -8181 : 1),
     971        (406 : 8180 : 1), (816 : -23310 : 1), (816 : 23309 : 1)]
     972       
     973    Example 1 from [SMART]_::
     974   
     975        sage: K.<a> = NumberField(x^2+5)
     976        sage: E = EllipticCurve(K,[-112,400])
     977        sage: E.integral_points(L = [E((4,-4)),E((0,-20)),E((-4,-28)),E((-89/5,-637*a/25))])  # long time
     978        [(8 : 4 : 1), (148 : 1796 : 1), (1 : 17 : 1), (84 : -764 : 1), (1368 : -50596 : 1), (9 : -11 : 1), (16
     979        : 52 : 1), (12 : -28 : 1), (32 : 172 : 1), (-8 : 28 : 1), (-12 : 4 : 1), (208 : 2996 : 1), (25 : -115 :
     980        1), (44 : 284 : 1), (56 : -412 : 1), (-7 : -29 : 1), (372 : -7172 : 1), (1624 : -65444 : 1), (3264 : 186476
     981        : 1), (-4 : -28 : 1), (0 : -20 : 1), (4 : -4 : 1)]
     982
     983       
     984    Example 2 from [SMART]_::
     985   
     986        sage: K.<a> = NumberField(x^2+2)
     987        sage: E = EllipticCurve(K,[-16,16])
     988        sage: E.integral_points([E((0,4)),E((2,-2*a))])
     989        [(4*a : 8*a - 12 : 1), (-4*a : -8*a - 12 : 1), (4 : 4 : 1), (4*a - 4 : 20 : 1), (-4*a - 4 : 20 : 1), (-4 : -4
     990        : 1), (-40*a + 60 : -480*a + 316 : 1), (40*a + 60 : 480*a + 316 : 1), (8 : -20 : 1), (1 : -1 : 1), (-4*a - 10
     991        : -18*a + 28 : 1), (4*a - 10 : 18*a + 28 : 1), (24 : 116 : 1), (2 : -2*a : 1), (0 : 4 : 1)]
     992
     993    Example 3 from [SMART]_, also an example where we can not compute generators using the simon two descent::
     994     
     995        sage: K.<a> = NumberField(x^3+3*x-1)
     996        sage: E = EllipticCurve(K,[3,0])
     997        sage: E.integral_points(L = [E((1,2)),E((a,1))])
     998        [(3*a^2 + 9 : -9*a^2 - 3*a - 27 : 1), (3 : -6 : 1), (6*a^2 + 2*a + 17 : -24*a^2 - 6*a - 74 : 1), (12 : 42 : 1),
     999        (204*a^2 + 65*a + 632 : -5286*a^2 - 1704*a - 16405 : 1), (0 : 0 : 1), (a : 1 : 1), (1 : 2 : 1)]
     1000        sage: E.integral_points()
     1001        Traceback (most recent call last):
     1002        ...
     1003        ValueError: Cannot compute a provable set of generators, please supply generators of the Mordel-Weil group.
     1004       
     1005       
     1006    REFERENCE:
     1007
     1008    .. [SMART] Smart, N. P. and Stephens, N. M., Integral Points on Elliptic Curves over Number Fields.
     1009       Math. Proc. Camb. Phil. Soc. (1997), 122, 9
     1010    """
     1011   
     1012    if L is None or L == [] or L == 'auto':
     1013        try:
     1014            L = self.gens()
     1015            r = self.rank()
     1016        except ValueError:
     1017            raise ValueError, "Cannot compute a provable set of generators, please supply generators of the Mordel-Weil group."
     1018        assert len(L) == r
     1019       
     1020    if algorithm == "old":
     1021        try:
     1022            return self._integral_points_old(L, both_signs)
     1023        except:
     1024            pass
     1025     
     1026    Q = _calculate_Q(self,L)
     1027    int_points = _integral_points_with_Q(self, L, Q)
     1028    if both_signs:
     1029        int_points += [-P for P in int_points]
     1030       
     1031    #need this so doc tests will always have points in the same order
     1032    int_points.sort()
     1033    return int_points
  • sage/schemes/elliptic_curves/ell_number_field.py

    diff --git a/sage/schemes/elliptic_curves/ell_number_field.py b/sage/schemes/elliptic_curves/ell_number_field.py
    a b  
    6565
    6666- Chris Wuthrich
    6767
     68- Alyson Deines and Jen Balakrishnan 2011: silverman_height_bounds, faltings_height, integral_points
     69
    6870REFERENCE:
    6971   
    7072- [Sil] Silverman, Joseph H. The arithmetic of elliptic curves. Second edition. Graduate Texts in   
     
    9698from sage.rings.ring import Ring
    9799from sage.rings.arith import gcd, prime_divisors
    98100from sage.misc.misc import prod
     101from sage.misc.all import n
    99102import sage.databases.cremona
    100103import ell_torsion
    101104from ell_generic import is_EllipticCurve
     
    107110from sage.rings.integer import Integer
    108111from sage.rings.infinity import Infinity # just for verbose output
    109112from sage.rings.arith import valuation
     113from ell_int_points import integral_points
    110114
    111115class EllipticCurve_number_field(EllipticCurve_field):
    112116    r"""
     
    276280        prob_gens = [self(P) for P in t[2]]
    277281        self._simon_two_descent_data[lim1,lim3,limtriv,maxprob,limbigprime] = (prob_rank, two_selmer_rank, prob_gens)
    278282        return prob_rank, two_selmer_rank, prob_gens
     283       
     284    integral_points = integral_points
     285               
     286    def faltings_height(self):
     287        r"""
     288        Faltings's height of the elliptic curve `E`.
     289        This is double the formula in Silverman Math Comp 1990 page 725.
     290   
     291        EXAMPLE::
     292   
     293            sage: K.<a> = NumberField(x^2+5)
     294            sage: E = EllipticCurve(K,[-112,400])
     295            sage: E.faltings_height()
     296            4.29484994110149
     297
     298        """
     299        RR=RealField()
     300        K = self.base_field()
     301        c = RR(2).log()
     302        if self.b2().is_zero():
     303            c = 0
     304        from sage.schemes.elliptic_curves.ell_int_points import abs_log_height
     305        h1 = abs_log_height([K(self.discriminant()), K(1)])/6
     306        h2 = K(self.j_invariant()).global_height_arch()/6
     307        h3 = K(self.b2()/12).global_height_arch()
     308        return n(h1 + h2/2 + h3/2 + c)
     309               
     310    def silverman_height_bounds(self):
     311        r"""
     312        Return the Silverman height bound.  This is a positive real
     313        (floating point) number B such that for all points `P` on the
     314        curve over any number field, `|h(P) - \hat{h}(P)| \leq B`,
     315        where `h(P)` is the naive logarithmic height of `P` and
     316        `\hat{h}(P)` is the canonical height.
     317
     318        NOTES:
     319
     320           - Silverman's paper is 'The Difference Between the Weil
     321             Height and the Canonical Height on Elliptic Curves',
     322             Math. Comp., Volume 55, Number 192, pages 723-743.  We
     323             use a correction by Bremner with 0.973 replaced by 0.961,
     324             as explained in the source code to mwrank (htconst.cc).
     325             
     326        EXAMPLES::
     327       
     328            sage: E = EllipticCurve('37a1')
     329            sage: E.silverman_height_bounds()
     330            (-4.82540075818092, 4.07560050545395)
     331                       
     332        """
     333        from sage.rings.all import CC, QQ
     334        Delta   = self.discriminant()
     335        j       = self.j_invariant()
     336        b2      = self.b2()
     337        K       = self.base_field()
     338        r,s     = K.signature()
     339        infplaces = K.embeddings(CC) if K is QQ else K.places()
     340        nvs = [1 if i < r else 2 for i,_ in enumerate(infplaces)]
     341       
     342        twostar = 2 if b2 else 1
     343        from math import log
     344        from ell_int_points import abs_log_height
     345        def h(x):
     346            return abs_log_height([K(x),K(1)])
     347        def h_oo(x):
     348            return sum(nvs[i]*log(max(abs(v(x)),1)) for i,v in enumerate(infplaces))/K.degree()
     349        mu    = h(Delta)/12 + h_oo(j)/12 + h_oo(b2/12)/2 + log(twostar)/2
     350        lower = 2*(-h(j)/24 - mu - 0.961)
     351        upper = 2*(mu + 1.07)
     352        return lower, upper
    279353               
    280354    def height_pairing_matrix(self, points=None, precision=None):
    281355        r"""
  • sage/schemes/elliptic_curves/ell_rational_field.py

    diff --git a/sage/schemes/elliptic_curves/ell_rational_field.py b/sage/schemes/elliptic_curves/ell_rational_field.py
    a b  
    2727
    2828- Christian Wuthrich (2010-01): moved Galois reps and modular
    2929  parametrization in a separate file
     30 
    3031
    3132"""
    3233
     
    49964997    def height(self, precision=None):
    49974998        """
    49984999        Returns the real height of this elliptic curve. This is used in
    4999         integral_points()
     5000        _integral_points_old()
    50005001       
    50015002        INPUT:
    50025003       
     
    52165217        return ans
    52175218
    52185219    prove_BSD = BSD.prove_BSD
    5219 
    5220     def integral_points(self, mw_base='auto', both_signs=False, verbose=False):
     5220   
     5221    def _integral_points_old(self, mw_base='auto', both_signs=False, verbose=False):
    52215222        """
    52225223        Computes all integral points (up to sign) on this elliptic curve.
    52235224       
     5225        This has been replaced by a newer and more correct, but slower and more general implementation
     5226        in ell_int_points.py.
     5227       
    52245228        INPUT:
    52255229       
    52265230       
     
    52535257       
    52545258            sage: E=EllipticCurve([0,0,1,-7,6])
    52555259            sage: P1=E.point((2,0)); P2=E.point((-1,3)); P3=E.point((4,6))
    5256             sage: a=E.integral_points([P1,P2,P3]); a 
     5260            sage: a=E._integral_points_old([P1,P2,P3]); a 
    52575261            [(-3 : 0 : 1), (-2 : 3 : 1), (-1 : 3 : 1), (0 : 2 : 1), (1 : 0 : 1), (2 : 0 : 1), (3 : 3 : 1), (4 : 6 : 1), (8 : 21 : 1), (11 : 35 : 1), (14 : 51 : 1), (21 : 95 : 1), (37 : 224 : 1), (52 : 374 : 1), (93 : 896 : 1), (342 : 6324 : 1), (406 : 8180 : 1), (816 : 23309 : 1)]
    52585262       
    52595263        ::
    52605264       
    5261             sage: a = E.integral_points([P1,P2,P3], verbose=True)
     5265            sage: a = E._integral_points_old([P1,P2,P3], verbose=True)
    52625266            Using mw_basis  [(2 : 0 : 1), (3 : -4 : 1), (8 : -22 : 1)]
    52635267            e1,e2,e3:  -3.0124303725933... 1.0658205476962... 1.94660982489710
    52645268            Minimal eigenvalue of height pairing matrix:  0.63792081458500...
     
    52785282        ::
    52795283       
    52805284            sage: E=EllipticCurve([0,0,1,-7,6])
    5281             sage: a=E.integral_points(both_signs=True); a
     5285            sage: a=E._integral_points_old(both_signs=True); a
    52825286            [(-3 : -1 : 1), (-3 : 0 : 1), (-2 : -4 : 1), (-2 : 3 : 1), (-1 : -4 : 1), (-1 : 3 : 1), (0 : -3 : 1), (0 : 2 : 1), (1 : -1 : 1), (1 : 0 : 1), (2 : -1 : 1), (2 : 0 : 1), (3 : -4 : 1), (3 : 3 : 1), (4 : -7 : 1), (4 : 6 : 1), (8 : -22 : 1), (8 : 21 : 1), (11 : -36 : 1), (11 : 35 : 1), (14 : -52 : 1), (14 : 51 : 1), (21 : -96 : 1), (21 : 95 : 1), (37 : -225 : 1), (37 : 224 : 1), (52 : -375 : 1), (52 : 374 : 1), (93 : -897 : 1), (93 : 896 : 1), (342 : -6325 : 1), (342 : 6324 : 1), (406 : -8181 : 1), (406 : 8180 : 1), (816 : -23310 : 1), (816 : 23309 : 1)]
    52835287       
    52845288        An example with negative discriminant::
    52855289       
    5286             sage: EllipticCurve('900d1').integral_points()
     5290            sage: EllipticCurve('900d1')._integral_points_old()
    52875291            [(-11 : 27 : 1), (-4 : 34 : 1), (4 : 18 : 1), (16 : 54 : 1)]
    52885292       
    52895293        Another example with rank 5 and no torsion points::
     
    52925296            sage: P1=E.point((540,1188)); P2=E.point((576,1836))
    52935297            sage: P3=E.point((468,3132)); P4=E.point((612,3132))
    52945298            sage: P5=E.point((432,4428))
    5295             sage: a=E.integral_points([P1,P2,P3,P4,P5]); len(a)  # long time (34s on sage.math, 2011)
     5299            sage: a=E._integral_points_old([P1,P2,P3,P4,P5]); len(a)  # long time (34s on sage.math, 2011)
    52965300            54
    52975301       
    52985302        TESTS:
    52995303
    53005304        The bug reported on trac #4525 is now fixed::
    53015305       
    5302             sage: EllipticCurve('91b1').integral_points()
     5306            sage: EllipticCurve('91b1')._integral_points_old()
    53035307            [(-1 : 3 : 1), (1 : 0 : 1), (3 : 4 : 1)]
    53045308       
    53055309        ::
    53065310       
    5307             sage: [len(e.integral_points(both_signs=False)) for e in cremona_curves([11..100])]  # long time (15s on sage.math, 2011)
     5311            sage: [len(e._integral_points_old(both_signs=False)) for e in cremona_curves([11..100])]  # long time (15s on sage.math, 2011)
    53085312            [2, 0, 2, 3, 2, 1, 3, 0, 2, 4, 2, 4, 3, 0, 0, 1, 2, 1, 2, 0, 2, 1, 0, 1, 3, 3, 1, 1, 4, 2, 3, 2, 0, 0, 5, 3, 2, 2, 1, 1, 1, 0, 1, 3, 0, 1, 0, 1, 1, 3, 6, 1, 2, 2, 2, 0, 0, 2, 3, 1, 2, 2, 1, 1, 0, 3, 2, 1, 0, 1, 0, 1, 3, 3, 1, 1, 5, 1, 0, 1, 1, 0, 1, 2, 0, 2, 0, 1, 1, 3, 1, 2, 2, 4, 4, 2, 1, 0, 0, 5, 1, 0, 1, 2, 0, 2, 2, 0, 0, 0, 1, 0, 3, 1, 5, 1, 2, 4, 1, 0, 1, 0, 1, 0, 1, 0, 2, 2, 0, 0, 1, 0, 1, 1, 4, 1, 0, 1, 1, 0, 4, 2, 0, 1, 1, 2, 3, 1, 1, 1, 1, 6, 2, 1, 1, 0, 2, 0, 6, 2, 0, 4, 2, 2, 0, 0, 1, 2, 0, 2, 1, 0, 3, 1, 2, 1, 4, 6, 3, 2, 1, 0, 2, 2, 0, 0, 5, 4, 1, 0, 0, 1, 0, 2, 2, 0, 0, 2, 3, 1, 3, 1, 1, 0, 1, 0, 0, 1, 2, 2, 0, 2, 0, 0, 1, 2, 0, 0, 4, 1, 0, 1, 1, 0, 1, 2, 0, 1, 4, 3, 1, 2, 2, 1, 1, 1, 1, 6, 3, 3, 3, 3, 1, 1, 1, 1, 1, 0, 7, 3, 0, 1, 3, 2, 1, 0, 3, 2, 1, 0, 2, 2, 6, 0, 0, 6, 2, 2, 3, 3, 5, 5, 1, 0, 6, 1, 0, 3, 1, 1, 2, 3, 1, 2, 1, 1, 0, 1, 0, 1, 0, 5, 5, 2, 2, 0, 0, 1, 0, 0, 0, 0, 1, 1, 1, 1]
    53095313       
    53105314        The bug reported at #4897 is now fixed::
    53115315
    5312             sage: [P[0] for P in EllipticCurve([0,0,0,-468,2592]).integral_points()]
     5316            sage: [P[0] for P in EllipticCurve([0,0,0,-468,2592])._integral_points_old()]
    53135317            [-24, -18, -14, -6, -3, 4, 6, 18, 21, 24, 36, 46, 102, 168, 186, 381, 1476, 2034, 67246]
    53145318
    53155319        .. note::
     
    53325336        #####################################################################
    53335337        # INPUT CHECK #######################################################
    53345338        if not self.is_integral():
    5335             raise ValueError, "integral_points() can only be called on an integral model"
     5339            raise ValueError, "_integral_points_old() can only be called on an integral model"
    53365340       
    53375341        if mw_base=='auto':
    53385342            mw_base = self.gens()
     
    56045608        if verbose:
    56055609            print 'Total number of integral points:',len(int_points)
    56065610        return int_points
    5607        
     5611
    56085612    def S_integral_points(self, S, mw_base='auto', both_signs=False, verbose=False, proof=None):
    56095613        """
    56105614        Computes all S-integral points (up to sign) on this elliptic curve.
     
    57585762        try:
    57595763            len_S = len(S)
    57605764            if len_S == 0:
    5761                 return self.integral_points(mw_base, both_signs, verbose)
     5765                return self.integral_points(L = mw_base,both_signs = both_signs)
    57625766            if not all([s.is_prime() for s in S]):
    57635767                raise ValueError, "All elements of S must be prime"
    57645768            S.sort()
     
    60226026        Qx = rings.PolynomialRing(RationalField(),'x')
    60236027        pol = Qx([-54*c6,-27*c4,0,1])
    60246028        if disc > 0: # two real component -> 3 roots in RR
    6025             # it is possible that only one root is found with default precision! (see integral_points())
     6029            # it is possible that only one root is found with default precision! (see _integral_points_old())
    60266030            RR = R
    60276031            prec = RR.precision()
    60286032            ei = pol.roots(RR,multiplicities=False)
     
    60706074        if verbose:
    60716075            print 'k1,k2,k3,k4',k1,k2,k3,k4
    60726076            sys.stdout.flush()
    6073         #H_q -> [PZGH]:N_0 (due to consistency to integral_points())
     6077        #H_q -> [PZGH]:N_0 (due to consistency to _integral_points_old())
    60746078        H_q = R(((k1/2+k2)/lamda).sqrt())
    60756079
    60766080        #computation of logs
     
    61496153            m_gram = m_LLL.gram_schmidt()[0]
    61506154            b1_norm = R(m_LLL.row(0).norm())
    61516155   
    6152             #compute constant c1_LLL (cf. integral_points())
     6156            #compute constant c1_LLL (cf. _integral_points_old())
    61536157            c1_LLL = -1
    61546158            for i in range(n):
    61556159                tmp = R(b1_norm/(m_gram.row(i).norm()))