# 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
# 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 23.2000000000000 sage: QQ(a, 10) 116/5 Here's a nice example involving elliptic curves:: sage: E = EllipticCurve('11a') sage: L = E.lseries().at1(300)[0]; L 0.253841860855911 0.2538418608559106843377589233... sage: O = E.period_lattice().omega(); O 1.26920930427955 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 sage: R. = QQ[] sage: C = Curve(x^3+y^3+z^3) sage: P = C(1,-1,0) sage: E = EllipticCurve_from_plane_curve(C,P);  E sage: E = EllipticCurve_from_plane_curve(C,P); E  # long time (3s on sage.math, 2013) doctest:...: DeprecationWarning: use Jacobian(C) instead See http://trac.sagemath.org/3416 for details. 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 Using the "proof=False" flag suppresses these warnings. EXAMPLES: We find all elliptic curves with good reduction outside 2, listing the label of each: listing the label of each:: sage: [e.label() for e in EllipticCurves_with_good_reduction_outside_S([2])] sage: [e.label() for e in EllipticCurves_with_good_reduction_outside_S([2])]  # long time (5s on sage.math, 2013) ['32a1', '32a2', '32a3', [] sage: [e.label() for e in egros_from_j_0([3])] ['27a1', '27a3', '243a1', '243a2', '243b1', '243b2'] sage: len(egros_from_j_0([2,3,5])) sage: len(egros_from_j_0([2,3,5]))  # long time (8s on sage.math, 2013) 432 """ Elist=[] Elist=[] if not 3 in S: return Elist no2 = not 2 in S for ei in xmrange([2] + [6]*len(S)): for ei in xmrange([2] + [6]*len(S)): u = prod([p**e for p,e in zip([-1]+S,ei)],QQ(1)) if no2: u*=16 ## make sure 12|val(D,2) sage: from sage.schemes.elliptic_curves.ell_egros import egros_get_j sage: egros_get_j([]) [1728] sage: egros_get_j([2]) sage: egros_get_j([2])  # long time (3s on sage.math, 2013) [128, 432, -864, 1728, 3375/2, -3456, 6912, 8000, 10976, -35937/4, 287496, -784446336, -189613868625/128] sage: egros_get_j([3]) sage: egros_get_j([3])  # long time (3s on sage.math, 2013) [0, -576, 1536, 1728, -5184, -13824, 21952/9, -41472, 140608/3, -12288000] sage: jlist=egros_get_j([2,3]); len(jlist) # long time (30s) 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 A curve with 2-torsion:: sage: K. = NumberField(x^2 + 7, 'a') sage: K. = NumberField(x^2 + 7) sage: E = EllipticCurve(K, '15a') sage: v = E.simon_two_descent(); v  # long time (about 10 seconds), points can vary sage: E.simon_two_descent()  # long time (3s on sage.math, 2013), points can vary (1, 3, [...]) A failure in the PARI/GP script ell.gp (VERSION 25/03/2009) is reported:: sage: K = CyclotomicField(43).subfields(3)[0][0] sage: E = EllipticCurve(K, '37') sage: E.simon_two_descent() sage: E.simon_two_descent()  # long time (4s on sage.math, 2013) Traceback (most recent call last): ... RuntimeError: RuntimeError: ***   at top-level: ans=bnfellrank(K,[0,0,1, ***                     ^-------------------- ***   in function bnfellrank: ...eqtheta,rnfeq,bbnf];rang= sage: E = EllipticCurve('11a1') sage: K.=NumberField(x^4 + x^3 + 11*x^2 + 41*x + 101) sage: EK=E.base_extend(K) sage: tor = EK.torsion_subgroup() sage: tor sage: EK = E.base_extend(K) sage: tor = EK.torsion_subgroup()  # long time (3s on sage.math, 2013) sage: tor  # long time 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 sage: tor.gens() sage: tor.gens()  # long time ((16 : 60 : 1), (t : 1/11*t^3 + 6/11*t^2 + 19/11*t + 48/11 : 1)) :: [(0 : 1 : 0), (5 : -6 : 1), (5 : 5 : 1), (16 : -61 : 1), (16 : 60 : 1)] sage: K. = NumberField(x^4 + x^3 + 11*x^2 + 41*x + 101) sage: EK = E.base_extend(K) sage: EK.torsion_points() sage: EK.torsion_points()  # long time (3s on sage.math, 2013) [(16 : 60 : 1), (5 : 5 : 1), (5 : -6 : 1), 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 sage: rho.is_surjective(3) True sage: rho.is_surjective(5) sage: rho.is_surjective(5)  # long time (9s on sage.math, 2013) False sage: rho.non_surjective() [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 EXAMPLES:: sage: K = NumberField(x**2 - 29, 'a'); a = K.gen() sage: E = EllipticCurve([1, 0, ((5 + a)/2)**2, 0, 0]) sage: rho = E.galois_representation() sage: rho.is_surjective(29) # Cyclotomic character not surjective. False sage: rho.is_surjective(31) # See Section 5.10 of [Serre72]. True sage: rho.non_surjective() [3, 5, 29] sage: K = NumberField(x**2 - 29, 'a'); a = K.gen() sage: E = EllipticCurve([1, 0, ((5 + a)/2)**2, 0, 0]) sage: rho = E.galois_representation() sage: rho.is_surjective(29) # Cyclotomic character not surjective. False sage: rho.is_surjective(31) # See Section 5.10 of [Serre72]. True sage: rho.non_surjective()  # long time (8s on sage.math, 2013) [3, 5, 29] sage: E = EllipticCurve_from_j(1728).change_ring(K) # CM sage: E.galois_representation().non_surjective() [0] sage: E = EllipticCurve_from_j(1728).change_ring(K) # CM sage: E.galois_representation().non_surjective() [0] AUTHORS: sage: K = NumberField(x**2 + 3, 'a'); a = K.gen() sage: E = EllipticCurve([0, -1, 1, -10, -20]).change_ring(K) # X_0(11) sage: rho = E.galois_representation() sage: rho.non_surjective() sage: rho.non_surjective()  # long time (8s on sage.math, 2013) [3, 5] sage: K = NumberField(x**2 + 1, 'a'); a = K.gen() 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 rank 1), and we reduce it modulo several primes.:: sage: E = EllipticCurve('11a1'); P = E.kolyvagin_point(-7) sage: P.mod(3,70) sage: P.mod(3,70)  # long time (4s on sage.math, 2013) (1 : 2 : 1) sage: P.mod(5,70) (1 : 4 : 1) sage: N = 389; D = -7; ell = 5; c = 17; q = 3 sage: H = heegner_points(N).reduce_mod(ell) sage: k = H.rational_kolyvagin_divisor(D, c); k sage: k = H.rational_kolyvagin_divisor(D, c); k  # long time (5s on sage.math, 2013) (14, 16, 0, 0, ... 0, 0, 0) sage: V = H.modp_dual_elliptic_curve_factor(EllipticCurve('389a'), q, 2) sage: [b.dot_product(k.element().change_ring(GF(q))) for b in V.basis()] sage: [b.dot_product(k.element().change_ring(GF(q))) for b in V.basis()]  # long time [0, 0] sage: k = H.rational_kolyvagin_divisor(D, 59) 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 """ Complex Elliptic Curve L-series AUTHORS: - Jeroen Demeyer (2013-10-17): compute L series with arbitrary precision instead of floats. - William Stein et al. (2005 and later) """ #***************************************************************************** #       Copyright (C) 2005 William Stein #       Copyright (C) 2013 Jeroen Demeyer # #  Distributed under the terms of the GNU General Public License (GPL) #  as published by the Free Software Foundation; either version 2 of #  the License, or (at your option) any later version. #                  http://www.gnu.org/licenses/ #***************************************************************************** from sage.structure.sage_object import SageObject from sage.rings.all import ( RealField, RationalField, ComplexField) from math import sqrt, exp, ceil from sage.rings.all import RealField, RationalField from math import sqrt, exp, log, ceil import sage.functions.exp_integral as exp_integral R = RealField() Q = RationalField() C = ComplexField() import sage.misc.all as misc class Lseries_ell(SageObject): """ An elliptic curve $L$-series. EXAMPLES: """ def __init__(self, E): """ Create an elliptic curve $L$-series. EXAMPLES: EXAMPLES:: sage: EllipticCurve([1..5]).lseries() 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 """ def elliptic_curve(self): """ Return the elliptic curve that this L-series is attached to. EXAMPLES: EXAMPLES:: sage: E = EllipticCurve('389a') sage: L = E.lseries() sage: L.elliptic_curve () The output is a series in var, where you should view var as equal to s-a.  Thus this function returns the formal power series whose coefficients are L^{(n)}(a)/n!. EXAMPLES: EXAMPLES:: sage: E = EllipticCurve('389a') sage: L = E.lseries() sage: L.taylor_series(series_prec=3)      # random nearly 0 constant and linear terms -2.69129566562797e-23 + (1.52514901968783e-23)*z + 0.759316500288427*z^2 + O(z^3) sage: L.taylor_series(series_prec=3) -1.28158145675273e-23 + (7.26268290541182e-24)*z + 0.759316500288427*z^2 + O(z^3)  # 32-bit -2.69129566562797e-23 + (1.52514901968783e-23)*z + 0.759316500288427*z^2 + O(z^3)  # 64-bit sage: L.taylor_series(series_prec=3)[2:] 0.000000000000000 + 0.000000000000000*z + 0.759316500288427*z^2 + O(z^3) """ def _repr_(self): """ Return string representation of this L-series. EXAMPLES: EXAMPLES:: sage: E = EllipticCurve('37a') sage: L = E.lseries() sage: L._repr_() than bits and the object returned is a Magma L-series, which has different functionality from the Sage L-series.} EXAMPLES: EXAMPLES:: sage: E = EllipticCurve('37a') sage: L = E.lseries().dokchitser() sage: L(2) If the curve has too large a conductor, it isn't possible to compute with the L-series using this command.  Instead a RuntimeError is raised: RuntimeError is raised:: sage: e = EllipticCurve([1,1,0,-63900,-1964465932632]) sage: L = e.lseries().dokchitser(15) Traceback (most recent call last): where \code{} is replaced by your value of $n$.  This command takes a long time to run.} EXAMPLES: EXAMPLES:: sage: E = EllipticCurve('37a') sage: a = E.lseries().sympow(2,16)   # not tested - requires precomputing "sympow('-new_data 2')" sage: a                              # not tested minutes.  If this function fails it will indicate what commands have to be run.} EXAMPLES: EXAMPLES:: sage: E = EllipticCurve('37a') sage: print E.lseries().sympow_derivs(1,16,2)      # not tested -- requires precomputing "sympow('-new_data 2')" sympow 1.018 RELEASE  (c) Mark Watkins --- see README and COPYING for details Return the imaginary parts of the first $n$ nontrivial zeros on the critical line of the L-function in the upper half plane, as 32-bit reals. EXAMPLES: EXAMPLES:: sage: E = EllipticCurve('37a') sage: E.lseries().zeros(2) ***   Warning:...new stack size = ... zeros, is missed). Higher up the critical strip you should use a smaller stepsize so as not to miss zeros. EXAMPLES: EXAMPLES:: sage: E = EllipticCurve('37a') sage: E.lseries().zeros_in_interval(6, 10, 0.1)      # long time ***   Warning:...new stack size = ... equally spaced sampled points on the line from s0 to s1. EXAMPLES: sage: I = CC.0 EXAMPLES:: sage: E = EllipticCurve('37a') sage: E.lseries().values_along_line(1, 0.5+20*I, 5)     # long time and slightly random output [(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)] sage: E.lseries().values_along_line(1, 0.5 + 20*I, 5) ***   Warning:...new stack size = ... [(0.500000000, ...), (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)] """ from sage.lfunctions.lcalc import lcalc return lcalc.values_along_line(s0-RationalField()('1/2'), critical strip is 1.} INPUT: s -- complex numbers dmin -- integer dmax -- integer - s -- complex numbers - dmin -- integer - dmax -- integer OUTPUT: list -- list of pairs (d, L(E, s,chi_d)) EXAMPLES: list of pairs (d, L(E, s,chi_d)) EXAMPLES:: sage: E = EllipticCurve('37a') sage: E.lseries().twist_values(1, -12, -4)    # slightly random output depending on architecture [(-11, 1.4782434171), (-8, 0), (-7, 1.8530761916), (-4, 2.4513893817)] sage: vals = E.lseries().twist_values(1, -12, -4) ***   Warning:...new stack size = ... sage: vals  # abs tol 1e-17 [(-11, 1.47824342), (-8, 8.9590946e-18), (-7, 1.85307619), (-4, 2.45138938)] sage: F = E.quadratic_twist(-8) sage: F.rank() 1 dict -- keys are the discriminants $d$, and values are list of corresponding zeros. EXAMPLES: EXAMPLES:: sage: E = EllipticCurve('37a') sage: E.lseries().twist_zeros(3, -4, -3)         # long time ***   Warning:...new stack size = ... from sage.lfunctions.lcalc import lcalc return lcalc.twist_zeros(n, dmin, dmax, L=self.__E) def at1(self, k=0): def at1(self, k=None, prec=None): r""" Compute $L(E,1)$ using $k$ terms of the series for $L(E,1)$ as explained on page 406 of Henri Cohen's book"A Course in Computational Algebraic Number Theory".  If the argument $k$ is not specified, then it defaults to $\sqrt(N)$, where $N$ is the conductor. The real precision used in each step of the computation is the precision of machine floats. Compute L(E,1) using k terms of the series for L(E,1) as explained on page 406 of Henri Cohen's book "A Course in Computational Algebraic Number Theory".  If the argument k is not specified, then it defaults to \sqrt(N), where N is the conductor. INPUT: k -- (optional) an integer, defaults to sqrt(N). - k -- number of terms of the series. If zero or None, use k = \sqrt(N), where N is the conductor. - prec -- numerical precision in bits. If zero or None, use a reasonable automatic default. OUTPUT: float -- L(E,1) float -- a bound on the error in the approximation; this is a provably correct upper bound on the sum of the tail end of the series used to compute L(E,1). This function is disjoint from the PARI \code{elllseries} A tuple of real numbers (L, err) where L is an approximation for L(E,1) and err is a bound on the error in the approximation. This function is disjoint from the PARI elllseries command, which is for a similar purpose.  To use that command (via the PARI C library), simply type \code{E.pari_mincurve().elllseries(1)} E.pari_mincurve().elllseries(1). ALGORITHM: \begin{enumerate} \item Compute the root number eps.  If it is -1, return 0. \item Compute the Fourier coefficients a_n, for n up to and including k. \item Compute the sum $$2 * sum_{n=1}^{k} (a_n / n) * exp(-2*pi*n/Sqrt(N)),$$ where N is the conductor of E. \item Compute a bound on the tail end of the series, which is $$2 * e^(-2 * pi * (k+1) / sqrt(N)) / (1 - e^(-2*pi/sqrt(N))).$$ For a proof see [Grigov-Jorza-Patrascu-Patrikis-Stein]. \end{enumerate} EXAMPLES: - Compute the root number eps.  If it is -1, return 0. - Compute the Fourier coefficients a_n, for n up to and including k. - Compute the sum .. MATH:: 2 * sum_{n=1}^{k} (a_n / n) * exp(-2*pi*n/Sqrt(N)), where N is the conductor of E. - Compute a bound on the tail end of the series, which is .. MATH:: 2 e^{-2 \pi (k+1) / \sqrt{N}} / (1 - e^{-2 \pi/\sqrt{N}}). For a proof see [Grigov-Jorza-Patrascu-Patrikis-Stein]. EXAMPLES:: sage: L, err = EllipticCurve('11a1').lseries().at1() sage: L, err (0.253804, 0.000181444) sage: parent(L) Real Field with 24 bits of precision sage: E = EllipticCurve('37b') sage: E.lseries().at1() (0.7257177, 0.000800697) sage: E.lseries().at1(100) (0.725681061936153, 1.52437502288743e-45) (0.7256810619361527823362055410263965487367603361763, 1.52469e-45) sage: L,err = E.lseries().at1(100, prec=128) sage: L 0.72568106193615278233620554102639654873 sage: parent(L) Real Field with 128 bits of precision sage: err 1.70693e-37 sage: parent(err) Real Field with 24 bits of precision and rounding RNDU Rank 1 through 3 elliptic curves:: sage: E = EllipticCurve('37a1') sage: E.lseries().at1() (0.0000000, 0.000000) sage: E = EllipticCurve('389a1') sage: E.lseries().at1() (-0.001769566, 0.00911776) sage: E = EllipticCurve('5077a1') sage: E.lseries().at1() (0.0000000, 0.000000) """ sqrtN = sqrt(self.__E.conductor()) if k: k = int(k) else: k = int(ceil(sqrtN)) if prec: prec = int(prec) else: # Use the same precision as deriv_at1() below for # consistency prec = int(9.065*k/sqrtN + 1.443*log(k)) + 12 R = RealField(prec) # Compute error term with bounded precision of 24 bits and # round towards +infinity Rerror = RealField(24, rnd='RNDU') if self.__E.root_number() == -1: return 0 sqrtN = float(self.__E.conductor().sqrt()) k = int(k) if k == 0: k = int(ceil(sqrtN)) an = self.__E.anlist(k)           # list of Sage ints # Compute z = e^(-2pi/sqrt(N)) pi = 3.14159265358979323846 z = exp(-2*pi/sqrtN) return (R.zero(), Rerror.zero()) an = self.__E.anlist(k)  # list of Sage Integers pi = R.pi() sqrtN = R(self.__E.conductor()).sqrt() z = (-2*pi/sqrtN).exp() zpow = z s = 0.0 # Compute series sum and accumulate floating point errors L = R.zero() error = Rerror.zero() for n in xrange(1,k+1): s += (zpow * float(an[n]))/n term = (zpow * an[n])/n zpow *= z L += term # We express relative error in units of epsilon, where # epsilon is a number divided by 2^precision. # Instead of multiplying the error by 2 after the loop # (to account for L *= 2), we already multiply it now. # # For multiplication and division, the relative error # in epsilons is bounded by (1+e)^n - 1, where n is the # number of operations (assuming exact inputs). # exp(x) additionally multiplies this error by abs(x) and # adds one epsilon. The inputs pi and sqrtN each contribute # another epsilon. # Assuming that 2*pi/sqrtN <= 2, the relative error for z is # 7 epsilon. This implies a relative error of (8n-1) epsilon # for zpow. We add 2 for the computation of term and 1/2 to # compensate for the approximation (1+e)^n = 1+ne. # # The error of the addition is at most half an ulp of the # result. # # Multiplying everything by two gives: error += term.epsilon(Rerror)*(16*n + 3) + L.ulp(Rerror) L *= 2 error = 2*zpow / (1 - z) return R(2*s), R(error) # Add series error (we use (-2)/(z-1) instead of 2/(1-z) # because this causes 1/(1-z) to be rounded up) error += ((-2)*Rerror(zpow)) / Rerror(z - 1) return (L, error) def deriv_at1(self, k=0): def deriv_at1(self, k=None, prec=None): r""" Compute $L'(E,1)$ using$k$ terms of the series for $L'(E,1)$. Compute L'(E,1) using k terms of the series for L'(E,1), under the assumption that L(E,1) = 0. The algorithm used is from page 406 of Henri Cohen's book A Course in Computational Algebraic Number Theory.'' The real precision of the computation is the precision of Python floats. INPUT: INPUT: k -- int; number of terms of the series - k -- number of terms of the series. If zero or None, use k = \sqrt(N), where N is the conductor. - prec -- numerical precision in bits. If zero or None, use a reasonable automatic default. OUTPUT: real number -- an approximation for L'(E,1) real number -- a bound on the error in the approximation A tuple of real numbers (L1, err) where L1 is an approximation for L'(E,1) and err is a bound on the error in the approximation. .. WARNING:: This function only makes sense if L(E) has positive order of vanishing at 1, or equivalently if L(E,1) = 0. ALGORITHM: \begin{enumerate} \item Compute the root number eps.  If it is 1, return 0. \item Compute the Fourier coefficients $a_n$, for $n$ up to and including $k$. \item Compute the sum $$- Compute the root number eps. If it is 1, return 0. - Compute the Fourier coefficients a_n, for n up to and including k. - Compute the sum .. MATH:: 2 * \sum_{n=1}^{k} (a_n / n) * E_1(2 \pi n/\sqrt{N}),$$ where $N$ is the conductor of $E$, and $E_1$ is the exponential integral function. \item Compute a bound on the tail end of the series, which is $$2 * e^{-2 \pi (k+1) / \sqrt{N}} / (1 - e^{-2 \ pi/\sqrt{N}}).$$ For a proof see [Grigorov-Jorza-Patrascu-Patrikis-Stein].  This is exactly the same as the bound for the approximation to $L(E,1)$ produced by \code{E.lseries().at1}. \end{enumerate} where N is the conductor of E, and E_1 is the exponential integral function. - Compute a bound on the tail end of the series, which is .. MATH:: 2 e^{-2 \pi (k+1) / \sqrt{N}} / (1 - e^{-2 \pi/\sqrt{N}}). For a proof see [Grigorov-Jorza-Patrascu-Patrikis-Stein].  This is exactly the same as the bound for the approximation to L(E,1) produced by :meth:at1. EXAMPLES:: sage: E = EllipticCurve('37a') sage: E.lseries().deriv_at1() (0.305986660898516, 0.000800351433106958) (0.3059866, 0.000801045) sage: E.lseries().deriv_at1(100) (0.305999773834052, 1.52437502288740e-45) (0.3059997738340523018204836833216764744526377745903, 1.52493e-45) sage: E.lseries().deriv_at1(1000) (0.305999773834052, 0.000000000000000) (0.305999773834052301820483683321676474452637774590771998..., 2.75031e-449) With less numerical precision, the error is bounded by numerical accuracy:: sage: L,err = E.lseries().deriv_at1(100, prec=64) sage: L,err (0.305999773834052302, 5.55318e-18) sage: parent(L) Real Field with 64 bits of precision sage: parent(err) Real Field with 24 bits of precision and rounding RNDU Rank 2 and rank 3 elliptic curves:: sage: E = EllipticCurve('389a1') sage: E.lseries().deriv_at1() (0.0000000, 0.000000) sage: E = EllipticCurve((1, 0, 1, -131, 558))  # curve 59450i1 sage: E.lseries().deriv_at1() (-0.00010911444, 0.142428) sage: E.lseries().deriv_at1(4000) (6.9902290...e-50, 1.31318e-43) """ if self.__E.root_number() == 1: return 0 k = int(k) sqrtN = float(self.__E.conductor().sqrt()) if k == 0: k = int(ceil(sqrtN)) an = self.__E.anlist(k)           # list of Sage Integers # Compute z = e^(-2pi/sqrt(N)) pi = 3.14159265358979323846 sqrtN = sqrt(self.__E.conductor()) if k: k = int(k) else: k = int(ceil(sqrtN)) if prec: prec = int(prec) else: # Estimate number of bits for the computation, based on error # estimate below (the denominator of that error is close enough # to 1 that we can ignore it). # 9.065 = 2*Pi/log(2) # 1.443 = 1/log(2) # 12 is an arbitrary extra number of bits (it is chosen #    such that the precision is 24 bits when the conductor #    equals 11 and k is the default value 4) prec = int(9.065*k/sqrtN + 1.443*log(k)) + 12 R = RealField(prec) # Compute error term with bounded precision of 24 bits and # round towards +infinity Rerror = RealField(24, rnd='RNDU') if self.__E.root_number() == 1: # Order of vanishing at 1 of L(E) is even and assumed to be # positive, so L'(E,1) = 0. return (R.zero(), Rerror.zero()) an = self.__E.anlist(k)  # list of Sage Integers pi = R.pi() sqrtN = R(self.__E.conductor()).sqrt() v = exp_integral.exponential_integral_1(2*pi/sqrtN, k) L = 2*float(sum([ (v[n-1] * an[n])/n for n in xrange(1,k+1)])) error = 2*exp(-2*pi*(k+1)/sqrtN)/(1-exp(-2*pi/sqrtN)) return R(L), R(error) # Compute series sum and accumulate floating point errors L = R.zero() error = Rerror.zero() # Sum of |an[n]|/n sumann = Rerror.zero() for n in xrange(1,k+1): term = (v[n-1] * an[n])/n L += term error += term.epsilon(Rerror)*5 + L.ulp(Rerror) sumann += Rerror(an[n].abs())/n L *= 2 # Add error term for exponential_integral_1() errors. # Absolute error for 2*v[i] is 4*max(1, v[0])*2^-prec if v[0] > 1.0: sumann *= Rerror(v[0]) error += (sumann >> (prec - 2)) # Add series error (we use (-2)/(z-1) instead of 2/(1-z) # because this causes 1/(1-z) to be rounded up) z = (-2*pi/sqrtN).exp() zpow = ((-2*(k+1))*pi/sqrtN).exp() error += ((-2)*Rerror(zpow)) / Rerror(z - 1) return (L, error) def __call__(self, s): r""" Returns the value of the L-series of the elliptic curve E at s, where s must be a real number. Use self.extended for s complex. .. NOTE:: \note{If the conductor of the curve is large, say $>10^{12}$, then this function will take a very long time, since it uses an $O(\sqrt{N})$ algorithm.} If the conductor of the curve is large, say >10^{12}, then this function will take a very long time, since it uses an O(\sqrt{N}) algorithm. EXAMPLES: EXAMPLES:: sage: E = EllipticCurve([1,2,3,4,5]) sage: L = E.lseries() sage: L(1) """ return self.dokchitser()(s) #def extended(self, s, prec): #    r""" #    Returns the value of the L-series of the elliptic curve E at s #    can be any complex number using prec terms of the power series #    expansion. # # #    WARNING: This may be slow.  Consider using \code{dokchitser()} #    instead. # #    INPUT: #        s -- complex number #        prec -- integer # #    EXAMPLES: #        sage: E = EllipticCurve('389a') #        sage: E.lseries().extended(1 + I, 50) #        -0.638409959099589 + 0.715495262192901*I #        sage: E.lseries().extended(1 + 0.1*I, 50) #        -0.00761216538818315 + 0.000434885704670107*I # #    NOTE: You might also want to use Tim Dokchitser's #    L-function calculator, which is available by typing #    L = E.lseries().dokchitser(), then evaluating L.  It #    gives the same information but is sometimes much faster. # #    """ #    try: #        s = C(s) #    except TypeError: #        raise TypeError, "Input argument %s must be coercible to a complex number"%s #    prec = int(prec) #    if abs(s.imag()) < R(0.0000000000001): #        return self(s.real()) #    N = self.__E.conductor() #    from sage.symbolic.constants import pi #    pi = R(pi) #    Gamma = transcendental.gamma #    Gamma_inc = transcendental.gamma_inc #    a = self.__E.anlist(prec) #    eps = self.__E.root_number() #    sqrtN = float(N.sqrt()) #    def F(n, t): #        return Gamma_inc(t+1, 2*pi*n/sqrtN) * C(sqrtN/(2*pi*n))**(t+1) #    return C(N)**(-s/2) * C(2*pi)**s * Gamma(s)**(-1)\ #           * sum([a[n]*(F(n,s-1) + eps*F(n,1-s)) for n in xrange(1,prec+1)]) def L1_vanishes(self): """ Returns whether or not L(E,1) = 0. The result is provably Returns whether or not L(E,1) = 0. The result is provably correct if the Manin constant of the associated optimal quotient is <= 2.  This hypothesis on the Manin constant is true for all curves of conductor <= 40000 (by Cremona) and all semistable curves (i.e., squarefree conductor). EXAMPLES: ALGORITHM: see :meth:L_ratio. EXAMPLES:: sage: E = EllipticCurve([0, -1, 1, -10, -20])   # 11A  = X_0(11) sage: E.lseries().L1_vanishes() False sage: E.lseries().L1_vanishes() False WARNING: It's conceivable that machine floats are not large enough precision for the computation; if this could be the case a RuntimeError is raised.  The curve's real period would have to be very small for this to occur. ALGORITHM: Compute the root number.  If it is -1 then L(E,s) vanishes to odd order at 1, hence vanishes.  If it is +1, use a result about modular symbols and Mazur's "Rational Isogenies" paper to determine a provably correct bound (assuming Manin constant is <= 2) so that we can determine whether L(E,1) = 0. AUTHOR: William Stein, 2005-04-20. """ return self.L_ratio() == 0 def L_ratio(self): r""" Returns the ratio $L(E,1)/\Omega$ as an exact rational Returns the ratio L(E,1)/\Omega as an exact rational number. The result is \emph{provably} correct if the Manin constant of the associated optimal quotient is $\leq 2$.  This constant of the associated optimal quotient is \leq 2.  This hypothesis on the Manin constant is true for all semistable curves (i.e., squarefree conductor), by a theorem of Mazur from his \emph{Rational Isogenies of Prime Degree} paper. EXAMPLES: EXAMPLES:: sage: E = EllipticCurve([0, -1, 1, -10, -20])   # 11A  = X_0(11) sage: E.lseries().L_ratio() 1/5 sage: E.lseries().L_ratio() 2 # See trac #3651: See :trac:3651 and :trac:15299:: sage: EllipticCurve([0,0,0,-193^2,0]).sha().an() 4 WARNING: It's conceivable that machine floats are not large enough precision for the computation; if this could be the case a RuntimeError is raised.  The curve's real period would have to be very small for this to occur. sage: EllipticCurve([1, 0, 1, -131, 558]).sha().an()  # long time 1.00000000000000 ALGORITHM: Compute the root number.  If it is -1 then L(E,s) vanishes to odd order at 1, hence vanishes.  If it is +1, use self.__lratio = self.__E.minimal_model().lseries().L_ratio() return self.__lratio QQ = RationalField() if self.__E.root_number() == -1: self.__lratio = Q(0) self.__lratio = QQ.zero() return self.__lratio # Even root number.  Decide if L(E,1) = 0.  If E is a modular d = self.__E._multiple_of_degree_of_isogeny_to_optimal_curve() C = 8*d*t eps = omega / C # coercion of 10**(-15) to our real field is needed to # make unambiguous comparison if eps < R(10**(-15)):  # liberal bound on precision of float raise RuntimeError, "Insufficient machine precision (=%s) for computation."%eps sqrtN = 2*int(sqrt(self.__E.conductor())) k = sqrtN + 10 while True: L1, error_bound = self.at1(k) if error_bound < eps: n = int(round(L1*C/omega)) quo = Q(n) / Q(C) quo = QQ((n,C)) self.__lratio = quo / self.__E.real_components() return self.__lratio 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 See :trac:1115:: sage: sha=EllipticCurve('37a1').sha() sage: [sha.an_numerical(prec) for prec in xrange(40,100,10)] sage: [sha.an_numerical(prec) for prec in xrange(40,100,10)]  # long time (3s on sage.math, 2013) [1.0000000000, 1.0000000000000, 1.0000000000000000,