Ticket #10973: trac_10973.3.patch

File trac_10973.3.patch, 62.9 KB (added by Alyson Deines, 11 years ago)

replaces previous patch

  • sage/schemes/elliptic_curves/all.py

    # HG changeset patch
    # User Author Justin C. Walker <justin@mac.com>
    # Date 1301015058 25200
    # Node ID 177c64d0bb7f004009d713e49eb47864171d5449
    # Parent  b2c0ed2e09a453751e4b2fe092667487f530eac3
    Trac 10973: small bug fix to last patch, all tests should now pass
    
    diff -r b2c0ed2e09a4 -r 177c64d0bb7f 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# Temp until we integrate this more fully:
     45from ell_int_points import abs_log_height
  • new file sage/schemes/elliptic_curves/ell_int_points.py

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

    diff -r b2c0ed2e09a4 -r 177c64d0bb7f sage/schemes/elliptic_curves/ell_number_field.py
    a b  
    110110from sage.structure.element import RingElement
    111111from sage.rings.infinity import Infinity # just for verbose output
    112112from sage.rings.arith import valuation
     113from ell_int_points import integral_points
    113114
    114115class EllipticCurve_number_field(EllipticCurve_field):
    115116    r"""
     
    279280        prob_gens = [self(P) for P in t[2]]
    280281        self._simon_two_descent_data[lim1,lim3,limtriv,maxprob,limbigprime] = (prob_rank, two_selmer_rank, prob_gens)
    281282        return prob_rank, two_selmer_rank, prob_gens
     283       
     284    integral_points = integral_points
     285               
     286               
     287    def silverman_height_bounds(self):
     288        r"""
     289        Return the Silverman height bound.  This is a positive real
     290        (floating point) number B such that for all points `P` on the
     291        curve over any number field, `|h(P) - \hat{h}(P)| \leq B`,
     292        where `h(P)` is the naive logarithmic height of `P` and
     293        `\hat{h}(P)` is the canonical height.
     294
     295        NOTES:
     296
     297           - Silverman's paper is 'The Difference Between the Weil
     298             Height and the Canonical Height on Elliptic Curves',
     299             Math. Comp., Volume 55, Number 192, pages 723-743.  We
     300             use a correction by Bremner with 0.973 replaced by 0.961,
     301             as explained in the source code to mwrank (htconst.cc).
     302             
     303        EXAMPLES::
     304       
     305            sage: E=EllipticCurve('37a1')
     306           
     307                       
     308        """
     309        from sage.rings.all import CC
     310        Delta   = self.discriminant()
     311        j       = self.j_invariant()
     312        b2      = self.b2()
     313        K       = self.base_field()
     314        r,s     = K.signature()
     315        infplaces = K.embeddings(CC) if K is QQ else K.places()
     316        nvs = [1 if i < r else 2 for i,_ in enumerate(infplaces)]
     317       
     318        twostar = 2 if b2 else 1
     319        from math import log
     320        from ell_int_points import abs_log_height
     321        def h(x):
     322            return abs_log_height([K(x),K(1)])
     323        def h_oo(x):
     324            return sum(nvs[i]*log(max(abs(v(x)),1)) for i,v in enumerate(infplaces))/K.degree()
     325        mu    = h(Delta)/12 + h_oo(j)/12 + h_oo(b2/12)/2 + log(twostar)/2
     326        lower = 2*(-h(j)/24 - mu - 0.961)
     327        upper = 2*(mu + 1.07)
     328        return lower, upper
    282329               
    283330    def height_pairing_matrix(self, points=None, precision=None):
    284331        r"""
  • sage/schemes/elliptic_curves/ell_rational_field.py

    diff -r b2c0ed2e09a4 -r 177c64d0bb7f 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 
     31- Alyson Deines, R. Andrew Ohana, Jen Balakrishnan (2011): removed integral_points in
     32  preference for the more general number field method
    3033
    3134"""
    3235
     
    52215224
    52225225    prove_BSD = BSD.prove_BSD
    52235226
    5224     def integral_points(self, mw_base='auto', both_signs=False, verbose=False):
    5225         """
    5226         Computes all integral points (up to sign) on this elliptic curve.
    5227        
    5228         INPUT:
    5229        
    5230        
    5231         -  ``mw_base`` - list of EllipticCurvePoint generating
    5232            the Mordell-Weil group of E (default: 'auto' - calls self.gens())
    5233        
    5234         -  ``both_signs`` - True/False (default False): if
    5235            True the output contains both P and -P, otherwise only one of each
    5236            pair.
    5237        
    5238         -  ``verbose`` - True/False (default False): if True,
    5239            some details of the computation are output
    5240        
    5241        
    5242         OUTPUT: A sorted list of all the integral points on E (up to sign
    5243         unless both_signs is True)
    5244        
    5245         .. note::
    5246 
    5247            The complexity increases exponentially in the rank of curve
    5248            E. The computation time (but not the output!) depends on
    5249            the Mordell-Weil basis. If mw_base is given but is not a
    5250            basis for the Mordell-Weil group (modulo torsion), integral
    5251            points which are not in the subgroup generated by the given
    5252            points will almost certainly not be listed.
    5253        
    5254         EXAMPLES: A curve of rank 3 with no torsion points
    5255        
    5256         ::
    5257        
    5258             sage: E=EllipticCurve([0,0,1,-7,6])
    5259             sage: P1=E.point((2,0)); P2=E.point((-1,3)); P3=E.point((4,6))
    5260             sage: a=E.integral_points([P1,P2,P3]); a 
    5261             [(-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)]
    5262        
    5263         ::
    5264        
    5265             sage: a = E.integral_points([P1,P2,P3], verbose=True)
    5266             Using mw_basis  [(2 : 0 : 1), (3 : -4 : 1), (8 : -22 : 1)]
    5267             e1,e2,e3:  -3.0124303725933... 1.0658205476962... 1.94660982489710
    5268             Minimal eigenvalue of height pairing matrix:  0.63792081458500...
    5269             x-coords of points on compact component with  -3 <=x<= 1
    5270             [-3, -2, -1, 0, 1]
    5271             x-coords of points on non-compact component with  2 <=x<= 6
    5272             [2, 3, 4]
    5273             starting search of remaining points using coefficient bound  5
    5274             x-coords of extra integral points:
    5275             [2, 3, 4, 8, 11, 14, 21, 37, 52, 93, 342, 406, 816]
    5276             Total number of integral points: 18
    5277        
    5278         It is not necessary to specify mw_base; if it is not provided,
    5279         then the Mordell-Weil basis must be computed, which may take
    5280         much longer.
    5281 
    5282         ::
    5283        
    5284             sage: E=EllipticCurve([0,0,1,-7,6])
    5285             sage: a=E.integral_points(both_signs=True); a
    5286             [(-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)]
    5287        
    5288         An example with negative discriminant::
    5289        
    5290             sage: EllipticCurve('900d1').integral_points()
    5291             [(-11 : 27 : 1), (-4 : 34 : 1), (4 : 18 : 1), (16 : 54 : 1)]
    5292        
    5293         Another example with rank 5 and no torsion points::
    5294        
    5295             sage: E=EllipticCurve([-879984,319138704])
    5296             sage: P1=E.point((540,1188)); P2=E.point((576,1836))
    5297             sage: P3=E.point((468,3132)); P4=E.point((612,3132))
    5298             sage: P5=E.point((432,4428))
    5299             sage: a=E.integral_points([P1,P2,P3,P4,P5]); len(a)  # long time (34s on sage.math, 2011)
    5300             54
    5301        
    5302         TESTS:
    5303 
    5304         The bug reported on trac #4525 is now fixed::
    5305        
    5306             sage: EllipticCurve('91b1').integral_points()
    5307             [(-1 : 3 : 1), (1 : 0 : 1), (3 : 4 : 1)]
    5308        
    5309         ::
    5310        
    5311             sage: [len(e.integral_points(both_signs=False)) for e in cremona_curves([11..100])]  # long time (15s on sage.math, 2011)
    5312             [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]
    5313        
    5314         The bug reported at #4897 is now fixed::
    5315 
    5316             sage: [P[0] for P in EllipticCurve([0,0,0,-468,2592]).integral_points()]
    5317             [-24, -18, -14, -6, -3, 4, 6, 18, 21, 24, 36, 46, 102, 168, 186, 381, 1476, 2034, 67246]
    5318 
    5319         .. note::
    5320 
    5321            This function uses the algorithm given in [Co1].
    5322 
    5323         REFERENCES:
    5324 
    5325         - [Co1] Cohen H., Number Theory Vol I: Tools and Diophantine
    5326           Equations GTM 239, Springer 2007
    5327 
    5328         AUTHORS:
    5329 
    5330         - Michael Mardaus (2008-07)
    5331 
    5332         - Tobias Nagel (2008-07)
    5333 
    5334         - John Cremona (2008-07)
    5335         """
    5336         #####################################################################
    5337         # INPUT CHECK #######################################################
    5338         if not self.is_integral():
    5339             raise ValueError, "integral_points() can only be called on an integral model"
    5340        
    5341         if mw_base=='auto':
    5342             mw_base = self.gens()
    5343             r = len(mw_base)
    5344         else:
    5345             try:
    5346                 r = len(mw_base)
    5347             except TypeError:
    5348                 raise TypeError, 'mw_base must be a list'
    5349             if not all([P.curve() is self for P in mw_base]):
    5350                 raise ValueError, "points are not on the correct curve"
    5351                
    5352         tors_points = self.torsion_points()
    5353 
    5354         # END INPUT-CHECK####################################################
    5355         #####################################################################
    5356                
    5357         #####################################################################
    5358         # INTERNAL FUNCTIONS ################################################
    5359    
    5360         ############################## begin ################################
    5361         def point_preprocessing(free,tor):
    5362             r"""
    5363             Transforms the mw_basis ``free`` into a `\ZZ`-basis for
    5364             `E(\QQ)\cap E^0(`\RR)`. If there is a torsion point on the
    5365             "egg" we add it to any of the gens on the egg; otherwise
    5366             we replace the free generators with generators of a
    5367             subgroup of index 2.
    5368             """
    5369             r = len(free)
    5370             newfree = [Q for Q in free] # copy
    5371             tor_egg = [T for T in tor if not T.is_on_identity_component()]
    5372             free_id = [P.is_on_identity_component() for P in free]
    5373             if any(tor_egg):
    5374                 T = tor_egg[0]
    5375                 for i in range(r):
    5376                     if not free_id[i]:
    5377                         newfree[i] += T
    5378             else:
    5379                 if not all(free_id):
    5380                     i0 = free_id.index(False)
    5381                     P = free[i0]
    5382                     for i in range(r):
    5383                         if not free_id[i]:
    5384                             if i==i0:
    5385                                 newfree[i] = 2*newfree[i]
    5386                             else:
    5387                                 newfree[i] += P
    5388             return newfree
    5389         ##############################  end  ################################
    5390    
    5391         # END Internal functions #############################################
    5392         ######################################################################
    5393    
    5394         if (r == 0):
    5395             int_points = [P for P in tors_points if not P.is_zero()]
    5396             int_points = [P for P in int_points if P[0].is_integral()]
    5397             if not both_signs:
    5398                 xlist = set([P[0] for P in int_points])
    5399                 int_points = [self.lift_x(x) for x in xlist]
    5400             int_points.sort()
    5401             if verbose:
    5402                 print 'Total number of integral points:',len(int_points)
    5403             return int_points
    5404 
    5405         if verbose:
    5406             import sys  # so we can flush stdout for debugging
    5407        
    5408         g2 = self.c4()/12
    5409         g3 = self.c6()/216
    5410         disc = self.discriminant()
    5411         j = self.j_invariant()
    5412         b2 = self.b2()
    5413        
    5414         Qx = rings.PolynomialRing(RationalField(),'x')
    5415         pol = Qx([-self.c6()/216,-self.c4()/12,0,4])
    5416         if disc > 0: # two real component -> 3 roots in RR
    5417             #on curve 897e4, only one root is found with default precision!
    5418             RR = R
    5419             prec = RR.precision()
    5420             ei = pol.roots(RR,multiplicities=False)
    5421             while len(ei)<3:
    5422                 prec*=2
    5423                 RR=RealField(prec)
    5424                 ei = pol.roots(RR,multiplicities=False)
    5425             e1,e2,e3 = ei
    5426             if r >= 1: #preprocessing of mw_base only necessary if rank > 0
    5427                 mw_base = point_preprocessing(mw_base, tors_points)
    5428                   #at most one point in E^{egg}
    5429 
    5430         elif disc < 0: # one real component => 1 root in RR (=: e3),
    5431                        # 2 roots in C (e1,e2)
    5432             roots = pol.roots(C,multiplicities=False)
    5433             e3 = pol.roots(R,multiplicities=False)[0]
    5434             roots.remove(e3)
    5435             e1,e2 = roots
    5436 
    5437         from sage.all import pi
    5438         e = R(1).exp()
    5439         pi = R(pi)
    5440 
    5441         M = self.height_pairing_matrix(mw_base)
    5442         mw_base, U = self.lll_reduce(mw_base,M)
    5443         M = U.transpose()*M*U
    5444        
    5445         if verbose:
    5446             print "Using mw_basis ",mw_base
    5447             print "e1,e2,e3: ",e1,e2,e3
    5448             sys.stdout.flush()
    5449            
    5450         # Algorithm presented in [Co1]
    5451         h_E = self.height()
    5452         w1, w2 = self.period_lattice().basis()
    5453         mu = R(disc).abs().log() / 6
    5454         if j!=0:
    5455             mu += max(R(1),R(j).abs().log()) / 6
    5456         if b2!=0:
    5457             mu += max(R(1),R(b2).abs().log())
    5458             mu += log(R(2))
    5459         else:
    5460             mu += 1
    5461        
    5462         c1 = (mu + 2.14).exp()
    5463         c2 = min(M.charpoly ().roots(multiplicities=False))
    5464         if verbose:
    5465             print "Minimal eigenvalue of height pairing matrix: ", c2
    5466             sys.stdout.flush()
    5467 
    5468         c3 = (w1**2)*R(b2).abs()/48 + 8
    5469         c5 = (c1*c3).sqrt()
    5470         c7 = abs((3*pi)/((w1**2) * (w1/w2).imag()))
    5471        
    5472         mw_base_log = [] #contains \Phi(Q_i)
    5473         mod_h_list = []  #contains h_m(Q_i)
    5474         c9_help_list = []
    5475         for i in range(0,r):
    5476             mw_base_log.append(mw_base[i].elliptic_logarithm().abs())
    5477             mod_h_list.append(max(mw_base[i].height(),h_E,c7*mw_base_log[i]**2))
    5478             c9_help_list.append((mod_h_list[i]).sqrt()/mw_base_log[i])
    5479         c8 = max(e*h_E,max(mod_h_list))
    5480         c9 = e/c7.sqrt() * min(c9_help_list)
    5481         n=r+1
    5482         c10 = R(2 * 10**(8+7*n) * R((2/e)**(2 * n**2)) * (n+1)**(4 * n**2 + 10 * n) * log(c9)**(-2*n - 1) * misc.prod(mod_h_list))
    5483    
    5484         top = Z(128) #arbitrary first upper bound
    5485         bottom = Z(0)
    5486         log_c9=log(c9); log_c5=log(c5)
    5487         log_r_top = log(R(r*(10**top)))
    5488 #        if verbose:
    5489 #            print "[bottom,top] = ",[bottom,top]
    5490            
    5491         while R(c10*(log_r_top+log_c9)*(log(log_r_top)+h_E+log_c9)**(n+1)) > R(c2/2 * (10**top)**2 - log_c5):
    5492             #initial bound 'top' too small, upshift of search interval
    5493             bottom = top
    5494             top = 2*top
    5495         while top >= bottom: #binary-search like search for fitting exponent (bound)
    5496 #            if verbose:
    5497 #                print "[bottom,top] = ",[bottom,top]
    5498             bound = (bottom + (top - bottom)/2).floor()
    5499             log_r_bound = log(R(r*(10**bound)))
    5500             if R(c10*(log_r_bound+log_c9)*(log(log_r_bound)+h_E+log_c9)**(n+1)) > R(c2/2 * (10**bound)**2 - log_c5):
    5501                 bottom = bound + 1
    5502             else:
    5503                 top = bound - 1
    5504    
    5505         H_q = R(10)**bound
    5506         break_cond = 0 #at least one reduction step
    5507         #reduction via LLL
    5508         while break_cond < 0.9: #as long as the improvement of the new bound in comparison to the old is greater than 10%
    5509             c = R((H_q**n)*10)  #c has to be greater than H_q^n
    5510             M = matrix.MatrixSpace(Z,n)
    5511             m = copy(M.identity_matrix())
    5512             for i in range(r):
    5513                 m[i, r] = R(c*mw_base_log[i]).round()
    5514             m[r,r] = max(Z(1),R(c*w1).round()) #ensures that m isn't singular
    5515    
    5516             #LLL - implemented in sage - operates on rows not on columns
    5517             m_LLL = m.LLL()
    5518             m_gram = m_LLL.gram_schmidt()[0]
    5519             b1_norm = R(m_LLL.row(0).norm())
    5520    
    5521             #compute constant c1 ~ c1_LLL of Corollary 2.3.17 and hence d(L,0)^2 ~ d_L_0
    5522             c1_LLL = -1
    5523             for i in range(n):
    5524                 tmp = R(b1_norm/(m_gram.row(i).norm()))
    5525                 if tmp > c1_LLL:
    5526                     c1_LLL = tmp
    5527    
    5528             if c1_LLL < 0:
    5529                 raise RuntimeError, 'Unexpected intermediate result. Please try another Mordell-Weil base'
    5530            
    5531             d_L_0 = R(b1_norm**2 / c1_LLL)
    5532    
    5533             #Reducing of upper bound
    5534             Q = r * H_q**2
    5535             T = (1 + (3/2*r*H_q))/2
    5536             if d_L_0 < R(T**2+Q):
    5537                 d_L_0 = 10*(T**2*Q)
    5538             low_bound = R(((d_L_0 - Q).sqrt() - T)/c)
    5539    
    5540             #new bound according to low_bound and upper bound
    5541             #[c_5 exp((-c_2*H_q^2)/2)] provided by Corollary 8.7.3
    5542             if low_bound != 0:
    5543                 H_q_new = R((log(low_bound/c5)/(-c2/2))).sqrt()
    5544                 H_q_new = H_q_new.ceil()
    5545                 if H_q_new == 1:
    5546                     break_cond = 1 # stops reduction
    5547                 else:
    5548                     break_cond = R(H_q_new/H_q)
    5549                 H_q = H_q_new
    5550             else:
    5551                 break_cond = 1 # stops reduction, so using last H_q > 0
    5552             #END LLL-Reduction loop
    5553    
    5554         b2_12 = b2/12
    5555         if disc > 0:
    5556             ##Points in egg have X(P) between e1 and e2 [X(P)=x(P)+b2/12]:
    5557             x_int_points = self.integral_x_coords_in_interval((e1-b2_12).ceil(), (e2-b2_12).floor())
    5558             if verbose:
    5559                 print 'x-coords of points on compact component with ',(e1-b2_12).ceil(),'<=x<=',(e2-b2_12).floor()
    5560                 L = list(x_int_points) # to have the order
    5561                 L.sort()               # deterministic for doctests!
    5562                 print L
    5563                 sys.stdout.flush()
    5564         else:
    5565             x_int_points = set()
    5566            
    5567         ##Points in noncompact component with X(P)< 2*max(|e1|,|e2|,|e3|) , espec. X(P)>=e3
    5568         x0 = (e3-b2_12).ceil()
    5569         x1 = (2*max(abs(e1),abs(e2),abs(e3)) - b2_12).ceil()
    5570         x_int_points2 = self.integral_x_coords_in_interval(x0, x1)
    5571         x_int_points = x_int_points.union(x_int_points2)
    5572         if verbose:
    5573             print 'x-coords of points on non-compact component with ',x0,'<=x<=',x1-1
    5574             L = list(x_int_points2)
    5575             L.sort()
    5576             print L
    5577             sys.stdout.flush()
    5578        
    5579         if verbose:
    5580             print 'starting search of remaining points using coefficient bound ',H_q
    5581             sys.stdout.flush()
    5582         x_int_points3 = integral_points_with_bounded_mw_coeffs(self,mw_base,H_q)
    5583         x_int_points = x_int_points.union(x_int_points3)
    5584         if verbose:
    5585             print 'x-coords of extra integral points:'
    5586             L = list(x_int_points3)
    5587             L.sort()
    5588             print L
    5589             sys.stdout.flush()
    5590 
    5591         if len(tors_points)>1:
    5592             x_int_points_t = set()
    5593             for x in x_int_points:
    5594                 P = self.lift_x(x)
    5595                 for T in tors_points:
    5596                     Q = P+T
    5597                     if not Q.is_zero() and Q[0].is_integral():
    5598                         x_int_points_t = x_int_points_t.union([Q[0]])
    5599             x_int_points = x_int_points.union(x_int_points_t)
    5600            
    5601         # Now we have all the x-coordinates of integral points, and we
    5602         # construct the points, depending on the parameter both_signs:
    5603         if both_signs:
    5604             int_points = sum([self.lift_x(x,all=True) for x in x_int_points],[])
    5605         else:
    5606             int_points = [self.lift_x(x) for x in x_int_points]
    5607         int_points.sort()
    5608         if verbose:
    5609             print 'Total number of integral points:',len(int_points)
    5610         return int_points
    5611        
    56125227    def S_integral_points(self, S, mw_base='auto', both_signs=False, verbose=False, proof=None):
    56135228        """
    56145229        Computes all S-integral points (up to sign) on this elliptic curve.
     
    57625377        try:
    57635378            len_S = len(S)
    57645379            if len_S == 0:
    5765                 return self.integral_points(mw_base, both_signs, verbose)
     5380                if mw_base == 'auto':
     5381                    mw_base = self.gens(proof = proof)
     5382                return self.integral_points(L = mw_base, both_signs = both_signs)
    57665383            if not all([s.is_prime() for s in S]):
    57675384                raise ValueError, "All elements of S must be prime"
    57685385            S.sort()