Ticket #10973: trac_10973_v8.patch

File trac_10973_v8.patch, 51.5 KB (added by chapoton, 8 years ago)

replaces all previous patches ; rebased on 5.11

  • doc/en/reference/plane_curves/index.rst

    # HG changeset patch
    # User Author Justin C. Walker <justin@mac.com>
    # Date 1301015058 25200
    # Node ID 28f1b614c56e44ccb17ac2ad0aeef9c14135daf7
    # Parent  909180b6bbccd52134ebdcd60618a9cb6d10c926
    Trac #10973: integer points on elliptic curves over number fields
    
    diff --git a/doc/en/reference/plane_curves/index.rst b/doc/en/reference/plane_curves/index.rst
    a b Elliptic and Plane Curves 
    2626   sage/schemes/elliptic_curves/ell_number_field
    2727   sage/schemes/elliptic_curves/ell_finite_field
    2828   sage/schemes/elliptic_curves/ell_point
     29   sage/schemes/elliptic_curves/ell_int_points
    2930   sage/schemes/elliptic_curves/ell_torsion
    3031   sage/schemes/elliptic_curves/ell_local_data
    3132   sage/schemes/elliptic_curves/kodaira_symbol
  • sage/schemes/elliptic_curves/all.py

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