Ticket #15299: 15299_lseries_prec.patch

File 15299_lseries_prec.patch, 38.5 KB (added by Jeroen Demeyer, 9 years ago)
  • sage/rings/rational_field.py

    # HG changeset patch
    # User Jeroen Demeyer <jdemeyer@cage.ugent.be>
    # Date 1382305549 -7200
    # Node ID 2f0adb62fcae0bd548fa817cb72b9e4902febc49
    # Parent  75deaaeb130a2179081f90496e0d4a66d05d2fe7
    Use arbitrary precision reals instead of Python floats in L series
    
    diff --git a/sage/rings/rational_field.py b/sage/rings/rational_field.py
    a b  
    121121        23.2000000000000
    122122        sage: QQ(a, 10)
    123123        116/5
    124    
     124
    125125    Here's a nice example involving elliptic curves::
    126    
     126
    127127        sage: E = EllipticCurve('11a')
    128128        sage: L = E.lseries().at1(300)[0]; L
    129         0.253841860855911
     129        0.2538418608559106843377589233...
    130130        sage: O = E.period_lattice().omega(); O
    131131        1.26920930427955
    132132        sage: t = L/O; t
  • sage/schemes/elliptic_curves/constructor.py

    diff --git a/sage/schemes/elliptic_curves/constructor.py b/sage/schemes/elliptic_curves/constructor.py
    a b  
    975975        sage: R.<x,y,z> = QQ[]
    976976        sage: C = Curve(x^3+y^3+z^3)
    977977        sage: P = C(1,-1,0)
    978         sage: E = EllipticCurve_from_plane_curve(C,P);  E
     978        sage: E = EllipticCurve_from_plane_curve(C,P); E  # long time (3s on sage.math, 2013)
    979979        doctest:...: DeprecationWarning: use Jacobian(C) instead
    980980        See http://trac.sagemath.org/3416 for details.
    981981        Elliptic Curve defined by y^2 = x^3 - 27/4 over Rational Field
  • 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  
    3030Using the "proof=False" flag suppresses these warnings.
    3131
    3232EXAMPLES: We find all elliptic curves with good reduction outside 2,
    33 listing the label of each:
     33listing the label of each::
    3434
    35     sage: [e.label() for e in EllipticCurves_with_good_reduction_outside_S([2])]
     35    sage: [e.label() for e in EllipticCurves_with_good_reduction_outside_S([2])]  # long time (5s on sage.math, 2013)
    3636    ['32a1',
    3737    '32a2',
    3838    '32a3',
     
    217217        []
    218218        sage: [e.label() for e in egros_from_j_0([3])]
    219219        ['27a1', '27a3', '243a1', '243a2', '243b1', '243b2']
    220         sage: len(egros_from_j_0([2,3,5]))
     220        sage: len(egros_from_j_0([2,3,5]))  # long time (8s on sage.math, 2013)
    221221        432
    222222    """
    223     Elist=[]   
     223    Elist=[]
    224224    if not 3 in S:
    225225        return Elist
    226226    no2 = not 2 in S
    227     for ei in xmrange([2] + [6]*len(S)): 
     227    for ei in xmrange([2] + [6]*len(S)):
    228228        u = prod([p**e for p,e in zip([-1]+S,ei)],QQ(1))
    229229        if no2:
    230230            u*=16 ## make sure 12|val(D,2)
     
    368368        sage: from sage.schemes.elliptic_curves.ell_egros import egros_get_j
    369369        sage: egros_get_j([])
    370370        [1728]
    371         sage: egros_get_j([2])
     371        sage: egros_get_j([2])  # long time (3s on sage.math, 2013)
    372372        [128, 432, -864, 1728, 3375/2, -3456, 6912, 8000, 10976, -35937/4, 287496, -784446336, -189613868625/128]
    373         sage: egros_get_j([3])
     373        sage: egros_get_j([3])  # long time (3s on sage.math, 2013)
    374374        [0, -576, 1536, 1728, -5184, -13824, 21952/9, -41472, 140608/3, -12288000]
    375375        sage: jlist=egros_get_j([2,3]); len(jlist) # long time (30s)
    376376        83
  • 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  
    267267
    268268        A curve with 2-torsion::
    269269
    270             sage: K.<a> = NumberField(x^2 + 7, 'a')
     270            sage: K.<a> = NumberField(x^2 + 7)
    271271            sage: E = EllipticCurve(K, '15a')
    272             sage: v = E.simon_two_descent(); v  # long time (about 10 seconds), points can vary
     272            sage: E.simon_two_descent()  # long time (3s on sage.math, 2013), points can vary
    273273            (1, 3, [...])
    274274
    275275        A failure in the PARI/GP script ell.gp (VERSION 25/03/2009) is reported::
    276            
     276
    277277            sage: K = CyclotomicField(43).subfields(3)[0][0]
    278278            sage: E = EllipticCurve(K, '37')
    279             sage: E.simon_two_descent()
     279            sage: E.simon_two_descent()  # long time (4s on sage.math, 2013)
    280280            Traceback (most recent call last):
    281281            ...
    282             RuntimeError: 
     282            RuntimeError:
    283283              ***   at top-level: ans=bnfellrank(K,[0,0,1,
    284284              ***                     ^--------------------
    285285              ***   in function bnfellrank: ...eqtheta,rnfeq,bbnf];rang=
     
    15001500
    15011501            sage: E = EllipticCurve('11a1')
    15021502            sage: K.<t>=NumberField(x^4 + x^3 + 11*x^2 + 41*x + 101)
    1503             sage: EK=E.base_extend(K)
    1504             sage: tor = EK.torsion_subgroup()
    1505             sage: tor
     1503            sage: EK = E.base_extend(K)
     1504            sage: tor = EK.torsion_subgroup()  # long time (3s on sage.math, 2013)
     1505            sage: tor  # long time
    15061506            Torsion Subgroup isomorphic to Z/5 + Z/5 associated to the Elliptic Curve defined by y^2 + y = x^3 + (-1)*x^2 + (-10)*x + (-20) over Number Field in t with defining polynomial x^4 + x^3 + 11*x^2 + 41*x + 101
    1507             sage: tor.gens()
     1507            sage: tor.gens()  # long time
    15081508            ((16 : 60 : 1), (t : 1/11*t^3 + 6/11*t^2 + 19/11*t + 48/11 : 1))
    15091509
    15101510        ::
     
    15961596            [(0 : 1 : 0), (5 : -6 : 1), (5 : 5 : 1), (16 : -61 : 1), (16 : 60 : 1)]
    15971597            sage: K.<t> = NumberField(x^4 + x^3 + 11*x^2 + 41*x + 101)
    15981598            sage: EK = E.base_extend(K)
    1599             sage: EK.torsion_points()
     1599            sage: EK.torsion_points()  # long time (3s on sage.math, 2013)
    16001600            [(16 : 60 : 1),
    16011601             (5 : 5 : 1),
    16021602             (5 : -6 : 1),
     
    21682168            Compatible family of Galois representations associated to the Elliptic Curve defined by y^2 + y = x^3 + (-1)*x^2 + (-10)*x + (-20) over Number Field in a with defining polynomial x^2 + 1
    21692169            sage: rho.is_surjective(3)
    21702170            True
    2171             sage: rho.is_surjective(5)
     2171            sage: rho.is_surjective(5)  # long time (9s on sage.math, 2013)
    21722172            False
    21732173            sage: rho.non_surjective()
    21742174            [5]
  • sage/schemes/elliptic_curves/gal_reps_number_field.py

    diff --git a/sage/schemes/elliptic_curves/gal_reps_number_field.py b/sage/schemes/elliptic_curves/gal_reps_number_field.py
    a b  
    88
    99EXAMPLES::
    1010
    11 sage: K = NumberField(x**2 - 29, 'a'); a = K.gen()
    12 sage: E = EllipticCurve([1, 0, ((5 + a)/2)**2, 0, 0])
    13 sage: rho = E.galois_representation()
    14 sage: rho.is_surjective(29) # Cyclotomic character not surjective.
    15 False
    16 sage: rho.is_surjective(31) # See Section 5.10 of [Serre72].
    17 True
    18 sage: rho.non_surjective()
    19 [3, 5, 29]
     11    sage: K = NumberField(x**2 - 29, 'a'); a = K.gen()
     12    sage: E = EllipticCurve([1, 0, ((5 + a)/2)**2, 0, 0])
     13    sage: rho = E.galois_representation()
     14    sage: rho.is_surjective(29) # Cyclotomic character not surjective.
     15    False
     16    sage: rho.is_surjective(31) # See Section 5.10 of [Serre72].
     17    True
     18    sage: rho.non_surjective()  # long time (8s on sage.math, 2013)
     19    [3, 5, 29]
    2020
    21 sage: E = EllipticCurve_from_j(1728).change_ring(K) # CM
    22 sage: E.galois_representation().non_surjective()
    23 [0]
     21    sage: E = EllipticCurve_from_j(1728).change_ring(K) # CM
     22    sage: E.galois_representation().non_surjective()
     23    [0]
    2424
    2525AUTHORS:
    2626
     
    174174            sage: K = NumberField(x**2 + 3, 'a'); a = K.gen()
    175175            sage: E = EllipticCurve([0, -1, 1, -10, -20]).change_ring(K) # X_0(11)
    176176            sage: rho = E.galois_representation()
    177             sage: rho.non_surjective()
     177            sage: rho.non_surjective()  # long time (8s on sage.math, 2013)
    178178            [3, 5]
    179179            sage: K = NumberField(x**2 + 1, 'a'); a = K.gen()
    180180            sage: E = EllipticCurve_from_j(1728).change_ring(K) # CM
  • sage/schemes/elliptic_curves/heegner.py

    diff --git a/sage/schemes/elliptic_curves/heegner.py b/sage/schemes/elliptic_curves/heegner.py
    a b  
    42444244        rank 1), and we reduce it modulo several primes.::
    42454245
    42464246            sage: E = EllipticCurve('11a1'); P = E.kolyvagin_point(-7)
    4247             sage: P.mod(3,70)
     4247            sage: P.mod(3,70)  # long time (4s on sage.math, 2013)
    42484248            (1 : 2 : 1)
    42494249            sage: P.mod(5,70)
    42504250            (1 : 4 : 1)
     
    53865386
    53875387            sage: N = 389; D = -7; ell = 5; c = 17; q = 3
    53885388            sage: H = heegner_points(N).reduce_mod(ell)
    5389             sage: k = H.rational_kolyvagin_divisor(D, c); k
     5389            sage: k = H.rational_kolyvagin_divisor(D, c); k  # long time (5s on sage.math, 2013)
    53905390            (14, 16, 0, 0, ... 0, 0, 0)
    53915391            sage: V = H.modp_dual_elliptic_curve_factor(EllipticCurve('389a'), q, 2)
    5392             sage: [b.dot_product(k.element().change_ring(GF(q))) for b in V.basis()]
     5392            sage: [b.dot_product(k.element().change_ring(GF(q))) for b in V.basis()]  # long time
    53935393            [0, 0]
    53945394            sage: k = H.rational_kolyvagin_divisor(D, 59)
    53955395            sage: [b.dot_product(k.element().change_ring(GF(q))) for b in V.basis()]
  • sage/schemes/elliptic_curves/lseries_ell.py

    diff --git a/sage/schemes/elliptic_curves/lseries_ell.py b/sage/schemes/elliptic_curves/lseries_ell.py
    a b  
    11"""
    22Complex Elliptic Curve L-series
    33
     4AUTHORS:
     5
     6- Jeroen Demeyer (2013-10-17): compute L series with arbitrary precision
     7  instead of floats.
     8
     9- William Stein et al. (2005 and later)
     10
    411"""
     12#*****************************************************************************
     13#       Copyright (C) 2005 William Stein
     14#       Copyright (C) 2013 Jeroen Demeyer
     15#
     16#  Distributed under the terms of the GNU General Public License (GPL)
     17#  as published by the Free Software Foundation; either version 2 of
     18#  the License, or (at your option) any later version.
     19#                  http://www.gnu.org/licenses/
     20#*****************************************************************************
    521
    622from sage.structure.sage_object import SageObject
    7 from sage.rings.all import (
    8     RealField,
    9     RationalField,
    10     ComplexField)
    11 from math import sqrt, exp, ceil
     23from sage.rings.all import RealField, RationalField
     24from math import sqrt, exp, log, ceil
    1225import sage.functions.exp_integral as exp_integral
    13 R = RealField()
    14 Q = RationalField()
    15 C = ComplexField()
    1626import sage.misc.all as misc
    1727
    1828class Lseries_ell(SageObject):
    1929    """
    2030    An elliptic curve $L$-series.
    21 
    22     EXAMPLES:
    23    
    2431    """
    2532    def __init__(self, E):
    2633        """
    2734        Create an elliptic curve $L$-series.
    2835
    29         EXAMPLES:
     36        EXAMPLES::
     37
    3038            sage: EllipticCurve([1..5]).lseries()
    3139            Complex L-series of the Elliptic Curve defined by y^2 + x*y + 3*y = x^3 + 2*x^2 + 4*x + 5 over Rational Field
    3240        """
     
    3543    def elliptic_curve(self):
    3644        """
    3745        Return the elliptic curve that this L-series is attached to.
    38        
    39         EXAMPLES:
     46
     47        EXAMPLES::
     48
    4049            sage: E = EllipticCurve('389a')
    4150            sage: L = E.lseries()
    4251            sage: L.elliptic_curve ()
     
    5261        The output is a series in var, where you should view var as
    5362        equal to s-a.  Thus this function returns the formal power
    5463        series whose coefficients are L^{(n)}(a)/n!.
    55        
    56         EXAMPLES:
     64
     65        EXAMPLES::
     66
    5767            sage: E = EllipticCurve('389a')
    5868            sage: L = E.lseries()
    59             sage: L.taylor_series(series_prec=3)      # random nearly 0 constant and linear terms
    60             -2.69129566562797e-23 + (1.52514901968783e-23)*z + 0.759316500288427*z^2 + O(z^3)
     69            sage: L.taylor_series(series_prec=3)
     70            -1.28158145675273e-23 + (7.26268290541182e-24)*z + 0.759316500288427*z^2 + O(z^3)  # 32-bit
     71            -2.69129566562797e-23 + (1.52514901968783e-23)*z + 0.759316500288427*z^2 + O(z^3)  # 64-bit
    6172            sage: L.taylor_series(series_prec=3)[2:]
    6273            0.000000000000000 + 0.000000000000000*z + 0.759316500288427*z^2 + O(z^3)
    6374        """
     
    6778    def _repr_(self):
    6879        """
    6980        Return string representation of this L-series.
    70        
    71         EXAMPLES:
     81
     82        EXAMPLES::
     83
    7284            sage: E = EllipticCurve('37a')
    7385            sage: L = E.lseries()
    7486            sage: L._repr_()
     
    96108        than bits and the object returned is a Magma L-series, which has
    97109        different functionality from the Sage L-series.}
    98110
    99         EXAMPLES:
     111        EXAMPLES::
     112
    100113            sage: E = EllipticCurve('37a')
    101114            sage: L = E.lseries().dokchitser()
    102115            sage: L(2)
     
    107120
    108121        If the curve has too large a conductor, it isn't possible to
    109122        compute with the L-series using this command.  Instead a
    110         RuntimeError is raised:
     123        RuntimeError is raised::
     124
    111125            sage: e = EllipticCurve([1,1,0,-63900,-1964465932632])
    112126            sage: L = e.lseries().dokchitser(15)
    113127            Traceback (most recent call last):
     
    159173        where \code{<n>} is replaced by your value of $n$.  This
    160174        command takes a long time to run.}
    161175
    162         EXAMPLES:
     176        EXAMPLES::
     177
    163178            sage: E = EllipticCurve('37a')
    164179            sage: a = E.lseries().sympow(2,16)   # not tested - requires precomputing "sympow('-new_data 2')"
    165180            sage: a                              # not tested
     
    188203        minutes.  If this function fails it will indicate what
    189204        commands have to be run.}
    190205
    191         EXAMPLES:
     206        EXAMPLES::
     207
    192208            sage: E = EllipticCurve('37a')
    193209            sage: print E.lseries().sympow_derivs(1,16,2)      # not tested -- requires precomputing "sympow('-new_data 2')"
    194210            sympow 1.018 RELEASE  (c) Mark Watkins --- see README and COPYING for details
     
    220236        Return the imaginary parts of the first $n$ nontrivial zeros
    221237        on the critical line of the L-function in the upper half
    222238        plane, as 32-bit reals.
    223        
    224         EXAMPLES:
     239
     240        EXAMPLES::
     241
    225242            sage: E = EllipticCurve('37a')
    226243            sage: E.lseries().zeros(2)
    227244              ***   Warning:...new stack size = ...
     
    255272        zeros, is missed). Higher up the critical strip you should use
    256273        a smaller stepsize so as not to miss zeros.
    257274
    258         EXAMPLES:
     275        EXAMPLES::
     276
    259277            sage: E = EllipticCurve('37a')
    260278            sage: E.lseries().zeros_in_interval(6, 10, 0.1)      # long time
    261279              ***   Warning:...new stack size = ...
     
    282300                    equally spaced sampled points on the line from
    283301                    s0 to s1.
    284302
    285         EXAMPLES:
    286             sage: I = CC.0
     303        EXAMPLES::
     304
    287305            sage: E = EllipticCurve('37a')
    288             sage: E.lseries().values_along_line(1, 0.5+20*I, 5)     # long time and slightly random output
    289             [(0.500000000, 0), (0.400000000 + 4.00000000*I, 3.31920245 - 2.60028054*I), (0.300000000 + 8.00000000*I, -0.886341185 - 0.422640337*I), (0.200000000 + 12.0000000*I, -3.50558936 - 0.108531690*I), (0.100000000 + 16.0000000*I, -3.87043288 - 1.88049411*I)]
     306            sage: E.lseries().values_along_line(1, 0.5 + 20*I, 5)
     307              ***   Warning:...new stack size = ...
     308            [(0.500000000, ...),
     309             (0.400000000 + 4.00000000*I, 3.31920245 - 2.60028054*I),
     310             (0.300000000 + 8.00000000*I, -0.886341185 - 0.422640337*I),
     311             (0.200000000 + 12.0000000*I, -3.50558936 - 0.108531690*I),
     312             (0.100000000 + 16.0000000*I, -3.87043288 - 1.88049411*I)]
     313
    290314        """
    291315        from sage.lfunctions.lcalc import lcalc
    292316        return lcalc.values_along_line(s0-RationalField()('1/2'),
     
    302326        critical strip is 1.}
    303327
    304328        INPUT:
    305             s -- complex numbers
    306             dmin -- integer
    307             dmax -- integer
     329
     330        - ``s`` -- complex numbers
     331        - ``dmin`` -- integer
     332        - ``dmax`` -- integer
    308333
    309334        OUTPUT:
    310             list -- list of pairs (d, L(E, s,chi_d))
    311335
    312         EXAMPLES:
     336        list of pairs (d, L(E, s,chi_d))
     337
     338        EXAMPLES::
     339
    313340            sage: E = EllipticCurve('37a')
    314             sage: E.lseries().twist_values(1, -12, -4)    # slightly random output depending on architecture
    315             [(-11, 1.4782434171), (-8, 0), (-7, 1.8530761916), (-4, 2.4513893817)]
     341            sage: vals = E.lseries().twist_values(1, -12, -4)
     342              ***   Warning:...new stack size = ...
     343            sage: vals  # abs tol 1e-17
     344            [(-11, 1.47824342), (-8, 8.9590946e-18), (-7, 1.85307619), (-4, 2.45138938)]
    316345            sage: F = E.quadratic_twist(-8)
    317346            sage: F.rank()
    318347            1
     
    341370            dict -- keys are the discriminants $d$, and
    342371                    values are list of corresponding zeros.
    343372
    344         EXAMPLES:
     373        EXAMPLES::
     374
    345375            sage: E = EllipticCurve('37a')
    346376            sage: E.lseries().twist_zeros(3, -4, -3)         # long time
    347377              ***   Warning:...new stack size = ...
     
    350380        from sage.lfunctions.lcalc import lcalc
    351381        return lcalc.twist_zeros(n, dmin, dmax, L=self.__E)
    352382
    353     def at1(self, k=0):
     383    def at1(self, k=None, prec=None):
    354384        r"""
    355         Compute $L(E,1)$ using $k$ terms of the series for $L(E,1)$ as
    356         explained on page 406 of Henri Cohen's book"A Course in Computational
    357         Algebraic Number Theory".  If the argument $k$ is not specified,
    358         then it defaults to $\sqrt(N)$, where $N$ is the conductor.
    359 
    360         The real precision used in each step of the computation is the
    361         precision of machine floats.
     385        Compute `L(E,1)` using `k` terms of the series for `L(E,1)` as
     386        explained on page 406 of Henri Cohen's book "A Course in Computational
     387        Algebraic Number Theory".  If the argument `k` is not specified,
     388        then it defaults to `\sqrt(N)`, where `N` is the conductor.
    362389
    363390        INPUT:
    364             k -- (optional) an integer, defaults to sqrt(N).
    365            
     391
     392        - ``k`` -- number of terms of the series. If zero or ``None``,
     393          use `k = \sqrt(N)`, where `N` is the conductor.
     394
     395        - ``prec`` -- numerical precision in bits. If zero or ``None``,
     396          use a reasonable automatic default.
     397
    366398        OUTPUT:
    367             float -- L(E,1)
    368             float -- a bound on the error in the approximation; this
    369                      is a provably correct upper bound on the sum
    370                      of the tail end of the series used to compute L(E,1).
    371399
    372         This function is disjoint from the PARI \code{elllseries}
     400        A tuple of real numbers ``(L, err)`` where ``L`` is an
     401        approximation for `L(E,1)` and ``err`` is a bound on the error
     402        in the approximation.
     403
     404        This function is disjoint from the PARI ``elllseries``
    373405        command, which is for a similar purpose.  To use that command
    374406        (via the PARI C library), simply type
    375                 \code{E.pari_mincurve().elllseries(1)}
     407        ``E.pari_mincurve().elllseries(1)``.
    376408
    377409        ALGORITHM:
    378         \begin{enumerate}
    379             \item Compute the root number eps.  If it is -1, return 0.
    380            
    381             \item Compute the Fourier coefficients a_n, for n up to and
    382                including k.
    383                
    384             \item Compute the sum
    385             $$
    386                  2 * sum_{n=1}^{k} (a_n / n) * exp(-2*pi*n/Sqrt(N)),
    387             $$
    388                where N is the conductor of E.
    389                
    390             \item Compute a bound on the tail end of the series, which is
    391             $$           
    392                  2 * e^(-2 * pi * (k+1) / sqrt(N)) / (1 - e^(-2*pi/sqrt(N))).
    393             $$     
    394                For a proof see [Grigov-Jorza-Patrascu-Patrikis-Stein].
    395         \end{enumerate}
    396410
    397         EXAMPLES:
     411        - Compute the root number eps.  If it is -1, return 0.
     412
     413        - Compute the Fourier coefficients `a_n`, for `n` up to and
     414          including `k`.
     415
     416        - Compute the sum
     417
     418          .. MATH::
     419
     420              2 * sum_{n=1}^{k} (a_n / n) * exp(-2*pi*n/Sqrt(N)),
     421
     422          where `N` is the conductor of `E`.
     423
     424        - Compute a bound on the tail end of the series, which is
     425
     426          .. MATH::
     427
     428                 2 e^{-2 \pi (k+1) / \sqrt{N}} / (1 - e^{-2 \pi/\sqrt{N}}).
     429
     430          For a proof see [Grigov-Jorza-Patrascu-Patrikis-Stein].
     431
     432        EXAMPLES::
     433
     434            sage: L, err = EllipticCurve('11a1').lseries().at1()
     435            sage: L, err
     436            (0.253804, 0.000181444)
     437            sage: parent(L)
     438            Real Field with 24 bits of precision
    398439            sage: E = EllipticCurve('37b')
     440            sage: E.lseries().at1()
     441            (0.7257177, 0.000800697)
    399442            sage: E.lseries().at1(100)
    400             (0.725681061936153, 1.52437502288743e-45)
     443            (0.7256810619361527823362055410263965487367603361763, 1.52469e-45)
     444            sage: L,err = E.lseries().at1(100, prec=128)
     445            sage: L
     446            0.72568106193615278233620554102639654873
     447            sage: parent(L)
     448            Real Field with 128 bits of precision
     449            sage: err
     450            1.70693e-37
     451            sage: parent(err)
     452            Real Field with 24 bits of precision and rounding RNDU
     453
     454        Rank 1 through 3 elliptic curves::
     455
     456            sage: E = EllipticCurve('37a1')
     457            sage: E.lseries().at1()
     458            (0.0000000, 0.000000)
     459            sage: E = EllipticCurve('389a1')
     460            sage: E.lseries().at1()
     461            (-0.001769566, 0.00911776)
     462            sage: E = EllipticCurve('5077a1')
     463            sage: E.lseries().at1()
     464            (0.0000000, 0.000000)
    401465        """
     466        sqrtN = sqrt(self.__E.conductor())
     467        if k:
     468            k = int(k)
     469        else:
     470            k = int(ceil(sqrtN))
     471
     472        if prec:
     473            prec = int(prec)
     474        else:
     475            # Use the same precision as deriv_at1() below for
     476            # consistency
     477            prec = int(9.065*k/sqrtN + 1.443*log(k)) + 12
     478        R = RealField(prec)
     479        # Compute error term with bounded precision of 24 bits and
     480        # round towards +infinity
     481        Rerror = RealField(24, rnd='RNDU')
     482
    402483        if self.__E.root_number() == -1:
    403             return 0
    404         sqrtN = float(self.__E.conductor().sqrt())
    405         k = int(k)
    406         if k == 0: k = int(ceil(sqrtN))
    407         an = self.__E.anlist(k)           # list of Sage ints
    408         # Compute z = e^(-2pi/sqrt(N))
    409         pi = 3.14159265358979323846
    410         z = exp(-2*pi/sqrtN)
     484           return (R.zero(), Rerror.zero())
     485
     486        an = self.__E.anlist(k)  # list of Sage Integers
     487        pi = R.pi()
     488        sqrtN = R(self.__E.conductor()).sqrt()
     489
     490        z = (-2*pi/sqrtN).exp()
    411491        zpow = z
    412         s = 0.0
     492        # Compute series sum and accumulate floating point errors
     493        L = R.zero()
     494        error = Rerror.zero()
     495
    413496        for n in xrange(1,k+1):
    414             s += (zpow * float(an[n]))/n
     497            term = (zpow * an[n])/n
    415498            zpow *= z
     499            L += term
     500            # We express relative error in units of epsilon, where
     501            # epsilon is a number divided by 2^precision.
     502            # Instead of multiplying the error by 2 after the loop
     503            # (to account for L *= 2), we already multiply it now.
     504            #
     505            # For multiplication and division, the relative error
     506            # in epsilons is bounded by (1+e)^n - 1, where n is the
     507            # number of operations (assuming exact inputs).
     508            # exp(x) additionally multiplies this error by abs(x) and
     509            # adds one epsilon. The inputs pi and sqrtN each contribute
     510            # another epsilon.
     511            # Assuming that 2*pi/sqrtN <= 2, the relative error for z is
     512            # 7 epsilon. This implies a relative error of (8n-1) epsilon
     513            # for zpow. We add 2 for the computation of term and 1/2 to
     514            # compensate for the approximation (1+e)^n = 1+ne.
     515            #
     516            # The error of the addition is at most half an ulp of the
     517            # result.
     518            #
     519            # Multiplying everything by two gives:
     520            error += term.epsilon(Rerror)*(16*n + 3) + L.ulp(Rerror)
     521        L *= 2
    416522
    417         error = 2*zpow / (1 - z)
    418        
    419         return R(2*s), R(error)
     523        # Add series error (we use (-2)/(z-1) instead of 2/(1-z)
     524        # because this causes 1/(1-z) to be rounded up)
     525        error += ((-2)*Rerror(zpow)) / Rerror(z - 1)
     526        return (L, error)
    420527
    421     def deriv_at1(self, k=0):
     528    def deriv_at1(self, k=None, prec=None):
    422529        r"""
    423         Compute $L'(E,1)$ using$ k$ terms of the series for $L'(E,1)$.
     530        Compute `L'(E,1)` using `k` terms of the series for `L'(E,1)`,
     531        under the assumption that `L(E,1) = 0`.
    424532
    425533        The algorithm used is from page 406 of Henri Cohen's book ``A
    426534        Course in Computational Algebraic Number Theory.''
    427535
    428         The real precision of the computation is the precision of
    429         Python floats.
     536        INPUT:
    430537
    431         INPUT:
    432             k -- int; number of terms of the series
     538        - ``k`` -- number of terms of the series. If zero or ``None``,
     539          use `k = \sqrt(N)`, where `N` is the conductor.
     540
     541        - ``prec`` -- numerical precision in bits. If zero or ``None``,
     542          use a reasonable automatic default.
    433543
    434544        OUTPUT:
    435             real number -- an approximation for L'(E,1)
    436             real number -- a bound on the error in the approximation
     545
     546        A tuple of real numbers ``(L1, err)`` where ``L1`` is an
     547        approximation for `L'(E,1)` and ``err`` is a bound on the error
     548        in the approximation.
     549
     550        .. WARNING::
     551
     552            This function only makes sense if `L(E)` has positive order
     553            of vanishing at 1, or equivalently if `L(E,1) = 0`.
    437554
    438555        ALGORITHM:
    439         \begin{enumerate}
    440             \item Compute the root number eps.  If it is 1, return 0.
    441556
    442             \item Compute the Fourier coefficients $a_n$, for $n$ up to and
    443                   including $k$.
    444                
    445             \item Compute the sum
    446                $$
     557        - Compute the root number eps.  If it is 1, return 0.
     558
     559        - Compute the Fourier coefficients `a_n`, for `n` up to and
     560          including `k`.
     561
     562        - Compute the sum
     563
     564          .. MATH::
     565
    447566                 2 * \sum_{n=1}^{k} (a_n / n) * E_1(2 \pi n/\sqrt{N}),
    448                $$ 
    449                where $N$ is the conductor of $E$, and $E_1$ is the
    450                exponential integral function.
    451567
    452             \item Compute a bound on the tail end of the series, which is
    453                $$
    454                  2 * e^{-2 \pi (k+1) / \sqrt{N}} / (1 - e^{-2 \ pi/\sqrt{N}}).
    455                $$ 
    456                For a proof see [Grigorov-Jorza-Patrascu-Patrikis-Stein].  This
    457                is exactly the same as the bound for the approximation to
    458                $L(E,1)$ produced by \code{E.lseries().at1}.
    459         \end{enumerate}
     568          where `N` is the conductor of `E`, and `E_1` is the
     569          exponential integral function.
     570
     571        - Compute a bound on the tail end of the series, which is
     572
     573          .. MATH::
     574
     575                 2 e^{-2 \pi (k+1) / \sqrt{N}} / (1 - e^{-2 \pi/\sqrt{N}}).
     576
     577          For a proof see [Grigorov-Jorza-Patrascu-Patrikis-Stein].  This
     578          is exactly the same as the bound for the approximation to
     579          `L(E,1)` produced by :meth:`at1`.
    460580
    461581        EXAMPLES::
    462582
    463583            sage: E = EllipticCurve('37a')
    464584            sage: E.lseries().deriv_at1()
    465             (0.305986660898516, 0.000800351433106958)
     585            (0.3059866, 0.000801045)
    466586            sage: E.lseries().deriv_at1(100)
    467             (0.305999773834052, 1.52437502288740e-45)
     587            (0.3059997738340523018204836833216764744526377745903, 1.52493e-45)
    468588            sage: E.lseries().deriv_at1(1000)
    469             (0.305999773834052, 0.000000000000000)
     589            (0.305999773834052301820483683321676474452637774590771998..., 2.75031e-449)
     590
     591        With less numerical precision, the error is bounded by numerical accuracy::
     592
     593            sage: L,err = E.lseries().deriv_at1(100, prec=64)
     594            sage: L,err
     595            (0.305999773834052302, 5.55318e-18)
     596            sage: parent(L)
     597            Real Field with 64 bits of precision
     598            sage: parent(err)
     599            Real Field with 24 bits of precision and rounding RNDU
     600
     601        Rank 2 and rank 3 elliptic curves::
     602
     603            sage: E = EllipticCurve('389a1')
     604            sage: E.lseries().deriv_at1()
     605            (0.0000000, 0.000000)
     606            sage: E = EllipticCurve((1, 0, 1, -131, 558))  # curve 59450i1
     607            sage: E.lseries().deriv_at1()
     608            (-0.00010911444, 0.142428)
     609            sage: E.lseries().deriv_at1(4000)
     610            (6.9902290...e-50, 1.31318e-43)
    470611        """
    471         if self.__E.root_number() == 1: return 0
    472         k = int(k)
    473         sqrtN = float(self.__E.conductor().sqrt())
    474         if k == 0: k = int(ceil(sqrtN))
    475         an = self.__E.anlist(k)           # list of Sage Integers
    476         # Compute z = e^(-2pi/sqrt(N))
    477         pi = 3.14159265358979323846
     612        sqrtN = sqrt(self.__E.conductor())
     613        if k:
     614            k = int(k)
     615        else:
     616            k = int(ceil(sqrtN))
     617
     618        if prec:
     619            prec = int(prec)
     620        else:
     621            # Estimate number of bits for the computation, based on error
     622            # estimate below (the denominator of that error is close enough
     623            # to 1 that we can ignore it).
     624            # 9.065 = 2*Pi/log(2)
     625            # 1.443 = 1/log(2)
     626            # 12 is an arbitrary extra number of bits (it is chosen
     627            #    such that the precision is 24 bits when the conductor
     628            #    equals 11 and k is the default value 4)
     629            prec = int(9.065*k/sqrtN + 1.443*log(k)) + 12
     630        R = RealField(prec)
     631        # Compute error term with bounded precision of 24 bits and
     632        # round towards +infinity
     633        Rerror = RealField(24, rnd='RNDU')
     634
     635        if self.__E.root_number() == 1:
     636           # Order of vanishing at 1 of L(E) is even and assumed to be
     637           # positive, so L'(E,1) = 0.
     638           return (R.zero(), Rerror.zero())
     639
     640        an = self.__E.anlist(k)  # list of Sage Integers
     641        pi = R.pi()
     642        sqrtN = R(self.__E.conductor()).sqrt()
    478643        v = exp_integral.exponential_integral_1(2*pi/sqrtN, k)
    479         L = 2*float(sum([ (v[n-1] * an[n])/n for n in xrange(1,k+1)]))
    480         error = 2*exp(-2*pi*(k+1)/sqrtN)/(1-exp(-2*pi/sqrtN))
    481         return R(L), R(error)
     644
     645        # Compute series sum and accumulate floating point errors
     646        L = R.zero()
     647        error = Rerror.zero()
     648        # Sum of |an[n]|/n
     649        sumann = Rerror.zero()
     650
     651        for n in xrange(1,k+1):
     652            term = (v[n-1] * an[n])/n
     653            L += term
     654            error += term.epsilon(Rerror)*5 + L.ulp(Rerror)
     655            sumann += Rerror(an[n].abs())/n
     656        L *= 2
     657
     658        # Add error term for exponential_integral_1() errors.
     659        # Absolute error for 2*v[i] is 4*max(1, v[0])*2^-prec
     660        if v[0] > 1.0:
     661            sumann *= Rerror(v[0])
     662        error += (sumann >> (prec - 2))
     663
     664        # Add series error (we use (-2)/(z-1) instead of 2/(1-z)
     665        # because this causes 1/(1-z) to be rounded up)
     666        z = (-2*pi/sqrtN).exp()
     667        zpow = ((-2*(k+1))*pi/sqrtN).exp()
     668        error += ((-2)*Rerror(zpow)) / Rerror(z - 1)
     669        return (L, error)
    482670
    483671    def __call__(self, s):
    484672        r"""
    485673        Returns the value of the L-series of the elliptic curve E at s, where s
    486674        must be a real number.
    487675
    488         Use self.extended for s complex.
     676        .. NOTE::
    489677
    490         \note{If the conductor of the curve is large, say $>10^{12}$,
    491         then this function will take a very long time, since it uses
    492         an $O(\sqrt{N})$ algorithm.}
     678            If the conductor of the curve is large, say `>10^{12}`,
     679            then this function will take a very long time, since it
     680            uses an `O(\sqrt{N})` algorithm.
    493681
    494         EXAMPLES:
     682        EXAMPLES::
     683
    495684            sage: E = EllipticCurve([1,2,3,4,5])
    496685            sage: L = E.lseries()
    497686            sage: L(1)
     
    503692        """
    504693        return self.dokchitser()(s)
    505694
    506     #def extended(self, s, prec):
    507     #    r"""
    508     #    Returns the value of the L-series of the elliptic curve E at s
    509     #    can be any complex number using prec terms of the power series
    510     #    expansion.
    511     #
    512     #
    513     #    WARNING: This may be slow.  Consider using \code{dokchitser()}
    514     #    instead.
    515     #   
    516     #    INPUT:
    517     #        s -- complex number
    518     #        prec -- integer
    519     #
    520     #    EXAMPLES:
    521     #        sage: E = EllipticCurve('389a')
    522     #        sage: E.lseries().extended(1 + I, 50)
    523     #        -0.638409959099589 + 0.715495262192901*I
    524     #        sage: E.lseries().extended(1 + 0.1*I, 50)
    525     #        -0.00761216538818315 + 0.000434885704670107*I
    526     #
    527     #    NOTE: You might also want to use Tim Dokchitser's
    528     #    L-function calculator, which is available by typing
    529     #    L = E.lseries().dokchitser(), then evaluating L.  It
    530     #    gives the same information but is sometimes much faster.
    531     #   
    532     #    """
    533     #    try:
    534     #        s = C(s)
    535     #    except TypeError:
    536     #        raise TypeError, "Input argument %s must be coercible to a complex number"%s
    537     #    prec = int(prec)
    538     #    if abs(s.imag()) < R(0.0000000000001):
    539     #        return self(s.real())
    540     #    N = self.__E.conductor()
    541     #    from sage.symbolic.constants import pi
    542     #    pi = R(pi)
    543     #    Gamma = transcendental.gamma
    544     #    Gamma_inc = transcendental.gamma_inc
    545     #    a = self.__E.anlist(prec)
    546     #    eps = self.__E.root_number()
    547     #    sqrtN = float(N.sqrt())
    548     #    def F(n, t):
    549     #        return Gamma_inc(t+1, 2*pi*n/sqrtN) * C(sqrtN/(2*pi*n))**(t+1)
    550     #    return C(N)**(-s/2) * C(2*pi)**s * Gamma(s)**(-1)\
    551     #           * sum([a[n]*(F(n,s-1) + eps*F(n,1-s)) for n in xrange(1,prec+1)])
    552 
    553695    def L1_vanishes(self):
    554696        """
    555         Returns whether or not L(E,1) = 0. The result is provably
     697        Returns whether or not `L(E,1) = 0`. The result is provably
    556698        correct if the Manin constant of the associated optimal
    557699        quotient is <= 2.  This hypothesis on the Manin constant
    558700        is true for all curves of conductor <= 40000 (by Cremona) and
    559701        all semistable curves (i.e., squarefree conductor).
    560702
    561         EXAMPLES:
     703        ALGORITHM: see :meth:`L_ratio`.
     704
     705        EXAMPLES::
     706
    562707            sage: E = EllipticCurve([0, -1, 1, -10, -20])   # 11A  = X_0(11)
    563708            sage: E.lseries().L1_vanishes()
    564709            False
     
    578723            sage: E.lseries().L1_vanishes()
    579724            False
    580725
    581         WARNING: It's conceivable that machine floats are not large
    582         enough precision for the computation; if this could be the
    583         case a RuntimeError is raised.  The curve's real period would
    584         have to be very small for this to occur. 
    585 
    586         ALGORITHM: Compute the root number.  If it is -1 then L(E,s)
    587         vanishes to odd order at 1, hence vanishes.  If it is +1, use
    588         a result about modular symbols and Mazur's "Rational Isogenies"
    589         paper to determine a provably correct bound (assuming Manin
    590         constant is <= 2) so that we can determine whether L(E,1) = 0.
    591 
    592726        AUTHOR: William Stein, 2005-04-20.
    593727        """
    594728        return self.L_ratio() == 0
    595729
    596730    def L_ratio(self):
    597731        r"""
    598         Returns the ratio $L(E,1)/\Omega$ as an exact rational
     732        Returns the ratio `L(E,1)/\Omega` as an exact rational
    599733        number. The result is \emph{provably} correct if the Manin
    600         constant of the associated optimal quotient is $\leq 2$.  This
     734        constant of the associated optimal quotient is `\leq 2`.  This
    601735        hypothesis on the Manin constant is true for all semistable
    602736        curves (i.e., squarefree conductor), by a theorem of Mazur
    603737        from his \emph{Rational Isogenies of Prime Degree} paper.
    604738
    605         EXAMPLES:
     739        EXAMPLES::
     740
    606741            sage: E = EllipticCurve([0, -1, 1, -10, -20])   # 11A  = X_0(11)
    607742            sage: E.lseries().L_ratio()
    608743            1/5
     
    625760            sage: E.lseries().L_ratio()
    626761            2
    627762
    628             # See trac #3651:
     763        See :trac:`3651` and :trac:`15299`::
     764
    629765            sage: EllipticCurve([0,0,0,-193^2,0]).sha().an()
    630766            4
    631 
    632         WARNING: It's conceivable that machine floats are not large
    633         enough precision for the computation; if this could be the
    634         case a RuntimeError is raised.  The curve's real period would
    635         have to be very small for this to occur.
     767            sage: EllipticCurve([1, 0, 1, -131, 558]).sha().an()  # long time
     768            1.00000000000000
    636769
    637770        ALGORITHM: Compute the root number.  If it is -1 then L(E,s)
    638771        vanishes to odd order at 1, hence vanishes.  If it is +1, use
     
    651784            self.__lratio = self.__E.minimal_model().lseries().L_ratio()
    652785            return self.__lratio
    653786
     787        QQ = RationalField()
    654788        if self.__E.root_number() == -1:
    655             self.__lratio = Q(0)
     789            self.__lratio = QQ.zero()
    656790            return self.__lratio
    657791
    658792        # Even root number.  Decide if L(E,1) = 0.  If E is a modular
     
    692826        d = self.__E._multiple_of_degree_of_isogeny_to_optimal_curve()
    693827        C = 8*d*t
    694828        eps = omega / C
    695         # coercion of 10**(-15) to our real field is needed to
    696         # make unambiguous comparison
    697         if eps < R(10**(-15)):  # liberal bound on precision of float
    698             raise RuntimeError, "Insufficient machine precision (=%s) for computation."%eps
     829
    699830        sqrtN = 2*int(sqrt(self.__E.conductor()))
    700831        k = sqrtN + 10
    701832        while True:
    702833            L1, error_bound = self.at1(k)
    703834            if error_bound < eps:
    704835                n = int(round(L1*C/omega))
    705                 quo = Q(n) / Q(C)
     836                quo = QQ((n,C))
    706837                self.__lratio = quo / self.__E.real_components()
    707838                return self.__lratio
    708839            k += sqrtN
  • sage/schemes/elliptic_curves/sha_tate.py

    diff --git a/sage/schemes/elliptic_curves/sha_tate.py b/sage/schemes/elliptic_curves/sha_tate.py
    a b  
    237237        See :trac:`1115`::
    238238
    239239            sage: sha=EllipticCurve('37a1').sha()
    240             sage: [sha.an_numerical(prec) for prec in xrange(40,100,10)]
     240            sage: [sha.an_numerical(prec) for prec in xrange(40,100,10)]  # long time (3s on sage.math, 2013)
    241241            [1.0000000000,
    242242            1.0000000000000,
    243243            1.0000000000000000,