source: sage/schemes/elliptic_curves/lseries_ell.py @ 6785:7d093366483b

Revision 6785:7d093366483b, 25.3 KB checked in by William Stein <wstein@…>, 6 years ago (diff)

Improve elliptic curve L-series and their doctests.
Line 
1"""
2Complex Elliptic Curve L-series
3
4"""
5
6from sage.structure.sage_object import SageObject
7from sage.rings.all import (
8    RealField,
9    RationalField,
10    ComplexField)
11from math import sqrt, exp
12import sage.functions.transcendental as transcendental
13R = RealField()
14Q = RationalField()
15C = ComplexField()
16import sage.functions.constants as constants
17
18class Lseries_ell(SageObject):
19    """
20    An elliptic curve $L$-series.
21
22    EXAMPLES:
23   
24    """
25    def __init__(self, E):
26        """
27        Create an elliptic curve $L$-series.
28
29        EXAMPLES:
30            sage: EllipticCurve([1..5]).Lseries()
31            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
32        """
33        self.__E = E
34
35    def elliptic_curve(self):
36        """
37        Return the elliptic curve that this L-series is attached to.
38       
39        EXAMPLES:
40            sage: E = EllipticCurve('389a')
41            sage: L = E.Lseries()
42            sage: L.elliptic_curve ()
43            Elliptic Curve defined by y^2 + y = x^3 + x^2 - 2*x over Rational Field
44        """
45        return self.__E
46
47    def taylor_series(self, a=1, prec=53, series_prec=6, var='z'):
48        """
49        Return the Taylor series of this $L$-series about $a$ to
50        the given precision (in bits) and the number of terms.
51
52        The output is a series in var, where you should view var as
53        equal to s-a.  Thus this function returns the formal power
54        series whose coefficients are L^{(n)}(a)/n!.
55       
56        EXAMPLES:
57            sage: E = EllipticCurve('389a')
58            sage: L = E.Lseries()
59            sage: L.taylor_series(series_prec=3)
60            -1.28158145691931e-23 + (7.26268290635587e-24)*z + 0.759316500288427*z^2 + O(z^3)
61        """
62        D = self.dokchitser(prec)
63        return D.taylor_series(a, series_prec, var)
64
65    def _repr_(self):
66        """
67        Return string representation of this L-series.
68       
69        EXAMPLES:
70            sage: E = EllipticCurve('37a')
71            sage: L = E.Lseries()
72            sage: L._repr_()
73            'Complex L-series of the Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field'
74        """
75        return "Complex L-series of the %s"%self.__E
76
77    def dokchitser(self, prec=53,
78                   max_imaginary_part=0,
79                   max_asymp_coeffs=40,
80                   algorithm='gp'):
81        r"""
82        Return interface to Tim Dokchitser's program for computing
83        with the L-series of this elliptic curve; this provides a way
84        to compute Taylor expansions and higher derivatives of
85        $L$-series.
86
87        INPUT:
88            prec -- integer (bits precision)
89            max_imaginary_part -- real number
90            max_asymp_coeffs -- integer
91            algorithm -- string: 'gp' or 'magma'
92
93        \note{If algorithm='magma', then the precision is in digits rather
94        than bits and the object returned is a Magma L-series, which has
95        different functionality from the SAGE L-series.}
96
97        EXAMPLES:
98            sage: E = EllipticCurve('37a')
99            sage: L = E.Lseries().dokchitser()
100            sage: L(2)
101            0.381575408260711
102            sage.: L = E.Lseries().dokchitser(algorithm='magma')         # optional  (not auto tested)
103            sage.: L.Evaluate(2)                                       # optional  (not auto tested)
104            0.38157540826071121129371040958008663667709753398892116
105        """
106        if algorithm == 'magma':
107            from sage.interfaces.all import magma
108            return magma(self.__E).LSeries(Precision = prec)
109
110        from sage.lfunctions.all import Dokchitser
111        key = (prec, max_imaginary_part, max_asymp_coeffs)
112        try:
113            return self.__dokchitser[key]
114        except KeyError:
115            pass
116        except AttributeError:
117            self.__dokchitser = {}
118        L = Dokchitser(conductor = self.__E.conductor(),
119                       gammaV = [0,1],
120                       weight = 2,
121                       eps = self.__E.root_number(),
122                       poles = [],
123                       prec = prec)
124        gp = L.gp()
125        s = 'e = ellinit(%s);'%self.__E.minimal_model().a_invariants()
126        s += 'a(k) = ellak(e, k);'
127        L.init_coeffs('a(k)', 1, pari_precode = s,
128                      max_imaginary_part=max_imaginary_part,
129                      max_asymp_coeffs=max_asymp_coeffs)
130        L.rename('Dokchitser L-function associated to %s'%self.__E)
131        self.__dokchitser[key] = L
132        return L
133
134    def sympow(self, n, prec):
135        r"""
136        Return $L(\Sym^{(n)}(E, \text{edge}))$ to prec digits
137        of precision.
138
139        INPUT:
140            n -- integer
141            prec -- integer
142           
143        OUTPUT:
144            string -- real number to prec digits of precision as a string.
145
146        \note{Before using this function for the first time for
147        a given $n$, you may have to type \code{sympow('-new_data <n>')},
148        where \code{<n>} is replaced by your value of $n$.  This
149        command takes a long time to run.}
150
151        EXAMPLES:
152            sage: E = EllipticCurve('37a')
153            sage: a = E.Lseries().sympow(2,16)   # optional -- requires precomputing "sympow('-new_data 2')"
154            sage: a      # optional
155            '2.492262044273650E+00'
156            sage: RR(a)                      # optional
157            2.4922620442736498
158        """
159        from sage.lfunctions.sympow import sympow
160        return sympow.L(self.__E, n, prec)
161
162    def sympow_derivs(self, n, prec, d):
163        r"""
164        Return $0$th to $d$th derivatives of $L(\Sym^{(n)}(E,
165        \text{edge}))$ to prec digits of precision.
166
167        INPUT:
168            n -- integer
169            prec -- integer
170            d -- integer
171           
172        OUTPUT:
173            a string, exactly as output by sympow
174
175        \note{To use this function you may have to run a few commands
176        like \code{sympow('-new_data 1d2')}, each which takes a few
177        minutes.  If this function fails it will indicate what
178        commands have to be run.}
179
180        EXAMPLES:
181            sage: E = EllipticCurve('37a')   
182            sage: E.Lseries().sympow_derivs(1,16,2)      # optional -- requires precomputing "sympow('-new_data 2')"
183            ...
184            1n0: 3.837774351482055E-01
185            1w0: 3.777214305638848E-01
186            1n1: 3.059997738340522E-01
187            1w1: 3.059997738340524E-01
188            1n2: 1.519054910249753E-01
189            1w2: 1.545605024269432E-01
190        """
191        from sage.lfunctions.sympow import sympow
192        return sympow.Lderivs(self.__E, n, prec, d)
193
194    def zeros(self, n):
195        """
196        Return the imaginary parts of the first $n$ nontrivial zeros
197        on the critical line of the L-function in the upper half
198        plane, as 32-bit reals.
199       
200        EXAMPLES:
201            sage: E = EllipticCurve('37a')
202            sage: E.Lseries().zeros(2)                 
203            [0.000000000, 5.00317001]
204
205            sage: a = E.Lseries().zeros(20)             # long time
206            sage: point([(1,x) for x in a]).save()    # graph  (long time)
207
208        AUTHOR:
209            -- Uses Rubinstein's L-functions calculator.
210        """
211        from sage.lfunctions.lcalc import lcalc
212        return lcalc.zeros(n, L=self.__E)
213
214    def zeros_in_interval(self, x, y, stepsize):
215        r"""
216        Return the imaginary parts of (most of) the nontrivial zeros
217        on the critical line $\Re(s)=1$ with positive imaginary part
218        between $x$ and $y$, along with a technical quantity for each.
219
220        INPUT:
221            x, y, stepsize -- positive floating point numbers
222
223        OUTPUT:
224            list of pairs (zero, S(T)).
225
226        Rubinstein writes: The first column outputs the imaginary part
227        of the zero, the second column a quantity related to S(T) (it
228        increases roughly by 2 whenever a sign change, i.e. pair of
229        zeros, is missed). Higher up the critical strip you should use
230        a smaller stepsize so as not to miss zeros.
231
232        EXAMPLES:
233            sage: E = EllipticCurve('37a')
234            sage: E.Lseries().zeros_in_interval(6, 10, 0.1)      # long
235            [(6.87039122, 0.248922780), (8.01433081, -0.140168533), (9.93309835, -0.129943029)]
236        """
237        from sage.lfunctions.lcalc import lcalc
238        return lcalc.zeros_in_interval(x, y, stepsize, L=self.__E)
239
240    def values_along_line(self, s0, s1, number_samples):
241        """
242        Return values of $L(E, s)$ at \code{number_samples}
243        equally-spaced sample points along the line from $s_0$ to
244        $s_1$ in the complex plane.
245
246        \note{The L-series is normalized so that the center of the
247        critical strip is 1.}
248
249        INPUT:
250            s0, s1 -- complex numbers
251            number_samples -- integer
252           
253        OUTPUT:
254            list -- list of pairs (s, zeta(s)), where the s are
255                    equally spaced sampled points on the line from
256                    s0 to s1.
257
258        EXAMPLES:
259            sage: I = CC.0
260            sage: E = EllipticCurve('37a')
261            sage: E.Lseries().values_along_line(1, 0.5+20*I, 5)     # long time and slightly random output
262            [(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)]
263        """
264        from sage.lfunctions.lcalc import lcalc
265        return lcalc.values_along_line(s0-RationalField()('1/2'),
266                                       s1-RationalField()('1/2'),
267                                       number_samples, L=self.__E)
268
269    def twist_values(self, s, dmin, dmax):
270        r"""
271        Return values of $L(E, s, \chi_d)$ for each quadratic
272        character $\chi_d$ for $d_{\min} \leq d \leq d_{\max}$.
273
274        \note{The L-series is normalized so that the center of the
275        critical strip is 1.}
276
277        INPUT:
278            s -- complex numbers
279            dmin -- integer
280            dmax -- integer
281
282        OUTPUT:
283            list -- list of pairs (d, L(E, s,chi_d))
284
285        EXAMPLES:
286            sage: E = EllipticCurve('37a')
287            sage: E.Lseries().twist_values(1, -12, -4)    # slightly random output depending on architecture
288            [(-11, 1.4782434171), (-8, 0), (-7, 1.8530761916), (-4, 2.4513893817)]
289            sage: F = E.quadratic_twist(-8)
290            sage: F.rank()
291            1
292            sage: F = E.quadratic_twist(-7)
293            sage: F.rank()
294            0
295        """
296        from sage.lfunctions.lcalc import lcalc
297        return lcalc.twist_values(s - RationalField()('1/2'), dmin, dmax, L=self.__E)
298
299    def twist_zeros(self, n, dmin, dmax):
300        r"""
301        Return first $n$ real parts of nontrivial zeros of
302        $L(E,s,\chi_d)$ for each quadratic character $\chi_d$ with
303        $d_{\min} \leq d \leq d_{\max}$.
304
305        \note{The L-series is normalized so that the center of the
306        critical strip is 1.}
307
308        INPUT:
309            n -- integer
310            dmin -- integer
311            dmax -- integer
312
313        OUTPUT:
314            dict -- keys are the discriminants $d$, and
315                    values are list of corresponding zeros.
316
317        EXAMPLES:
318            sage: E = EllipticCurve('37a')
319            sage: E.Lseries().twist_zeros(3, -4, -3)         # long
320            {-4: [1.60813783, 2.96144840, 3.89751747], -3: [2.06170900, 3.48216881, 4.45853219]}
321        """
322        from sage.lfunctions.lcalc import lcalc
323        return lcalc.twist_zeros(n, dmin, dmax, L=self.__E)
324
325    def at1(self, k=0):
326        r"""
327        Compute $L(E,1)$ using $k$ terms of the series for $L(E,1)$ as
328        explained on page 406 of Henri Cohen's book"A Course in Computational
329        Algebraic Number Theory".  If the argument $k$ is not specified,
330        then it defaults to $\sqrt(N)$, where $N$ is the conductor.
331
332        The real precision used in each step of the computation is the
333        precision of machine floats.
334
335        INPUT:
336            k -- (optional) an integer, defaults to sqrt(N).
337           
338        OUTPUT:
339            float -- L(E,1)
340            float -- a bound on the error in the approximation; this
341                     is a proveably correct upper bound on the sum
342                     of the tail end of the series used to compute L(E,1).
343
344        This function is disjoint from the PARI \code{elllseries}
345        command, which is for a similar purpose.  To use that command
346        (via the PARI C library), simply type
347                \code{E.pari_mincurve().elllseries(1)}
348
349        ALGORITHM:
350        \begin{enumerate}
351            \item Compute the root number eps.  If it is -1, return 0.
352           
353            \item Compute the Fourier coefficients a_n, for n up to and
354               including k.
355               
356            \item Compute the sum
357            $$
358                 2 * sum_{n=1}^{k} (a_n / n) * exp(-2*pi*n/Sqrt(N)),
359            $$
360               where N is the conductor of E.
361               
362            \item Compute a bound on the tail end of the series, which is
363            $$           
364                 2 * e^(-2 * pi * (k+1) / sqrt(N)) / (1 - e^(-2*pi/sqrt(N))).
365            $$     
366               For a proof see [Grigov-Jorza-Patrascu-Patrikis-Stein].
367        \end{enumerate}
368
369        EXAMPLES:
370            sage: E = EllipticCurve('37b')
371            sage: E.Lseries().at1(100)
372            (0.725681061936153, 1.52437502288743e-45)
373        """
374        if self.__E.root_number() == -1:
375            return 0
376        sqrtN = float(self.__E.conductor().sqrt())
377        k = int(k)
378        if k == 0: k = int(math.ceil(sqrtN))
379        an = self.__E.anlist(k)           # list of SAGE ints
380        # Compute z = e^(-2pi/sqrt(N))
381        pi = 3.14159265358979323846
382        z = exp(-2*pi/sqrtN)
383        zpow = z
384        s = 0.0
385        for n in xrange(1,k+1):
386            s += (zpow * float(an[n]))/n
387            zpow *= z
388
389        error = 2*zpow / (1 - z)
390       
391        return R(2*s), R(error)
392
393    def deriv_at1(self, k=0):
394        r"""
395        Compute $L'(E,1)$ using$ k$ terms of the series for $L'(E,1)$.
396
397        The algorithm used is from page 406 of Henri Cohen's book ``A
398        Course in Computational Algebraic Number Theory.''
399
400        The real precision of the computation is the precision of
401        Python floats.
402
403        INPUT:
404            k -- int; number of terms of the series
405
406        OUTPUT:
407            real number -- an approximation for L'(E,1)
408            real number -- a bound on the error in the approximation
409
410        ALGORITHM:
411        \begin{enumerate}
412            \item Compute the root number eps.  If it is 1, return 0.
413
414            \item Compute the Fourier coefficients $a_n$, for $n$ up to and
415                  including $k$.
416               
417            \item Compute the sum
418               $$
419                 2 * \sum_{n=1}^{k} (a_n / n) * E_1(2 \pi n/\sqrt{N}),
420               $$ 
421               where $N$ is the conductor of $E$, and $E_1$ is the
422               exponential integral function.
423
424            \item Compute a bound on the tail end of the series, which is
425               $$
426                 2 * e^{-2 \pi (k+1) / \sqrt{N}} / (1 - e^{-2 \ pi/\sqrt{N}}).
427               $$ 
428               For a proof see [Grigorov-Jorza-Patrascu-Patrikis-Stein].  This
429               is exactly the same as the bound for the approximation to
430               $L(E,1)$ produced by \code{E.Lseries().at1}.
431        \end{enumerate}
432
433        EXAMPLES:
434            sage: E = EllipticCurve('37a')
435            sage: E.Lseries().deriv_at1(100)
436            (0.305999773834879, 1.52437502288740e-45)
437        """
438        if self.__E.root_number() == 1: return 0
439        k = int(k)
440        sqrtN = float(self.__E.conductor().sqrt())
441        if k == 0: k = int(math.ceil(sqrtN))
442        an = self.__E.anlist(k)           # list of C ints
443        # Compute z = e^(-2pi/sqrt(N))
444        pi = 3.14159265358979323846
445        v = transcendental.exponential_integral_1(2*pi/sqrtN, k)
446        L = 2*float(sum([ (v[n-1] * an[n])/n for n in xrange(1,k+1)]))
447        error = 2*exp(-2*pi*(k+1)/sqrtN)/(1-exp(-2*pi/sqrtN))
448        return R(L), R(error)
449
450    def __call__(self, s):
451        r"""
452        Returns the value of the L-series of the elliptic curve E at s, where s
453        must be a real number.
454
455        Use self.extended for s complex.
456
457        \note{If the conductor of the curve is large, say $>10^{12}$,
458        then this function will take a very long time, since it uses
459        an $O(\sqrt{N})$ algorithm.}
460
461        EXAMPLES:
462            sage: E = EllipticCurve([1,2,3,4,5])
463            sage: L = E.Lseries()
464            sage: L(1)
465            0
466            sage: L(1.1)
467            0.285491007678148
468            sage: L(1.1 + I)
469            0.174851377216615 + 0.816965038124457*I
470        """
471        return self.dokchitser()(s)
472
473    #def extended(self, s, prec):
474    #    r"""
475    #    Returns the value of the L-series of the elliptic curve E at s
476    #    can be any complex number using prec terms of the power series
477    #    expansion.
478    #
479    #
480    #    WARNING: This may be slow.  Consider using \code{dokchitser()}
481    #    instead.
482    #   
483    #    INPUT:
484    #        s -- complex number
485    #        prec -- integer
486    #
487    #    EXAMPLES:
488    #        sage: E = EllipticCurve('389a')
489    #        sage: E.Lseries().extended(1 + I, 50)
490    #        -0.638409959099589 + 0.715495262192901*I
491    #        sage: E.Lseries().extended(1 + 0.1*I, 50)
492    #        -0.00761216538818315 + 0.000434885704670107*I
493    #
494    #    NOTE: You might also want to use Tim Dokchitser's
495    #    L-function calculator, which is available by typing
496    #    L = E.Lseries().dokchitser(), then evaluating L.  It
497    #    gives the same information but is sometimes much faster.
498    #   
499    #    """
500    #    try:
501    #        s = C(s)
502    #    except TypeError:
503    #        raise TypeError, "Input argument %s must be coercible to a complex number"%s
504    #    prec = int(prec)
505    #    if abs(s.imag()) < R(0.0000000000001):
506    #        return self(s.real())
507    #    N = self.__E.conductor()
508    #    pi = R(constants.pi)
509    #    Gamma = transcendental.gamma
510    #    Gamma_inc = transcendental.gamma_inc
511    #    a = self.__E.anlist(prec)
512    #    eps = self.__E.root_number()
513    #    sqrtN = float(N.sqrt())
514    #    def F(n, t):
515    #        return Gamma_inc(t+1, 2*pi*n/sqrtN) * C(sqrtN/(2*pi*n))**(t+1)
516    #    return C(N)**(-s/2) * C(2*pi)**s * Gamma(s)**(-1)\
517    #           * sum([a[n]*(F(n,s-1) + eps*F(n,1-s)) for n in xrange(1,prec+1)])
518
519    def L1_vanishes(self):
520        """
521        Returns whether or not L(E,1) = 0. The result is provably
522        correct if the Manin constant of the associated optimal
523        quotient is <= 2.  This hypothesis on the Manin constant
524        is true for all curves of conductor <= 40000 (by Cremona) and
525        all semistable curves (i.e., squarefree conductor).
526
527        EXAMPLES:
528            sage: E = EllipticCurve([0, -1, 1, -10, -20])   # 11A  = X_0(11)
529            sage: E.Lseries().L1_vanishes()
530            False
531            sage: E = EllipticCurve([0, -1, 1, 0, 0])       # X_1(11)
532            sage: E.Lseries().L1_vanishes()
533            False
534            sage: E = EllipticCurve([0, 0, 1, -1, 0])       # 37A  (rank 1)
535            sage: E.Lseries().L1_vanishes()
536            True
537            sage: E = EllipticCurve([0, 1, 1, -2, 0])       # 389A (rank 2)
538            sage: E.Lseries().L1_vanishes()
539            True
540            sage: E = EllipticCurve([0, 0, 1, -38, 90])     # 361A (CM curve))
541            sage: E.Lseries().L1_vanishes()
542            True
543            sage: E = EllipticCurve([0,-1,1,-2,-1])         # 141C (13-isogeny)
544            sage: E.Lseries().L1_vanishes()
545            False
546
547        WARNING: It's conceivable that machine floats are not large
548        enough precision for the computation; if this could be the
549        case a RuntimeError is raised.  The curve's real period would
550        have to be very small for this to occur. 
551
552        ALGORITHM: Compute the root number.  If it is -1 then L(E,s)
553        vanishes to odd order at 1, hence vanishes.  If it is +1, use
554        a result about modular symbols and Mazur's "Rational Isogenies"
555        paper to determine a provably correct bound (assuming Manin
556        constant is <= 2) so that we can determine whether L(E,1) = 0.
557
558        AUTHOR: William Stein, 2005-04-20.
559        """
560        return self.L_ratio() == 0
561
562    def L_ratio(self):
563        r"""
564        Returns the ratio $L(E,1)/\Omega$ as an exact rational
565        number. The result is \emph{provably} correct if the Manin
566        constant of the associated optimal quotient is $\leq 2$.  This
567        hypothesis on the Manin constant is true for all semistable
568        curves (i.e., squarefree conductor), by a theorem of Mazur
569        from his \emph{Rational Isogenies of Prime Degree} paper.
570
571        EXAMPLES:
572            sage: E = EllipticCurve([0, -1, 1, -10, -20])   # 11A  = X_0(11)
573            sage: E.Lseries().L_ratio()
574            1/5
575            sage: E = EllipticCurve([0, -1, 1, 0, 0])       # X_1(11)
576            sage: E.Lseries().L_ratio()
577            1/25
578            sage: E = EllipticCurve([0, 0, 1, -1, 0])       # 37A  (rank 1)
579            sage: E.Lseries().L_ratio()
580            0
581            sage: E = EllipticCurve([0, 1, 1, -2, 0])       # 389A (rank 2)
582            sage: E.Lseries().L_ratio()
583            0
584            sage: E = EllipticCurve([0, 0, 1, -38, 90])     # 361A (CM curve))
585            sage: E.Lseries().L_ratio()
586            0
587            sage: E = EllipticCurve([0,-1,1,-2,-1])         # 141C (13-isogeny)
588            sage: E.Lseries().L_ratio()
589            1
590            sage: E = EllipticCurve(RationalField(), [1, 0, 0, 1/24624, 1/886464])
591            sage: E.Lseries().L_ratio()
592            2
593
594        WARNING: It's conceivable that machine floats are not large
595        enough precision for the computation; if this could be the
596        case a RuntimeError is raised.  The curve's real period would
597        have to be very small for this to occur.
598
599        ALGORITHM: Compute the root number.  If it is -1 then L(E,s)
600        vanishes to odd order at 1, hence vanishes.  If it is +1, use
601        a result about modular symbols and Mazur's "Rational Isogenies"
602        paper to determine a provably correct bound (assuming Manin
603        constant is <= 2) so that we can determine whether L(E,1) = 0.
604
605        AUTHOR: William Stein, 2005-04-20.
606        """
607        try:
608            return self.__lratio
609        except AttributeError:
610            pass
611       
612        if not self.__E.is_minimal():
613            self.__lratio = self.__E.minimal_model().Lseries().L_ratio()
614            return self.__lratio
615
616        if self.__E.root_number() == -1:
617            self.__lratio = Q(0)
618            return self.__lratio
619
620        # Even root number.  Decide if L(E,1) = 0.  If E is a modular
621        # *OPTIMAL* quotient of J_0(N) elliptic curve, we know that T *
622        # L(E,1)/omega is an integer n, where T is the order of the
623        # image of the rational torsion point (0)-(oo) in E(Q), and
624        # omega is the least real Neron period.  (This is proved in my
625        # Ph.D. thesis, but is probably well known.)  We can easily
626        # compute omega to very high precision using AGM.  So to prove
627        # that L(E,1) = 0 we compute T/omega * L(E,1) to sufficient
628        # precision to determine it as an integer.  If eps is the
629        # error in computation of L(E,1), then the error in computing
630        # the product is (2T/Omega_E) * eps, and we need this to be
631        # less than 0.5, i.e.,
632        #          (2T/Omega_E) * eps < 0.5,
633        # so
634        #          eps < 0.5 * Omega_E / (2T) = Omega_E / (4*T).
635        #
636        # Since in general E need not be optimal, we have to choose
637        # eps = Omega_E/(8*t*B), where t is the exponent of E(Q)_tor,
638        # and is a multiple of the degree of an isogeny between E
639        # and the optimal curve.
640        #
641        # NOTES: We *do* have to worry about the Manin constant, since
642        # we are using the Neron model to compute omega, not the
643        # newform.  My theorem replaces the omega above by omega/c,
644        # where c is the Manin constant, and the bound must be
645        # correspondingly smaller.  If the level is square free, then
646        # the Manin constant is 1 or 2, so there's no problem (since
647        # we took 8 instead of 4 in the denominator).  If the level
648        # is divisible by a square, then the Manin constant could
649        # be a divisible by an arbitrary power of that prime, except
650        # that Edixhoven claims the primes that appear are <= 7.
651       
652        t = self.__E.torsion_subgroup().order()
653        omega = self.__E.period_lattice().basis()[0]
654        d = self.__E._multiple_of_degree_of_isogeny_to_optimal_curve()
655        C = 8*d*t
656        eps = omega / C
657        # coercion of 10**(-15) to our real field is needed to
658        # make unambiguous comparison
659        if eps < R(10**(-15)):  # liberal bound on precision of float
660            raise RuntimeError, "Insufficient machine precision (=%s) for computation."%eps
661        sqrtN = 2*int(sqrt(self.__E.conductor()))
662        k = sqrtN + 10
663        while True:
664            L1, error_bound = self.at1(k)
665            if error_bound < eps:
666                n = int(round(L1*C/omega))
667                quo = Q(n) / Q(C)
668                self.__lratio = quo / self.__E.real_components()
669                return self.__lratio
670            k += sqrtN
671            misc.verbose("Increasing precision to %s terms."%k)
Note: See TracBrowser for help on using the repository browser.