Ticket #3674: ell_rational_field.py

File ell_rational_field.py, 160.8 KB (added by cremona, 9 years ago)

Not a patch (yet)

Line 
1"""
2Elliptic curves over the rational numbers
3
4AUTHORS:
5   -- William Stein (2005): first version
6   -- William Stein (2006-02-26): fixed Lseries_extended which didn't work
7            because of changes elsewhere in SAGE.
8   -- David Harvey (2006-09): Added padic_E2, padic_sigma, padic_height,
9            padic_regulator methods.
10   -- David Harvey (2007-02): reworked padic-height related code
11   -- Christian Wuthrich (2007): added padic sha computation
12   -- David Roe (2007-9): moved sha, l-series and p-adic functionality to separate files.
13   -- John Cremona (2008-01)
14   -- Tobias Nagel & Michael Mardaus (2008-07): added integral_points
15"""
16
17#*****************************************************************************
18#       Copyright (C) 2005,2006,2007 William Stein <wstein@gmail.com>
19#
20#  Distributed under the terms of the GNU General Public License (GPL)
21#
22#    This code is distributed in the hope that it will be useful,
23#    but WITHOUT ANY WARRANTY; without even the implied warranty of
24#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
25#    General Public License for more details.
26#
27#  The full text of the GPL is available at:
28#
29#                  http://www.gnu.org/licenses/
30#*****************************************************************************
31
32import ell_point
33import formal_group
34import rational_torsion
35from ell_generic import EllipticCurve_generic
36from ell_number_field import EllipticCurve_number_field
37
38import sage.groups.all
39import sage.rings.arith as arith
40import sage.rings.all as rings
41import sage.rings.number_field.number_field as number_field
42import sage.misc.misc as misc
43from sage.misc.all import verbose
44import sage.functions.constants as constants
45import sage.modular.modform.constructor
46import sage.modular.modform.element
47from sage.misc.functional import log
48from sage.rings.padics.factory import Zp, Qp
49
50# Use some interval arithmetic to guarantee correctness.  We assume
51# that alpha is computed to the precision of a float.
52# IR = rings.RIF
53#from sage.rings.interval import IntervalRing; IR = IntervalRing()
54
55import sage.matrix.all as matrix
56import sage.databases.cremona
57from   sage.libs.pari.all import pari
58import sage.functions.transcendental as transcendental
59import math
60from sage.calculus.calculus import sqrt, arcsin, floor, ceil
61import sage.libs.mwrank.all as mwrank
62import constructor
63from sage.interfaces.all import gp
64
65import ell_modular_symbols
66import padic_lseries
67import padics
68
69from lseries_ell import Lseries_ell
70
71import mod5family
72
73from sage.rings.all import (
74    PowerSeriesRing, LaurentSeriesRing, O, 
75    infinity as oo,
76    Integer,
77    Integers,
78    IntegerRing, RealField,
79    ComplexField, RationalField)
80
81import gp_cremona
82import sea
83
84from gp_simon import simon_two_descent
85
86import ell_tate_curve
87
88factor = arith.factor
89exp = math.exp
90mul = misc.mul
91next_prime = arith.next_prime
92kronecker_symbol = arith.kronecker_symbol
93
94Q = RationalField()         
95C = ComplexField()
96R = RealField()
97Z = IntegerRing()
98IR = rings.RealIntervalField(20)
99
100_MAX_HEIGHT=21
101
102# CMJ is a dict of pairs (j,D) where j is a rational CM j-invariant
103# and D is the corresponding quadratic discriminant
104
105CMJ={ 0: -3, 54000: -12, -12288000: -27, 1728: -4, 287496: -16, 
106      -3375: -7, 16581375: -28, 8000: -8, -32768: -11, -884736: -19, 
107      -884736000: -43, -147197952000: -67, -262537412640768000: -163}
108
109
110
111class EllipticCurve_rational_field(EllipticCurve_number_field):
112    """
113    Elliptic curve over the Rational Field.
114    """
115    def __init__(self, ainvs, extra=None):
116        if extra != None:   # possibility of two arguments (the first would be the field)
117            ainvs = extra
118        if isinstance(ainvs, str):
119            label = ainvs
120            X = sage.databases.cremona.CremonaDatabase()[label]
121            EllipticCurve_number_field.__init__(self, Q, X.a_invariants())
122            for attr in ['rank', 'torsion_order', 'cremona_label', 'conductor',
123                         'modular_degree', 'gens', 'regulator']:
124                s = "_EllipticCurve_rational_field__"+attr
125                if hasattr(X,s):
126                    setattr(self, s, getattr(X, s))
127            return
128        EllipticCurve_number_field.__init__(self, Q, ainvs)
129        self.__np = {}
130        self.__gens = {}
131        self.__rank = {}
132        self.__regulator = {}
133        if self.base_ring() != Q:
134            raise TypeError, "Base field (=%s) must be the Rational Field."%self.base_ring()
135       
136    def _set_rank(self, r):
137        """
138        Internal function to set the cached rank of this elliptic curve to r.
139
140        WARNING:  No checking is done!  Not intended for use by users.
141
142        EXAMPLES:
143            sage: E=EllipticCurve('37a1')
144            sage: E._set_rank(99)  # bogus value -- not checked
145            sage: E.rank()         # returns bogus cached value
146            99
147            sage: E.gens()         # causes actual rank to be computed
148            [(0 : -1 : 1)]
149            sage: E.rank()         # the correct rank
150            1
151        """
152        self.__rank = {}
153        self.__rank[True] = Integer(r)
154
155    def _set_torsion_order(self, t):
156        """
157        Internal function to set the cached torsion order of this elliptic curve to t.
158
159        WARNING:  No checking is done!  Not intended for use by users.
160
161        EXAMPLES:
162            sage: E=EllipticCurve('37a1')
163            sage: E._set_torsion_order(99)  # bogus value -- not checked
164            sage: E.torsion_order()         # returns bogus cached value
165            99
166            sage: T = E.torsion_subgroup()  # causes actual torsion to be computed
167            sage: E.torsion_order()         # the correct value
168            1
169        """
170        self.__torsion_order = Integer(t)
171
172    def _set_cremona_label(self, L):
173        """
174        Internal function to set the cached label of this elliptic curve to L.
175
176        WARNING:  No checking is done!  Not intended for use by users.
177
178        EXAMPLES:
179            sage: E=EllipticCurve('37a1')
180            sage: E._set_cremona_label('bogus')
181            sage: E.label()
182            'bogus'
183            sage: E.database_curve().label()
184            '37a1'
185            sage: E.label() # no change
186            'bogus'
187            sage: E._set_cremona_label(E.database_curve().label())
188            sage: E.label() # now it is correct
189            '37a1'
190        """
191        self.__cremona_label = L
192
193    def _set_conductor(self, N):
194        """
195        Internal function to set the cached conductor of this elliptic curve to N.
196
197        WARNING: No checking is done!  Not intended for use by users.
198        Setting to the wrong value will cause strange problems (see
199        examples).
200       
201        EXAMPLES:
202            sage: E=EllipticCurve('37a1')
203            sage: E._set_conductor(99)      # bogus value -- not checked
204            sage: E.conductor()             # returns bogus cached value
205            99
206
207        This will not work since the conductor is used when searching
208        the database:
209            sage: E._set_conductor(E.database_curve().conductor())
210            Traceback (most recent call last):
211            ... 
212            RuntimeError: Elliptic curve ... not in the database.
213            sage: E._set_conductor(EllipticCurve(E.a_invariants()).database_curve().conductor())
214            sage: E.conductor()             # returns correct value
215            37
216        """
217        self.__conductor_pari = Integer(N)
218
219    def _set_modular_degree(self, deg):
220        """
221        Internal function to set the cached modular degree of this elliptic curve to deg.
222
223        WARNING:  No checking is done!
224
225        EXAMPLES:
226            sage: E=EllipticCurve('5077a1')
227            sage: E.modular_degree()
228            1984
229            sage: E._set_modular_degree(123456789)
230            sage: E.modular_degree()
231            123456789
232            sage: E._set_modular_degree(E.database_curve().modular_degree())
233            sage: E.modular_degree()
234            1984
235        """
236        self.__modular_degree = Integer(deg)
237       
238    def _set_gens(self, gens):
239        """
240        Internal function to set the cached generators of this elliptic curve to gens.
241
242        WARNING:  No checking is done!
243
244        EXAMPLES:
245            sage: E=EllipticCurve('5077a1')
246            sage: E.rank()
247            3
248            sage: E.gens() # random
249            [(-2 : 3 : 1), (-7/4 : 25/8 : 1), (1 : -1 : 1)]
250            sage: E._set_gens([]) # bogus list
251            sage: E.rank()        # unchanged
252            3
253            sage: E._set_gens(E.database_curve().gens())
254            sage: E.gens()
255            [(-2 : 3 : 1), (-7/4 : 25/8 : 1), (1 : -1 : 1)]
256        """
257        self.__gens = {}
258        self.__gens[True] = [self.point(x, check=True) for x in gens]
259        self.__gens[True].sort()
260
261    def is_integral(self):
262        """
263        Returns True if this elliptic curve has integral coefficients (in Z)
264
265        EXAMPLES:
266            sage: E=EllipticCurve(QQ,[1,1]); E
267            Elliptic Curve defined by y^2  = x^3 + x +1 over Rational Field
268            sage: E.is_integral()
269            True
270            sage: E2=E.change_weierstrass_model(2,0,0,0); E2
271            Elliptic Curve defined by y^2  = x^3 + 1/16*x + 1/64 over Rational Field
272            sage: E2.is_integral()
273            False
274         """
275        try:
276            return self.__is_integral
277        except AttributeError:
278            one = Integer(1)
279            self.__is_integral = bool(misc.mul([x.denominator() == 1 for x in self.ainvs()]))
280            return self.__is_integral
281           
282
283    def mwrank(self, options=''):
284        """
285        Run Cremona's mwrank program on this elliptic curve and
286        return the result as a string.
287
288        INPUT:
289            options -- string; passed when starting mwrank.  The format is
290        q p<precision> v<verbosity> b<hlim_q> x<naux>  c<hlim_c> l t o s d>]
291
292        OUTPUT:
293            string -- output of mwrank on this curve
294
295        NOTE: The output is a raw string and completely illegible
296            using automatic display, so it is recommended to use print
297            for legible outout
298
299        EXAMPLES:
300            sage: E = EllipticCurve('37a1')
301            sage: E.mwrank() #random
302            ...
303            sage: print E.mwrank()
304            Curve [0,0,1,-1,0] :        Basic pair: I=48, J=-432
305            disc=255744
306            ...
307            Generator 1 is [0:-1:1]; height 0.05111...
308            <BLANKLINE>
309            Regulator = 0.05111...
310            <BLANKLINE>
311            The rank and full Mordell-Weil basis have been determined unconditionally.
312            ...
313
314        Options to mwrank can be passed:
315            sage: E = EllipticCurve([0,0,0,877,0])
316
317        Run mwrank with 'verbose' flag set to 0 but list generators if found
318            sage: print E.mwrank('-v0 -l')
319            Curve [0,0,0,877,0] :   0 <= rank <= 1
320            Regulator = 1
321           
322        Run mwrank again, this time with a higher bound for point
323        searching on homogeneous spaces:
324
325            sage: print E.mwrank('-v0 -l -b11')
326            Curve [0,0,0,877,0] :   Rank = 1
327            Generator 1 is [29604565304828237474403861024284371796799791624792913256602210:-256256267988926809388776834045513089648669153204356603464786949:490078023219787588959802933995928925096061616470779979261000]; height 95.980371987964
328            Regulator = 95.980371987964
329        """
330        if options == "":
331            from sage.interfaces.all import mwrank
332        else:
333            from sage.interfaces.all import Mwrank
334            mwrank = Mwrank(options=options)
335        return mwrank(self.a_invariants())
336
337    def conductor(self, algorithm="pari"):
338        """
339        Returns the conductor of the elliptic curve.
340
341        INPUT:
342            algorithm -- str, (default: "pari")
343                   "pari"   -- use the PARI C-library ellglobalred
344                               implementation of Tate's algorithm
345                   "mwrank" -- use Cremona's mwrank implementation of
346                               Tate's algorithm; can be faster if the
347                               curve has integer coefficients (TODO:
348                               limited to small conductor until mwrank
349                               gets integer factorization)
350                   "gp" -- use the GP interpreter.
351                   "all" -- use both implementations, verify that the
352                            results are the same (or raise an error),
353                            and output the common value.
354                                     
355        EXAMPLE:
356            sage: E = EllipticCurve([1, -1, 1, -29372, -1932937])
357            sage: E.conductor(algorithm="pari")
358            3006
359            sage: E.conductor(algorithm="mwrank")
360            3006
361            sage: E.conductor(algorithm="gp")
362            3006
363            sage: E.conductor(algorithm="all")
364            3006
365
366        NOTE: The conductor computed using each algorithm is cached separately.
367        Thus calling E.conductor("pari"), then E.conductor("mwrank") and
368        getting the same result checks that both systems compute the same answer.
369        """
370
371        if algorithm == "pari":
372            try:
373                return self.__conductor_pari
374            except AttributeError:
375                self.__conductor_pari = Integer(self.pari_mincurve().ellglobalred()[0])
376            return self.__conductor_pari
377
378        elif algorithm == "gp":
379            try:
380                return self.__conductor_gp
381            except AttributeError:
382                self.__conductor_gp = Integer(gp.eval('ellglobalred(ellinit(%s,0))[1]'%self.a_invariants()))
383                return self.__conductor_gp
384
385        elif algorithm == "mwrank":
386            try:
387                return self.__conductor_mwrank
388            except AttributeError:
389                if self.is_integral():
390                    self.__conductor_mwrank = Integer(self.mwrank_curve().conductor())
391                else:
392                    self.__conductor_mwrank = Integer(self.minimal_model().mwrank_curve().conductor())
393            return self.__conductor_mwrank
394
395        elif algorithm == "all":
396            N1 = self.conductor("pari")
397            N2 = self.conductor("mwrank")
398            N3 = self.conductor("gp")
399            if N1 != N2 or N2 != N3:
400                raise ArithmeticError, "Pari, mwrank and gp compute different conductors (%s,%s,%s) for %s"%(
401                    N1, N2, N3, self)
402            return N1
403        else:
404            raise RuntimeError, "algorithm '%s' is not known."%algorithm
405
406    ####################################################################
407    #  Access to PARI curves related to this curve.
408    ####################################################################
409
410    def pari_curve(self, prec = None, factor = 1):
411        """
412        Return the PARI curve corresponding to this elliptic curve.
413
414        INPUT:
415            prec -- The precision of quantities calculated for the
416                    returned curve (in decimal digits).  if None, defaults
417                    to factor * the precision of the largest cached curve
418                    (or 10 if none yet computed)
419            factor -- the factor to increase the precision over the
420                      maximum previously computed precision.  Only used if
421                      prec (which gives an explicit precision) is None.
422               
423        EXAMPLES:
424            sage: E = EllipticCurve([0, 0,1,-1,0])
425            sage: e = E.pari_curve()
426            sage: type(e)
427            <type 'sage.libs.pari.gen.gen'>
428            sage: e.type()
429            't_VEC'
430            sage: e.ellan(10)
431            [1, -2, -3, 2, -2, 6, -1, 0, 6, 4]
432
433            sage: E = EllipticCurve(RationalField(), ['1/3', '2/3'])
434            sage: e = E.pari_curve()
435            sage: e.type()
436            't_VEC'
437            sage: e[:5]
438            [0, 0, 0, 1/3, 2/3]
439        """
440        if prec is None:
441            try:
442                L = self._pari_curve.keys()
443                L.sort()
444                if factor == 1:
445                    return self._pari_curve[L[-1]]
446                else:
447                    prec = int(factor * L[-1])
448                    self._pari_curve[prec] = pari(self.a_invariants()).ellinit(precision=prec)
449                    return self._pari_curve[prec]
450            except AttributeError:
451                pass
452        try:
453            return self._pari_curve[prec]
454        except AttributeError:
455            prec = 10
456            self._pari_curve = {}
457        except KeyError:
458            pass
459        self._pari_curve[prec] = pari(self.a_invariants()).ellinit(precision=prec)
460        return self._pari_curve[prec]
461
462    def pari_mincurve(self, prec = None):
463        """
464        Return the PARI curve corresponding to a minimal model
465        for this elliptic curve.
466
467        INPUT:
468        prec -- The precision of quantities calculated for the
469                returned curve (in decimal digits).  if None, defaults
470                to the precision of the largest cached curve (or 10 if
471                none yet computed)
472
473        EXAMPLES:
474            sage: E = EllipticCurve(RationalField(), ['1/3', '2/3'])
475            sage: e = E.pari_mincurve()
476            sage: e[:5]
477            [0, 0, 0, 27, 486]
478            sage: E.conductor()
479            47232
480            sage: e.ellglobalred()
481            [47232, [1, 0, 0, 0], 2]
482        """
483        if prec is None:
484            try:
485                L = self._pari_mincurve.keys()
486                L.sort()
487                return self._pari_mincurve[L[len(L) - 1]]
488            except AttributeError:
489                pass
490        try:
491            return self._pari_mincurve[prec]
492        except AttributeError:
493            prec = 10
494            self._pari_mincurve = {}
495        except KeyError:
496            pass
497        e = self.pari_curve(prec)
498        mc, change = e.ellminimalmodel()
499        self._pari_mincurve[prec] = mc
500        # self.__min_transform = change
501        return mc
502
503    def database_curve(self):
504        """
505        Return the curve in the elliptic curve database isomorphic to
506        this curve, if possible.  Otherwise raise a RuntimeError
507        exception. 
508
509        EXAMPLES:
510            sage: E = EllipticCurve([0,1,2,3,4])
511            sage: E.database_curve()
512            Elliptic Curve defined by y^2  = x^3 + x^2 + 3*x + 5 over Rational Field
513
514        NOTES: The model of the curve in the database can be different
515               than the Weierstrass model for this curve, e.g.,
516               database models are always minimal.
517        """
518        try:
519            return self.__database_curve
520        except AttributeError:
521            misc.verbose("Looking up %s in the database."%self)
522            D = sage.databases.cremona.CremonaDatabase()
523            ainvs = self.minimal_model().ainvs()
524            try:
525                self.__database_curve = D.elliptic_curve_from_ainvs(self.conductor(), ainvs)
526            except RuntimeError:
527                raise RuntimeError, "Elliptic curve %s not in the database."%self
528            return self.__database_curve
529           
530
531    def Np(self, p):
532        """
533        The number of points on E modulo p, where p is a prime, not
534        necessarily of good reduction.  (When p is a bad prime, also
535        counts the singular point.)
536
537        EXAMPLES:
538            sage: E = EllipticCurve([0, -1, 1, -10, -20])
539            sage: E.Np(2)
540            5
541            sage: E.Np(3)
542            5
543            sage: E.conductor()
544            11
545            sage: E.Np(11)
546            11
547        """
548        if self.conductor() % p == 0:
549            return p + 1 - self.ap(p)
550        #raise ArithmeticError, "p (=%s) must be a prime of good reduction"%p
551        if p < 1125899906842624:   # TODO: choose more wisely?
552            return p+1 - self.ap(p)
553        else:
554            return self.sea(p)
555
556    def sea(self, p, early_abort=False):
557        r"""
558        Return the number of points on $E$ over $\F_p$ computed using
559        the SEA algorithm, as implemented in PARI by Christophe Doche
560        and Sylvain Duquesne.
561
562        INPUT:
563            p -- a prime number
564            early_abort -- bool (default: Falst); if True an early abort technique
565                       is used and the computation is interrupted as soon
566                       as a small divisor of the order is detected.
567
568        \note{As of 2006-02-02 this function does not work on
569        Microsoft Windows under Cygwin (though it works under
570        vmware of course).}
571
572        EXAMPLES:
573            sage: E = EllipticCurve('37a')
574            sage: E.sea(next_prime(10^30))
575            1000000000000001426441464441649
576        """
577        p = rings.Integer(p)
578        return sea.ellsea(self.minimal_model().a_invariants(), p, early_abort=early_abort)
579
580    #def __pari_double_prec(self):
581    #    EllipticCurve_number_field._EllipticCurve__pari_double_prec(self)
582    #    try:
583    #        del self._pari_mincurve
584    #    except AttributeError:
585    #        pass
586       
587    ####################################################################
588    #  Access to mwrank
589    ####################################################################
590    def mwrank_curve(self, verbose=False):
591        """
592        Construct an mwrank_EllipticCurve from this elliptic curve
593
594        The resulting mwrank_EllipticCurve has available methods from
595        John Cremona's eclib library.
596
597        EXAMPLES:
598            sage: E=EllipticCurve('11a1')
599            sage: EE=E.mwrank_curve()
600            sage: EE
601            y^2+ y = x^3 - x^2 - 10*x - 20
602            sage: type(EE)
603            <class 'sage.libs.mwrank.interface.mwrank_EllipticCurve'>
604            sage: EE.isogeny_class()
605            ([[0, -1, 1, -10, -20], [0, -1, 1, -7820, -263580], [0, -1, 1, 0, 0]],
606            [[0, 5, 5], [5, 0, 0], [5, 0, 0]])
607        """
608        try:
609            return self.__mwrank_curve
610        except AttributeError:
611            pass
612        self.__mwrank_curve = mwrank.mwrank_EllipticCurve(
613            self.ainvs(), verbose=verbose)
614        return self.__mwrank_curve
615
616    def two_descent(self, verbose=True,
617                    selmer_only = False,
618                    first_limit = 20,
619                    second_limit = 8,
620                    n_aux = -1,
621                    second_descent = 1):
622        """
623        Compute 2-descent data for this curve.
624
625        INPUT:
626            verbose     -- (default: True) print what mwrank is doing
627                           If False, *no output* is
628            selmer_only -- (default: False) selmer_only switch
629            first_limit -- (default: 20) firstlim is bound on |x|+|z|
630            second_limit-- (default: 8)  secondlim is bound on log max {|x|,|z| },
631                                         i.e. logarithmic
632            n_aux       -- (default: -1) n_aux only relevant for general
633                           2-descent when 2-torsion trivial; n_aux=-1 causes default
634                           to be used (depends on method)
635            second_descent -- (default: True) second_descent only relevant for
636                           descent via 2-isogeny
637        OUTPUT:
638
639            Nothing -- nothing is returned (though much is printed
640                                            unless verbose=False)
641
642        EXAMPLES:
643            sage: E=EllipticCurve('37a1')
644            sage: E.two_descent(verbose=False) # no output
645        """
646        self.mwrank_curve().two_descent(verbose, selmer_only,
647                                        first_limit, second_limit,
648                                        n_aux, second_descent)
649
650
651    ####################################################################
652    #  Etc.
653    ####################################################################
654
655    def aplist(self, n, python_ints=False):
656        r"""
657        The Fourier coefficients $a_p$ of the modular form attached to
658        this elliptic curve, for all primes $p\leq n$.
659
660        INPUT:
661            n -- integer
662            python_ints -- bool (default: False); if True return a list of
663                      Python ints instead of SAGE integers.
664
665        OUTPUT:
666            -- list of integers
667
668        EXAMPLES:
669            sage: e = EllipticCurve('37a')
670            sage: e.aplist(1)
671            []
672            sage: e.aplist(2)
673            [-2]
674            sage: e.aplist(10)
675            [-2, -3, -2, -1]
676            sage: v = e.aplist(13); v
677            [-2, -3, -2, -1, -5, -2]
678            sage: type(v[0])
679            <type 'sage.rings.integer.Integer'>
680            sage: type(e.aplist(13, python_ints=True)[0])
681            <type 'int'>
682        """
683        # How is this result dependant on the real precision in pari?  At all?
684        e = self.pari_mincurve()
685        v = e.ellaplist(n, python_ints=True)
686        if python_ints:
687            return v
688        else:
689            return [Integer(a) for a in v]
690       
691       
692
693    def anlist(self, n, python_ints=False):
694        """
695        The Fourier coefficients up to and including $a_n$ of the
696        modular form attached to this elliptic curve.  The ith element
697        of the return list is a[i].
698
699        INPUT:
700            n -- integer
701            python_ints -- bool (default: False); if True return a list of
702                      Python ints instead of SAGE integers.
703
704        OUTPUT:
705            -- list of integers
706
707        EXAMPLES:
708            sage: E = EllipticCurve([0, -1, 1, -10, -20])
709            sage: E.anlist(3)
710            [0, 1, -2, -1]
711           
712            sage: E = EllipticCurve([0,1])
713            sage: E.anlist(20)
714            [0, 1, 0, 0, 0, 0, 0, -4, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 8, 0]
715        """
716        n = int(n)
717        # How is this result dependent on the real precision in Pari?  At all?
718        e = self.pari_mincurve()
719        if n >= 2147483648:
720            raise RuntimeError, "anlist: n (=%s) must be < 2147483648."%n
721
722        v = [0] + e.ellan(n, python_ints=True)
723        if not python_ints:
724            v = [Integer(x) for x in v]
725        return v
726       
727       
728        # There is some overheard associated with coercing the PARI
729        # list back to Python, but it's not bad.  It's better to do it
730        # this way instead of trying to eval the whole list, since the
731        # int conversion is done very sensibly.  NOTE: This would fail
732        # if a_n won't fit in a C int, i.e., is bigger than
733        # 2147483648; however, we wouldn't realistically compute
734        # anlist for n that large anyways.
735        #
736        # Some relevant timings:
737        #
738        # E <--> [0, 1, 1, -2, 0]   389A
739        #  E = EllipticCurve([0, 1, 1, -2, 0]);   // SAGE or MAGMA
740        #  e = E.pari_mincurve()
741        #  f = ellinit([0,1,1,-2,0]);
742        #
743        #  Computation                                              Time (1.6Ghz Pentium-4m laptop)
744        #  time v:=TracesOfFrobenius(E,10000);  // MAGMA            0.120
745        #  gettime;v=ellan(f,10000);gettime/1000                    0.046
746        #  time v=e.ellan (10000)                                   0.04
747        #  time v=E.anlist(10000)                                   0.07
748       
749        #  time v:=TracesOfFrobenius(E,100000);  // MAGMA           1.620
750        #  gettime;v=ellan(f,100000);gettime/1000                   0.676
751        #  time v=e.ellan (100000)                                  0.7
752        #  time v=E.anlist(100000)                                  0.83
753
754        #  time v:=TracesOfFrobenius(E,1000000);  // MAGMA          20.850
755        #  gettime;v=ellan(f,1000000);gettime/1000                  9.238
756        #  time v=e.ellan (1000000)                                 9.61
757        #  time v=E.anlist(1000000)                                 10.95  (13.171 in cygwin vmware)
758       
759        #  time v:=TracesOfFrobenius(E,10000000);  //MAGMA          257.850
760        #  gettime;v=ellan(f,10000000);gettime/1000      FAILS no matter how many allocatemem()'s!!
761        #  time v=e.ellan (10000000)                                139.37 
762        #  time v=E.anlist(10000000)                                136.32
763        #
764        #  The last SAGE comp retries with stack size 40MB,
765        #  80MB, 160MB, and succeeds last time.  It's very interesting that this
766        #  last computation is *not* possible in GP, but works in py_pari!
767        #
768
769    def q_expansion(self, prec):
770        r"""
771        Return the $q$-expansion to precision prec of the newform
772        attached to this elliptic curve.
773
774        INPUT:
775            prec -- an integer
776
777        OUTPUT:
778            a power series (in th evariable 'q')
779
780        NOTE: If you want the output to be a modular form and not just
781        a $q$-expansion, use \code{self.modular_form()}.
782
783        EXAMPLES:
784            sage: E=EllipticCurve('37a1')
785            sage: E.q_expansion(20)
786            q - 2*q^2 - 3*q^3 + 2*q^4 - 2*q^5 + 6*q^6 - q^7 + 6*q^9 + 4*q^10 - 5*q^11 - 6*q^12 - 2*q^13 + 2*q^14 + 6*q^15 - 4*q^16 - 12*q^18 + O(q^20)
787
788        """
789        return PowerSeriesRing(Q, 'q')(self.anlist(prec), prec, check=True)
790
791    def modular_form(self):
792        r"""
793        Return the cuspidal modular form associated to this elliptic curve.
794
795        EXAMPLES:
796            sage: E = EllipticCurve('37a')
797            sage: f = E.modular_form()
798            sage: f
799            q - 2*q^2 - 3*q^3 + 2*q^4 - 2*q^5 + O(q^6)
800           
801            If you need to see more terms in the $q$-expansion:
802            sage: f.q_expansion(20)
803            q - 2*q^2 - 3*q^3 + 2*q^4 - 2*q^5 + 6*q^6 - q^7 + 6*q^9 + 4*q^10 - 5*q^11 - 6*q^12 - 2*q^13 + 2*q^14 + 6*q^15 - 4*q^16 - 12*q^18 + O(q^20)
804
805           
806        NOTE: If you just want the $q$-expansion, use
807        \code{self.q_expansion(prec)}.
808        """
809        try:
810            return self.__modular_form
811        except AttributeError:
812            M = sage.modular.modform.constructor.ModularForms(self.conductor(),weight=2)
813            f = sage.modular.modform.element.ModularFormElement_elliptic_curve(M, self)
814            self.__modular_form = f
815            return f
816
817    def modular_symbol_space(self, sign=1, base_ring=Q, bound=None):
818        r"""
819        Return the space of cuspidal modular symbols associated to
820        this elliptic curve, with given sign and base ring.
821
822        INPUT:
823            sign -- 0, -1, or 1
824            base_ring -- a ring
825
826        EXAMPLES:
827            sage: f = EllipticCurve('37b')
828            sage: f.modular_symbol_space()
829            Modular Symbols subspace of dimension 1 of Modular Symbols space of dimension 3 for Gamma_0(37) of weight 2 with sign 1 over Rational Field
830            sage: f.modular_symbol_space(-1)
831            Modular Symbols subspace of dimension 1 of Modular Symbols space of dimension 2 for Gamma_0(37) of weight 2 with sign -1 over Rational Field
832            sage: f.modular_symbol_space(0, bound=3)
833            Modular Symbols subspace of dimension 2 of Modular Symbols space of dimension 5 for Gamma_0(37) of weight 2 with sign 0 over Rational Field
834       
835        NOTE: If you just want the $q$-expansion, use
836        \code{self.q_expansion(prec)}.
837        """
838        typ = (sign, base_ring)
839        try:
840            return self.__modular_symbol_space[typ]
841        except AttributeError:
842            self.__modular_symbol_space = {}
843        except KeyError:
844            pass
845        M = ell_modular_symbols.modular_symbol_space(self, sign, base_ring, bound=bound)
846        self.__modular_symbol_space[typ] = M
847        return M
848
849    def modular_symbol(self, sign=1, normalize=True):
850        r"""
851        Return the modular symbol associated to this elliptic curve,
852        with given sign and base ring.  This is the map that sends r/s
853        to a fixed multiple of 2*pi*I*f(z)dz from oo to r/s,
854        normalized so that all values of this map take values in QQ.
855
856        If sign=1, the normalization is such that the p-adic
857        L-function associated to this modular symbol is correct.
858        I.e., the normalization is the same as for the integral period
859        mapping divided by 2.
860
861        INPUT:
862            sign -- -1, or 1
863            base_ring -- a ring
864            normalize -- (default: True); if True, the modular symbol
865                is correctly normalized (up to possibly a factor of
866                -1 or 2).  If False, the modular symbol is almost certainly
867                not correctly normalized, i.e., all values will be a
868                fixed scalar multiple of what they should be.  But
869                the initial computation of the modular symbol is
870                much faster, though evaluation of it after computing
871                it won't be any faster.
872
873        EXAMPLES:
874            sage: E=EllipticCurve('37a1')
875            sage: M=E.modular_symbol(); M
876            Modular symbol with sign 1 over Rational Field attached to Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field
877            sage: M(1/2)
878            0
879            sage: M(1/5)
880            1/2
881        """
882        typ = (sign, normalize)
883        try:
884            return self.__modular_symbol[typ]
885        except AttributeError:
886            self.__modular_symbol = {}
887        except KeyError:
888            pass
889        M = ell_modular_symbols.ModularSymbol(self, sign, normalize)
890        self.__modular_symbol[typ] = M
891        return M
892
893    padic_lseries = padics.padic_lseries
894
895    def newform(self):
896        r"""
897        Same as \code{self.modular_form()}.
898
899        EXAMPLES:
900            sage: E=EllipticCurve('37a1')
901            sage: E.newform()
902            q - 2*q^2 - 3*q^3 + 2*q^4 - 2*q^5 + O(q^6)
903            sage: E.newform() == E.modular_form()
904            True
905        """
906        return self.modular_form()
907
908    def q_eigenform(self, prec):
909        r"""
910        Synonym for \code{self.q_expansion(prec)}.
911
912        EXAMPLES:
913            sage: E=EllipticCurve('37a1')
914            sage: E.q_eigenform(10)
915            q - 2*q^2 - 3*q^3 + 2*q^4 - 2*q^5 + 6*q^6 - q^7 + 6*q^9 + O(q^10)
916            sage: E.q_eigenform(10) == E.q_expansion(10)
917            True
918        """
919        return self.q_expansion(prec)
920   
921    def analytic_rank(self, algorithm="cremona"):
922        r"""
923        Return an integer that is \emph{probably} the analytic rank of
924        this elliptic curve. 
925
926        INPUT:
927            algorithm:
928                -- 'cremona' (default) --  Use the Buhler-Gross algorithm
929                    as implemented in GP by Tom Womack and John Cremona,
930                    who note that their implementation is practical for
931                    any rank and conductor $\leq 10^{10}$ in 10 minutes.
932
933                -- 'sympow' --use Watkins's program sympow
934
935                -- 'rubinstein' -- use Rubinstein's L-function C++ program lcalc.
936
937                -- 'magma' -- use MAGMA
938
939                -- 'all' -- compute with all other free algorithms, check that the
940                            answers agree, and return the common answer.
941
942        \note{If the curve is loaded from the large Cremona database,
943        then the modular degree is taken from the database.}
944
945        Of the three above, probably Rubinstein's is the most
946        efficient (in some limited testing I've done).
947
948        \note{It is an open problem to \emph{prove} that \emph{any}
949        particular elliptic curve has analytic rank $\geq 4$.}
950
951        EXAMPLES:
952            sage: E = EllipticCurve('389a')
953            sage: E.analytic_rank(algorithm='cremona')
954            2
955            sage: E.analytic_rank(algorithm='rubinstein')
956            2
957            sage: E.analytic_rank(algorithm='sympow')
958            2
959            sage: E.analytic_rank(algorithm='magma')    # optional
960            2
961            sage: E.analytic_rank(algorithm='all')
962            2
963        """
964        if algorithm == 'cremona':
965            return rings.Integer(gp_cremona.ellanalyticrank(self.minimal_model().a_invariants()))
966        elif algorithm == 'rubinstein':
967            from sage.lfunctions.lcalc import lcalc
968            return lcalc.analytic_rank(L=self)
969        elif algorithm == 'sympow':
970            from sage.lfunctions.sympow import sympow
971            return sympow.analytic_rank(self)[0]
972        elif algorithm == 'magma':
973            return rings.Integer(self._magma_().AnalyticRank())
974        elif algorithm == 'all':
975            S = list(set([self.analytic_rank('cremona'), 
976                     self.analytic_rank('rubinstein'), self.analytic_rank('sympow')]))
977            if len(S) != 1:
978                raise RuntimeError, "Bug in analytic rank; algorithms don't agree! (E=%s)"%E
979            return S[0]
980        else:
981            raise ValueError, "algorithm %s not defined"%algorithm
982           
983    def p_isogenous_curves(self, p=None):
984        r"""
985        Return a list of pairs $(p, L)$ where $p$ is a prime and $L$
986        is a list of the elliptic curves over $\Q$ that are
987        $p$-isogenous to this elliptic curve.
988
989        INPUT:
990            p -- prime or None (default: None); if a prime, returns
991                 a list of the p-isogenous curves.  Otherwise, returns
992                 a list of all prime-degree isogenous curves sorted
993                 by isogeny degree.
994
995        This is implemented using Cremona's GP script \code{allisog.gp}.
996
997        EXAMPLES:
998            sage: E = EllipticCurve([0,-1,0,-24649,1355209])   
999            sage: E.p_isogenous_curves()
1000            [(2, [Elliptic Curve defined by y^2  = x^3 - x^2 - 91809*x - 9215775 over Rational Field, Elliptic Curve defined by y^2  = x^3 - x^2 - 383809*x + 91648033 over Rational Field, Elliptic Curve defined by y^2  = x^3 - x^2 + 1996*x + 102894 over Rational Field])]
1001
1002        The isogeny class of the curve 11a2 has three curves in it.
1003        But \code{p_isogenous_curves} only returns one curves, since
1004        there is only one curve $5$-isogenous to 11a2.
1005            sage: E = EllipticCurve('11a2')
1006            sage: E.p_isogenous_curves()
1007            [(5, [Elliptic Curve defined by y^2 + y = x^3 - x^2 - 10*x - 20 over Rational Field])]
1008            sage: E.p_isogenous_curves(5)
1009            [Elliptic Curve defined by y^2 + y = x^3 - x^2 - 10*x - 20 over Rational Field]
1010            sage: E.p_isogenous_curves(3)
1011            []
1012
1013        In contrast, the curve 11a1 admits two $5$-isogenies:
1014            sage: E = EllipticCurve('11a1')
1015            sage: E.p_isogenous_curves(5)
1016            [Elliptic Curve defined by y^2 + y = x^3 - x^2 - 7820*x - 263580 over Rational Field,
1017             Elliptic Curve defined by y^2 + y = x^3 - x^2 over Rational Field]
1018        """
1019        if p is None:
1020            X = eval(gp_cremona.allisog(self.minimal_model().a_invariants()))
1021            Y = [(p, [constructor.EllipticCurve(ainvs) for ainvs in L]) for p, L in X]
1022            Y.sort()
1023            return Y
1024        else:
1025            X = eval(gp_cremona.p_isog(self.minimal_model().a_invariants(), p))
1026            Y = [constructor.EllipticCurve(ainvs) for ainvs in X]
1027            Y.sort()
1028            return Y
1029
1030    def simon_two_descent(self, verbose=0, lim1=5, lim3=50, limtriv=10, maxprob=20, limbigprime=30):
1031        r"""
1032        Given a curve with no 2-torsion, computes (probably) the rank
1033        of the Mordell-Weil group, with certainty the rank of the
1034        2-Selmer group, and a list of independent points on the curve.
1035
1036        INPUT:
1037            verbose -- integer, 0,1,2,3; (default: 0), the verbosity level
1038            lim1    -- (default: 5) limite des points triviaux sur les quartiques
1039            lim3    -- (default: 50) limite des points sur les quartiques ELS
1040            limtriv -- (default: 10) limite des points triviaux sur la
1041                                     courbe elliptique
1042            maxprob -- (default: 20)
1043            limbigprime -- (default: 30)  to distinguish between small and large prime
1044                                          numbers. Use probabilistic tests for large
1045                                          primes. If 0, don't any probabilistic tests.
1046                           
1047        OUTPUT:
1048            integer -- "probably" the rank of self
1049            integer -- the 2-rank of the Selmer group
1050            list    -- list of independent points on the curve.
1051
1052        IMPLEMENTATION: Uses {\bf Denis Simon's} GP/PARI scripts from
1053                         \url{http://www.math.unicaen.fr/~simon/}
1054
1055        EXAMPLES:
1056        These computations use pseudo-random numbers, so we set the
1057        seed for reproducible testing.
1058
1059        We compute the ranks of the curves of lowest known conductor up to rank $8$.
1060        Amazingly, each of these computations finishes almost instantly!
1061           
1062            sage: E = EllipticCurve('11a1')
1063            sage: set_random_seed(0)
1064            sage: E.simon_two_descent()
1065            (0, 0, [])
1066            sage: E = EllipticCurve('37a1')
1067            sage: set_random_seed(0)
1068            sage: E.simon_two_descent()
1069            (1, 1, [(0 : 0 : 1)])
1070            sage: E = EllipticCurve('389a1')
1071            sage: set_random_seed(0)
1072            sage: E.simon_two_descent()
1073            (2, 2, [(1 : 0 : 1), (-11/9 : -55/27 : 1)])
1074            sage: E = EllipticCurve('5077a1')
1075            sage: set_random_seed(0)
1076            sage: E.simon_two_descent()
1077            (3, 3, [(1 : 0 : 1), (2 : -1 : 1), (0 : 2 : 1)])
1078
1079
1080        In this example Simon's program does not find any points, though
1081        it does correctly compute the rank of the 2-Selmer group.
1082            sage: E = EllipticCurve([1, -1, 0, -751055859, -7922219731979])     # long (0.6 seconds)
1083            sage: set_random_seed(0)
1084            sage: E.simon_two_descent ()
1085            (1, 1, [])           
1086
1087        The rest of these entries were taken from Tom Womack's page
1088        \url{http://tom.womack.net/maths/conductors.htm}
1089
1090            sage: E = EllipticCurve([1, -1, 0, -79, 289])
1091            sage: set_random_seed(0)
1092            sage: E.simon_two_descent()
1093            (4, 4, [(6 : -5 : 1), (4 : 3 : 1), (5 : -3 : 1), (8 : -15 : 1)])
1094            sage: E = EllipticCurve([0, 0, 1, -79, 342])
1095            sage: set_random_seed(0)
1096            sage: E.simon_two_descent()
1097            (5, 5, [(5 : 8 : 1), (10 : 23 : 1), (3 : 11 : 1), (4 : -10 : 1), (0 : 18 : 1)])
1098            sage: E = EllipticCurve([1, 1, 0, -2582, 48720])
1099            sage: set_random_seed(0)
1100            sage: r, s, G = E.simon_two_descent(); r,s
1101            (6, 6)
1102            sage: E = EllipticCurve([0, 0, 0, -10012, 346900])
1103            sage: set_random_seed(0)
1104            sage: r, s, G = E.simon_two_descent(); r,s
1105            (7, 7)
1106            sage: E = EllipticCurve([0, 0, 1, -23737, 960366])   
1107            sage: set_random_seed(0)
1108            sage: r, s, G = E.simon_two_descent(); r,s
1109            (8, 8)
1110        """
1111        t = simon_two_descent(self, verbose=verbose, lim1=lim1, lim3=lim3, limtriv=limtriv,
1112                              maxprob=maxprob, limbigprime=limbigprime)
1113        prob_rank = rings.Integer(t[0])
1114        two_selmer_rank = rings.Integer(t[1])
1115        prob_gens = [self(P) for P in t[2]]
1116        return prob_rank, two_selmer_rank, prob_gens
1117
1118    two_descent_simon = simon_two_descent
1119
1120    def three_selmer_rank(self, bound=0, method=2):
1121        r"""
1122        Return the 3-selmer rank of this elliptic curve, computed
1123        using Magma.
1124
1125        This is not implemented for all curves; a NotImplementedError
1126        exception is raised when this function is called on curves
1127        for which 3-descent isn't implemented.
1128       
1129        \note{Use a slightly modified version of Michael Stoll's MAGMA
1130        file \code{3descent.m}.  You must have Magma to use this
1131        function.}
1132
1133        EXAMPLES:
1134            sage: EllipticCurve('37a').three_selmer_rank()  # optional & long -- Magma
1135            1
1136
1137            sage: EllipticCurve('14a1').three_selmer_rank()      # optional
1138            Traceback (most recent call last):
1139            ...
1140            NotImplementedError:  Currently, only the case with irreducible phi3 is implemented.
1141        """
1142        import magma_3descent
1143        try:
1144            return magma_3descent.three_selmer_rank(self, bound, method)
1145        except RuntimeError, msg:
1146            msg = str(msg)
1147            i = msg.rfind(':')
1148            raise NotImplementedError, msg[i+1:]
1149
1150
1151    def rank(self, use_database=False, verbose=False,
1152                   only_use_mwrank=True,
1153                   algorithm='mwrank_shell',
1154                   proof=None):
1155        """
1156        Return the rank of this elliptic curve, assuming no conjectures.
1157
1158        If we fail to provably compute the rank, raises a RuntimeError
1159        exception.
1160
1161        INPUT:
1162            use_database (bool) -- (default: False), if True, try to
1163                  look up the regulator in the Cremona database.
1164            verbose -- (default: None), if specified changes the
1165                       verbosity of mwrank computations.
1166            algorithm -- 'mwrank_shell' -- call mwrank shell command
1167                      -- 'mwrank_lib' -- call mwrank c library
1168            only_use_mwrank -- (default: True) if False try using
1169                       analytic rank methods first.
1170            proof -- bool or None (default: None, see proof.elliptic_curve or
1171                       sage.structure.proof).  Note that results
1172                       obtained from databases are considered proof = True
1173           
1174        OUTPUT:
1175            rank (int) -- the rank of the elliptic curve.
1176
1177        IMPLEMENTATION: Uses L-functions, mwrank, and databases.
1178
1179        EXAMPLES:
1180            sage: EllipticCurve('11a').rank()
1181            0
1182            sage: EllipticCurve('37a').rank()
1183            1
1184            sage: EllipticCurve('389a').rank()
1185            2
1186            sage: EllipticCurve('5077a').rank()
1187            3
1188            sage: EllipticCurve([1, -1, 0, -79, 289]).rank()   # long time.  This will use the default proof behavior of True.
1189            4
1190            sage: EllipticCurve([0, 0, 1, -79, 342]).rank(proof=False)  # long time -- but under a minute
1191            5
1192            sage: EllipticCurve([0, 0, 1, -79, 342]).simon_two_descent()[0]  # much faster -- almost instant.
1193            5
1194
1195        Examples with denominators in defining equations:
1196            sage: E = EllipticCurve( [0, 0, 0, 0, -675/4])
1197            sage: E.rank()
1198            0
1199            sage: E = EllipticCurve( [0, 0, 1/2, 0, -1/5])
1200            sage: E.rank()
1201            1
1202            sage: E.minimal_model().rank()
1203            1
1204        """
1205        if proof is None:
1206            from sage.structure.proof.proof import get_flag
1207            proof = get_flag(proof, "elliptic_curve")
1208        else:
1209            proof = bool(proof)
1210        try:
1211            return self.__rank[proof]
1212        except KeyError:
1213            if proof is False and self.__rank.has_key(True):
1214                return self.__rank[True]
1215        if use_database:
1216            try:
1217                self.__rank[True] = self.database_curve().rank()
1218                return self.__rank[True]
1219            except (AttributeError, RuntimeError):
1220                pass
1221        if not only_use_mwrank:
1222            N = self.conductor()
1223            prec = int(4*float(sqrt(N))) + 10
1224            if self.root_number() == 1:
1225                L, err = self.lseries().at1(prec)           
1226                if abs(L) > err + R(0.0001):  # definitely doesn't vanish
1227                    misc.verbose("rank 0 because L(E,1)=%s"%L)
1228                    self.__rank[proof] = 0
1229                    return self.__rank[proof]
1230            else:
1231                Lprime, err = self.lseries().deriv_at1(prec)
1232                if abs(Lprime) > err + R(0.0001):  # definitely doesn't vanish
1233                    misc.verbose("rank 1 because L'(E,1)=%s"%Lprime)
1234                    self.__rank[proof] = 1
1235                    return self.__rank[proof]
1236
1237        if algorithm == 'mwrank_lib':
1238            misc.verbose("using mwrank lib")
1239            C = self.mwrank_curve()
1240            C.set_verbose(verbose)
1241            r = C.rank()
1242            if not C.certain():
1243                del self.__mwrank_curve
1244                raise RuntimeError, "Unable to compute the rank with certainty (lower bound=%s).  This could be because Sha(E/Q)[2] is nontrivial."%C.rank() + "\nTrying calling something like two_descent(second_limit=13) on the curve then trying this command again.  You could also try rank with only_use_mwrank=False."
1245            self.__rank[proof] = r
1246        elif algorithm == 'mwrank_shell':
1247            misc.verbose("using mwrank shell")
1248            X = self.mwrank()
1249            if not 'The rank and full Mordell-Weil basis have been determined unconditionally' in X:
1250                if proof:
1251                    raise RuntimeError, '%s\nRank not provably correct.'%X
1252                else:
1253                    misc.verbose("Warning -- rank not provably correct", level=1)
1254            elif proof is False:
1255                proof = True #since we actually provably found the rank
1256            i = X.find('Rank = ')
1257            assert i != -1
1258            j = i + X[i:].find('\n')
1259            self.__rank[proof] = Integer(X[i+7:j])
1260        return self.__rank[proof]
1261
1262    def gens(self, verbose=False, rank1_search=10,
1263             algorithm='mwrank_shell',
1264             only_use_mwrank=True,
1265             proof = None):
1266        """
1267        Compute and return generators for the Mordell-Weil group E(Q)
1268        *modulo* torsion.
1269
1270        HINT: If you would like to control the height bounds used in
1271        the 2-descent, first call the two_descent function with those
1272        height bounds. However that function, while it displays a lot
1273        of output, returns no values.
1274
1275        TODO: (1) Right now this function assumes that the input curve
1276        is in minimal Weierstrass form.  This restriction will be
1277        removed in the future.  This function raises a
1278        NotImplementedError if a non-minimal curve is given as input.
1279
1280              (2) Allow passing of command-line parameters to mwrank.
1281
1282        WARNING: If the program fails to give a provably correct
1283        result, it prints a warning message, but does not raise an
1284        exception.  Use the gens_certain command to find out if
1285        this warning message was printed.
1286       
1287        INPUT:
1288            verbose -- (default: None), if specified changes the
1289                       verbosity of mwrank computations.
1290            rank1_search -- (default: 16), if the curve has analytic
1291                       rank 1, try to find a generator by a direct
1292                       search up to this logarithmic height.  If this
1293                       fails the usual mwrank procedure is called.
1294            algorithm -- 'mwrank_shell' (default) -- call mwrank shell command
1295                      -- 'mwrank_lib' -- call mwrank c library
1296            only_use_mwrank -- bool (default True) if false, attempts to
1297                       first use more naive, natively implemented methods.
1298            proof -- bool or None (default None, see proof.elliptic_curve or
1299                       sage.structure.proof). 
1300        OUTPUT:
1301            generators -- List of generators for the Mordell-Weil group.
1302
1303        IMPLEMENTATION: Uses Cremona's mwrank C library.
1304
1305        EXAMPLES:
1306            sage: E = EllipticCurve('389a')
1307            sage: E.gens()                 # random output   
1308            [(-1 : 1 : 1), (0 : 0 : 1)]
1309        """
1310        if proof is None:
1311            from sage.structure.proof.proof import get_flag
1312            proof = get_flag(proof, "elliptic_curve")
1313        else:
1314            proof = bool(proof)
1315        try:
1316            return list(self.__gens[proof])  # return copy so not changed
1317        except KeyError:
1318            if proof is False and self.__gens.has_key(True):
1319                return self.__gens[True]
1320        if self.conductor() > 10**7:
1321            only_use_mwrank = True
1322
1323        if not only_use_mwrank:
1324            try:
1325                misc.verbose("Trying to compute rank.")
1326                r = self.rank(only_use_mwrank = False)
1327                misc.verbose("Got r = %s."%r)
1328                if r == 0:
1329                    misc.verbose("Rank = 0, so done.")
1330                    self.__gens[True] = []
1331                    self.__regulator[True] = R(1)
1332                    return self.__gens[True]
1333                if r == 1 and rank1_search:
1334                    misc.verbose("Rank = 1, so using direct search.")
1335                    h = 6
1336                    while h <= rank1_search:
1337                        misc.verbose("Trying direct search up to height %s"%h)
1338                        G = self.point_search(h, verbose)
1339                        G = [P for P in G if P.order() == oo]
1340                        if len(G) > 0:
1341                            misc.verbose("Direct search succeeded.")
1342                            G, _, reg = self.saturation(G, verbose=verbose)
1343                            misc.verbose("Computed saturation.")
1344                            self.__gens[True] = G
1345                            self.__regulator[True] = reg
1346                            return self.__gens[True]
1347                        h += 2
1348                    misc.verbose("Direct search FAILED.")
1349            except RuntimeError:
1350                pass
1351        # end if (not_use_mwrank)
1352        if not self.is_integral():
1353            raise NotImplementedError, "gens via mwrank only implemented for curves with integer coefficients."
1354        if algorithm == "mwrank_lib":
1355            misc.verbose("Calling mwrank C++ library.")
1356            C = self.mwrank_curve(verbose)
1357            if not (verbose is None):
1358                C.set_verbose(verbose)
1359            G = C.gens()
1360            if proof is True and C.certain() is False:
1361                del self.__mwrank_curve
1362                raise RuntimeError, "Unable to compute the rank, hence generators, with certainty (lower bound=%s).  This could be because Sha(E/Q)[2] is nontrivial."%C.rank() + \
1363                      "\nTrying calling something like two_descent(second_limit=13) on the curve then trying this command again."
1364            else:
1365                proof = C.certain()
1366        else:
1367            X = self.mwrank()
1368            misc.verbose("Calling mwrank shell.")
1369            if not 'The rank and full Mordell-Weil basis have been determined unconditionally' in X:
1370                msg = 'Generators not provably computed.'
1371                if proof:
1372                    raise RuntimeError, '%s\n%s'%(X,msg)
1373                else:
1374                    misc.verbose("Warning -- %s"%msg, level=1)
1375            elif proof is False:
1376                proof = True
1377            G = []
1378            i = X.find('Generator ')
1379            while i != -1:
1380                j = i + X[i:].find(';')
1381                k = i + X[i:].find('[')
1382                G.append(eval(X[k:j].replace(':',',')))
1383                X = X[j:]
1384                i = X.find('Generator ')
1385            i = X.find('Regulator = ')
1386            j = i + X[i:].find('\n')
1387            self.__regulator[proof] = R(X[i+len('Regulator = '):j])
1388        ####
1389        self.__gens[proof] = [self.point(x, check=True) for x in G]
1390        self.__gens[proof].sort()
1391        self.__rank[proof] = len(self.__gens[proof])
1392        return self.__gens[proof]
1393
1394    def gens_certain(self):
1395        """
1396        Return True if the generators have been proven correct.
1397
1398        EXAMPLES:
1399            sage: E=EllipticCurve('37a1')
1400            sage: E.gens()
1401            [(0 : -1 : 1)]
1402            sage: E.gens_certain()
1403            True
1404       """
1405        return self.__gens.has_key(True)
1406
1407    def ngens(self, proof = None):
1408        """
1409        Return the number of generators of this elliptic curve.
1410
1411        NOTE: See gens() for further documentation.  The function
1412        ngens() calls gens() if not already done, but only with
1413        default parameters.  Better results may be obtained by calling
1414        mwrank() with carefully chosen parameters.
1415
1416        EXAMPLES:
1417            sage: E=EllipticCurve('37a1')
1418            sage: E.ngens()
1419            1
1420
1421        TO DO: This example should not cause a run-time error.
1422            sage: E=EllipticCurve([0,0,0,877,0])
1423            sage: # E.ngens()  ######## causes run-time error
1424
1425            sage: print E.mwrank('-v0 -b12 -l')
1426            Curve [0,0,0,877,0] :   Rank = 1
1427            Generator 1 is [29604565304828237474403861024284371796799791624792913256602210:-256256267988926809388776834045513089648669153204356603464786949:490078023219787588959802933995928925096061616470779979261000]; height 95.980371987964
1428            Regulator = 95.980...
1429
1430
1431        """
1432        return len(self.gens(proof = proof))
1433
1434    def regulator(self, use_database=True, verbose=None, proof=None):
1435        """
1436        Returns the regulator of this curve, which must be defined
1437        over Q.
1438
1439        INPUT:
1440            use_database -- bool (default: False), if True, try to
1441                  look up the regulator in the Cremona database.
1442            verbose -- (default: None), if specified changes the
1443                  verbosity of mwrank computations.
1444            proof -- bool or None (default: None, see proof.[tab] or
1445                       sage.structure.proof).  Note that results from
1446                       databases are considered proof = True
1447
1448        EXAMPLES:
1449            sage: E = EllipticCurve([0, 0, 1, -1, 0])
1450            sage: E.regulator()              # long time (1 second)
1451            0.0511114082399688
1452            sage: EllipticCurve('11a').regulator()
1453            1.00000000000000
1454            sage: EllipticCurve('37a').regulator()
1455            0.0511114082399688
1456            sage: EllipticCurve('389a').regulator()
1457            0.152460177943144
1458            sage: EllipticCurve('5077a').regulator()    # random low order bit
1459            0.417143558758385
1460            sage: EllipticCurve([1, -1, 0, -79, 289]).regulator()  # long time (seconds)
1461            1.50434488827528
1462            sage: EllipticCurve([0, 0, 1, -79, 342]).regulator(proof=False)  # long time (seconds)
1463            14.7905275701310
1464        """
1465        if proof is None:
1466            from sage.structure.proof.proof import get_flag
1467            proof = get_flag(proof, "elliptic_curve")
1468        else:
1469            proof = bool(proof)
1470        try:
1471            return self.__regulator[proof]
1472        except KeyError:
1473            if proof is False and self.__regulator.has_key(True):
1474                return self.__regulator[True]
1475        if use_database:
1476            try:
1477                self.__regulator[True] = R(self.database_curve().db_extra[3])
1478                return self.__regulator[True]
1479            except (AttributeError, RuntimeError):
1480                pass
1481        G = self.gens(proof=proof)
1482        try:  # in some cases self.gens() efficiently computes regulator.
1483            return self.__regulator[proof]
1484        except KeyError:
1485            if proof is False and self.__regulator.has_key(True):
1486                return self.__regulator[True]
1487        C = self.mwrank_curve()
1488        reg = R(C.regulator())
1489        if proof is True and not C.certain():
1490            raise RuntimeError, "Unable to compute the rank, hence regulator, with certainty (lower bound=%s)."%C.rank()
1491        proof = C.certain()
1492        self.__regulator[proof] = reg
1493        return self.__regulator[proof]
1494
1495    def saturation(self, points, verbose=False, max_prime=0, odd_primes_only=False):
1496        """
1497        Given a list of rational points on E, compute the saturation
1498        in E(Q) of the subgroup they generate.
1499
1500        INPUT:
1501            points (list) -- list of points on E
1502            verbose  (bool) -- (default: False), if True, give verbose output
1503            max_prime (int) -- (default: 0), saturation is performed
1504                               for all primes up to max_prime.  If max_prime==0,
1505                               perform saturation at *all* primes, i.e., compute
1506                               the true saturation.
1507            odd_primes_only (bool) -- only do saturation at odd primes
1508
1509        OUTPUT:
1510            saturation (list) -- points that form a basis for the saturation
1511            index (int) -- the index of the group generated by points
1512                           in their saturation
1513            regulator (float) -- regulator of saturated points.
1514
1515        IMPLEMENTATION:
1516            Uses Cremona's mwrank package.  With max_prime=0, we call
1517            mwrank with successively larger prime bounds until the
1518            full saturation is provably found.  The results of
1519            saturation at the previous primes is stored in each case,
1520            so this should be reasonably fast.
1521
1522        EXAMPLES:
1523            sage: E=EllipticCurve('37a1')
1524            sage: P=E.gens()[0]
1525            sage: Q=5*P; Q
1526            (1/4 : -3/8 : 1)
1527            sage: E.saturation([Q])
1528            ([(0 : -1 : 1)], '5', 0.0511114075779915)
1529        """
1530        if not isinstance(points, list):
1531            raise TypeError, "points (=%s) must be a list."%points
1532
1533        v = []
1534        for P in points:
1535            if not isinstance(P, ell_point.EllipticCurvePoint_field):
1536                P = self(P)
1537            elif P.curve() != self:
1538                raise ArithmeticError, "point (=%s) must be %s."%(P,self)
1539            x, y = P.xy()
1540            d = x.denominator().lcm(y.denominator())
1541            v.append((x*d, y*d, d))
1542        c = self.mwrank_curve()
1543        mw = mwrank.mwrank_MordellWeil(c, verbose)
1544        mw.process(v)
1545        if max_prime == 0:
1546            repeat_until_saturated = True
1547            max_prime = 97
1548        while True:
1549            ok, index, unsat = mw.saturate(max_prime=max_prime, odd_primes_only = odd_primes_only)
1550            reg = mw.regulator()
1551            if not ok and repeat_until_saturated:
1552                max_prime = arith.next_prime(max_prime + 100)
1553                ok, index, unsat = mw.saturate(max_prime=max_prime, odd_primes_only = odd_primes_only)
1554                reg = mw.regulator()
1555            else:
1556                break
1557        sat = mw.points()
1558        sat = [self(P) for P in sat]
1559        return sat, index, R(reg)
1560
1561    def CPS_height_bound(self):
1562        """
1563        Return the Cremona-Prickett-Siksek height bound.  This is a
1564        floating point number B such that if P is a point on the curve,
1565        then the naive logarithmetic height of P is off from the
1566        canonical height by at most B.
1567
1568        EXAMPLES:
1569            sage: E = EllipticCurve("11a")
1570            sage: E.CPS_height_bound()
1571            2.8774743273580445
1572            sage: E = EllipticCurve("5077a")
1573            sage: E.CPS_height_bound()
1574            0.0
1575            sage: E = EllipticCurve([1,2,3,4,1])
1576            sage: E.CPS_height_bound()
1577            Traceback (most recent call last):
1578            ...
1579            RuntimeError: curve must be minimal.
1580            sage: F = E.quadratic_twist(-19)
1581            sage: F
1582            Elliptic Curve defined by y^2 + x*y + y = x^3 - x^2 + 1376*x - 130 over Rational Field
1583            sage: F.CPS_height_bound()
1584            0.65551583769728516
1585
1586        IMPLEMENTATION:
1587            Call the corresponding mwrank C++ library function.
1588        """
1589        if not self.is_minimal():
1590            raise RuntimeError, "curve must be minimal."
1591        return self.mwrank_curve().CPS_height_bound()
1592
1593
1594    def silverman_height_bound(self):
1595        r""" Return the Silverman height bound.  This is a positive real
1596        (floating point) number B such that for all rational points
1597        $P$ on the curve,
1598        $$
1599             h(P) \le \hat{h}(P) + B
1600        $$
1601        where h(P) is the logarithmic height of $P$ and $\hat{h}(P)$
1602        is the canonical height.
1603
1604        Note that the CPS_height_bound is often better (i.e. smaller)
1605        than the Silverman bound.
1606
1607        EXAMPLES:
1608            sage: E=EllipticCurve('37a1')
1609            sage: E.silverman_height_bound()
1610            4.8254007581809182
1611            sage: E.CPS_height_bound()
1612            0.16397076103046915
1613        """
1614        return self.mwrank_curve().silverman_bound()
1615       
1616
1617    def point_search(self, height_limit, verbose=True):
1618        """
1619        Search for points on a curve up to an input bound on the naive
1620        logarithmic height.
1621
1622        INPUT:
1623
1624            height_limit (float) -- bound on naive height (at most 21,
1625                                    or mwrank overflows -- see below)
1626
1627            verbose (bool) -- (default: True)
1628
1629                              If True, report on each point as found
1630                              together with linear relations between
1631                              the points found and the saturation
1632                              process.
1633
1634                              If False, just return the result.
1635
1636        OUTPUT: points (list) -- list of independent points which
1637            generate the subgroup of the Mordell-Weil group generated
1638            by the points found and then p-saturated for p<20.
1639
1640        WARNING: 
1641            height_limit is logarithmic, so increasing by 1 will cause
1642            the running time to increase by a factor of approximately
1643            4.5 (=exp(1.5)).  The limit of 21 is to prevent overflow,
1644            but in any case using height_limit=20 takes rather a long
1645            time!
1646
1647        IMPLEMENTATION: Uses Cremona's mwrank package. At the heart of
1648        this function is Cremona's port of Stoll's ratpoints program
1649        (version 1.4).
1650
1651        EXAMPLES:
1652            sage: E=EllipticCurve('389a1')
1653            sage: E.point_search(5, verbose=False)
1654            [(0 : -1 : 1), (-1 : 1 : 1)]
1655
1656        Increasing the height_limit takes longer, but finds no more points:
1657            sage: E.point_search(10, verbose=False)
1658            [(0 : -1 : 1), (-1 : 1 : 1)]
1659
1660        In fact this curve has rank 2 so no more than 2 points will
1661        ever be output, but we are not using this fact.
1662
1663            sage: E.saturation(_)
1664            ([(0 : -1 : 1), (-1 : 1 : 1)], '1', 0.152460172772408)
1665
1666        What this shows is that if the rank is 2 then the points
1667        listed do generate the Mordell-Weil group (mod torsion).
1668        Finally,
1669
1670            sage: E.rank()
1671            2
1672
1673        """
1674        height_limit = float(height_limit)
1675        if height_limit > _MAX_HEIGHT:
1676            raise OverflowError, "height_limit (=%s) must be at most %s."%(height_limit,_MAX_HEIGHT)
1677        c = self.mwrank_curve()
1678        mw = mwrank.mwrank_MordellWeil(c, verbose)
1679        mw.search(height_limit, verbose=verbose)
1680        v = mw.points()
1681        return [self(P) for P in v]
1682   
1683    def two_torsion_rank(self):
1684        r"""
1685        Return the dimension of the 2-torsion subgroup of $E(\Q)$.
1686
1687        This will be 0, 1 or 2.
1688
1689        NOTE: As s side-effect of calling this function, the full
1690        torsion subgroup of the curve is computed (if not already
1691        cached).  A simpler implementation of this function would be
1692        possible (by counting the roots of the 2-division polynomial),
1693        but the full torsion subgroup computation is not expensive.
1694       
1695        EXAMPLES:
1696            sage: EllipticCurve('11a1').two_torsion_rank()
1697            0
1698            sage: EllipticCurve('14a1').two_torsion_rank()
1699            1
1700            sage: EllipticCurve('15a1').two_torsion_rank()
1701            2
1702        """
1703        A = self.torsion_subgroup().invariants()
1704        if len(A) == 2:
1705            return rings.Integer(2)
1706        elif len(A) == 1 and A[0] % 2 == 0:
1707            return rings.Integer(1)
1708        else:
1709            return rings.Integer(0)
1710
1711    def selmer_rank_bound(self):
1712        """
1713        Bound on the rank of the curve, computed using the
1714        2-selmer group.  This is the rank of the curve
1715        minus the rank of the 2-torsion, minus a number
1716        determined by whatever mwrank was able to determine
1717        related to Sha[2].  Thus in many cases, this is
1718        the actual rank of the curve.
1719
1720        EXAMPLE:
1721        The following is the curve 960D1, which has rank 0,
1722        but Sha of order 4.
1723            sage: E = EllipticCurve([0, -1, 0, -900, -10098])
1724            sage: E.selmer_rank_bound()
1725            0
1726           
1727        It gives 0 instead of 2, because it knows Sha is nontrivial.
1728        In contrast, for the curve 571A, also with rank 0 and Sha
1729        of order 4, we get a worse bound:
1730            sage: E = EllipticCurve([0, -1, 1, -929, -10595])
1731            sage: E.selmer_rank_bound()
1732            2
1733            sage: E.rank(only_use_mwrank=False)   # uses L-function
1734            0
1735        """
1736        try:
1737            return self.__selmer_rank_bound
1738        except AttributeError:
1739            C = self.mwrank_curve()
1740            self.__selmer_rank_bound = C.selmer_rank_bound()
1741            return self.__selmer_rank_bound
1742   
1743           
1744    def an(self, n):
1745        """
1746        The n-th Fourier coefficient of the modular form corresponding
1747        to this elliptic curve, where n is a positive integer.
1748
1749        EXAMPLES:
1750            sage: E=EllipticCurve('37a1')
1751            sage: [E.an(n) for n in range(20) if n>0]
1752            [1, -2, -3, 2, -2, 6, -1, 0, 6, 4, -5, -6, -2, 2, 6, -4, 0, -12, 0]
1753        """
1754        return Integer(self.pari_mincurve().ellak(n))
1755
1756    def ap(self, p):
1757        """
1758        The p-th Fourier coefficient of the modular form corresponding
1759        to this elliptic curve, where p is prime.
1760
1761        EXAMPLES:
1762            sage: E=EllipticCurve('37a1')
1763            sage: [E.ap(p) for p in prime_range(50)]
1764            [-2, -3, -2, -1, -5, -2, 0, 0, 2, 6, -4, -1, -9, 2, -9]
1765
1766        """
1767        if not arith.is_prime(p):
1768            raise ArithmeticError, "p must be prime"
1769        return Integer(self.pari_mincurve().ellap(p))
1770
1771    def quadratic_twist(self, D):
1772        """
1773        Return the quadratic twist of this elliptic curve by D.
1774
1775        D must be a nonzero rational number.
1776
1777        NOTE: This function overrides the generic quadratic_twist()
1778            function for elliptic curves, returning a minimal model
1779
1780        EXAMPLES:
1781            sage: E=EllipticCurve('37a1')
1782            sage: E2=E.quadratic_twist(2); E2
1783            Elliptic Curve defined by y^2  = x^3 - 4*x + 2 over Rational Field
1784            sage: E2.conductor()
1785            2368
1786            sage: E2.quadratic_twist(2) == E
1787            True
1788       """
1789        return EllipticCurve_number_field.quadratic_twist(self, D).minimal_model()
1790   
1791    def minimal_model(self):
1792        r"""
1793        Return the unique minimal Weierstrass equation for this
1794        elliptic curve.  This is the model with minimal discriminant
1795        and $a_1,a_2,a_3 \in \{0,\pm 1\}$.
1796
1797        EXAMPLES:
1798            sage: E=EllipticCurve([10,100,1000,10000,1000000])
1799            sage: E.minimal_model()
1800            Elliptic Curve defined by y^2 + x*y + y = x^3 + x^2 + x +1 over Rational Field
1801        """
1802        try:
1803            return self.__minimal_model
1804        except AttributeError:
1805            F = self.pari_mincurve()
1806            self.__minimal_model = EllipticCurve_rational_field([Q(F[i]) for i in range(5)])
1807            return self.__minimal_model
1808
1809    def is_minimal(self):
1810        r"""
1811        Return True iff this elliptic curve is a reduced minimal model.
1812
1813        the unique minimal Weierstrass equation for this
1814        elliptic curve.  This is the model with minimal discriminant
1815        and $a_1,a_2,a_3 \in \{0,\pm 1\}$.
1816
1817        TO DO: This is not very efficient since it just computes the
1818        minimal model and compares.  A better implementation using the
1819        Kraus conditions would be preferable.
1820
1821        EXAMPLES:
1822            sage: E=EllipticCurve([10,100,1000,10000,1000000])
1823            sage: E.is_minimal()
1824            False
1825            sage: E=E.minimal_model()
1826            sage: E.is_minimal()
1827            True
1828        """
1829        return self.ainvs() == self.minimal_model().ainvs()
1830
1831#    Duplicate!
1832#
1833#    def is_integral(self):
1834#        for n in self.ainvs():
1835#            if n.denominator() != 1:
1836#                return False
1837#        return True
1838
1839    def kodaira_type(self, p):
1840        """
1841        Local Kodaira type of the elliptic curve at $p$.
1842
1843        INPUT:
1844           -- p, an integral prime
1845        OUTPUT:
1846           -- the kodaira type of this elliptic curve at p, as a KodairaSymbol.
1847
1848        EXAMPLES:
1849            sage: E = EllipticCurve('124a')
1850            sage: E.kodaira_type(2)
1851            IV
1852        """
1853        if not arith.is_prime(p):
1854            raise ArithmeticError, "p must be prime"
1855        try:
1856            self.__kodaira_type
1857        except AttributeError:
1858            self.__kodaira_type = {}
1859            self.__tamagawa_number = {}
1860        if not self.__kodaira_type.has_key(p):
1861            v = self.pari_mincurve().elllocalred(p)
1862            from kodaira_symbol import KodairaSymbol
1863            self.__kodaira_type[p] = KodairaSymbol(v[1])
1864            self.__tamagawa_number[p] = Integer(v[3])
1865        return self.__kodaira_type[p]
1866 
1867    def tamagawa_number(self, p):
1868        """
1869        The Tamagawa number of the elliptic curve at $p$.
1870
1871        EXAMPLES:
1872            sage: E = EllipticCurve('11a')
1873            sage: E.tamagawa_number(11)
1874            5
1875            sage: E = EllipticCurve('37b')
1876            sage: E.tamagawa_number(37)
1877            3
1878        """
1879        if not arith.is_prime(p):
1880            raise ArithmeticError, "p must be prime"
1881        try:
1882            return self.__tamagawa_number[p]
1883        except (AttributeError, KeyError):
1884            self.kodaira_type(p)
1885            return self.__tamagawa_number[p]
1886
1887    def tamagawa_numbers(self):
1888        """
1889        Return a list of all Tamagawa numbers for all prime divisors of
1890        the conductor (in order).
1891
1892        EXAMPLES:
1893            sage: e = EllipticCurve('30a1')
1894            sage: e.tamagawa_numbers()
1895            [2, 3, 1]
1896            sage: vector(e.tamagawa_numbers())
1897            (2, 3, 1)
1898        """
1899        return [self.tamagawa_number(p) for p in arith.prime_divisors(self.conductor())]
1900
1901    def tamagawa_product(self):
1902        """
1903        Returns the product of the Tamagawa numbers.
1904
1905        EXAMPLES:
1906            sage: E = EllipticCurve('54a')
1907            sage: E.tamagawa_product ()
1908            3
1909        """
1910        try:
1911            return self.__tamagawa_product
1912        except AttributeError:
1913            self.__tamagawa_product = self.pari_mincurve().ellglobalred()[2].python()
1914            return self.__tamagawa_product
1915
1916    def real_components(self):
1917        """
1918        Returns 1 if there is 1 real component and 2 if there are 2.
1919
1920        EXAMPLES:
1921            sage: E = EllipticCurve('37a')
1922            sage: E.real_components ()
1923            2
1924            sage: E = EllipticCurve('37b')
1925            sage: E.real_components ()
1926            2
1927            sage: E = EllipticCurve('11a')
1928            sage: E.real_components ()
1929            1
1930        """
1931        invs = self.short_weierstrass_model().ainvs()
1932        x = rings.polygen(self.base_ring())
1933        f = x**3 + invs[3]*x + invs[4]
1934        if f.discriminant() > 0:
1935            return 2
1936        else:
1937            return 1
1938
1939    def period_lattice(self):
1940        r"""
1941        Returns the period lattice of the elliptic curve.
1942
1943        EXAMPLES:
1944        sage: E = EllipticCurve('37a')
1945        sage: E.period_lattice()
1946        Period lattice associated to Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field
1947        """
1948        from sage.schemes.elliptic_curves.period_lattice import PeriodLattice_ell
1949        return PeriodLattice_ell(self)
1950
1951    def lseries(self):
1952        """
1953        Returns the L-series of this elliptic curve.
1954
1955        Further documentation is available for the functions which
1956        apply to the L-series.
1957
1958        EXAMPLES:
1959            sage: E=EllipticCurve('37a1')
1960            sage: E.lseries()
1961            Complex L-series of the Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field
1962        """
1963        try:
1964            return self.__lseries
1965        except AttributeError:
1966            from lseries_ell import Lseries_ell
1967            self.__lseries = Lseries_ell(self)
1968            return self.__lseries
1969
1970    def Lambda(self, s, prec):
1971        r"""
1972        Returns the value of the Lambda-series of the elliptic curve E
1973        at s, where s can be any complex number.
1974
1975        IMPLEMENTATION: Fairly *slow* computation using the definitions
1976        and implemented in Python.
1977
1978        Uses prec terms of the power series.
1979
1980        EXAMPLES:
1981            sage: E = EllipticCurve('389a')
1982            sage: E.Lambda(1.4+0.5*I, 50)
1983            -0.354172680517... + 0.874518681720...*I
1984        """
1985        s = C(s)
1986        N = self.conductor()
1987        pi = R(constants.pi)
1988        Gamma = transcendental.gamma
1989        Gamma_inc = transcendental.gamma_inc
1990        a = self.anlist(prec)
1991        eps = self.root_number()
1992        sqrtN = float(N.sqrt())
1993        def _F(n, t):
1994            return Gamma_inc(t+1, 2*pi*n/sqrtN) * C(sqrtN/(2*pi*n))**(t+1)
1995        return sum([a[n]*(_F(n,s-1) + eps*_F(n,1-s)) for n in xrange(1,prec+1)])
1996               
1997    def is_local_integral_model(self,*p):
1998        r"""
1999        Tests if self is integral at the prime $p$, or at all the
2000        primes if $p$ is a list or tuple of primes
2001       
2002        EXAMPLES:
2003            sage: E=EllipticCurve([1/2,1/5,1/5,1/5,1/5])
2004            sage: [E.is_local_integral_model(p) for p in (2,3,5)]
2005            [False, True, False]
2006            sage: E.is_local_integral_model(2,3,5)
2007            False
2008            sage: Eint2=E.local_integral_model(2)
2009            sage: Eint2.is_local_integral_model(2)
2010            True
2011        """
2012        if len(p)==1: p=p[0]
2013        if isinstance(p,(tuple,list)):
2014            return misc.forall(p, lambda x : self.is_local_integral_model(x))[0]
2015        assert p.is_prime(), "p must be prime in is_local_integral_model()"
2016        return misc.forall(self.ainvs(), lambda x : x.valuation(p) >= 0)[0]
2017
2018    def local_integral_model(self,p):
2019        r"""
2020        Return a model of self which is integral at the prime $p$
2021       
2022        EXAMPLES:
2023            sage: E=EllipticCurve([0, 0, 1/216, -7/1296, 1/7776])
2024            sage: E.local_integral_model(2)
2025             Elliptic Curve defined by y^2 + 1/27*y = x^3 - 7/81*x + 2/243 over Rational Field
2026            sage: E.local_integral_model(3)
2027             Elliptic Curve defined by y^2 + 1/8*y = x^3 - 7/16*x + 3/32 over Rational Field
2028            sage: E.local_integral_model(2).local_integral_model(3) == EllipticCurve('5077a1')
2029            True
2030        """
2031        assert p.is_prime(), "p must be prime in local_integral_model()"
2032        ai = self.a_invariants()
2033        e  = min([(ai[i].valuation(p)/[1,2,3,4,6][i]) for i in range(5)]).floor()
2034        return constructor.EllipticCurve([ai[i]/p**(e*[1,2,3,4,6][i]) for i in range(5)])
2035
2036    def is_global_integral_model(self):
2037        r"""
2038        Return true iff self is integral at all primes
2039       
2040        EXAMPLES:
2041            sage: E=EllipticCurve([1/2,1/5,1/5,1/5,1/5])
2042            sage: E.is_global_integral_model()
2043            False
2044            sage: Emin=E.global_integral_model()
2045            sage: Emin.is_global_integral_model()
2046            True
2047        """
2048        return self.is_integral()
2049
2050    def global_integral_model(self):
2051        r"""
2052        Return a model of self which is integral at all primes
2053       
2054        EXAMPLES:
2055            sage: E = EllipticCurve([0, 0, 1/216, -7/1296, 1/7776])
2056            sage: F = E.global_integral_model(); F
2057            Elliptic Curve defined by y^2 + y = x^3 - 7*x + 6 over Rational Field
2058            sage: F == EllipticCurve('5077a1')
2059            True
2060        """
2061        ai = self.a_invariants()
2062        for a in ai:
2063            if not a.is_integral():
2064               for p, _ in a.denom().factor():
2065                  e  = min([(ai[i].valuation(p)/[1,2,3,4,6][i]) for i in range(5)]).floor()
2066                  ai = [ai[i]/p**(e*[1,2,3,4,6][i]) for i in range(5)]
2067        for z in ai:
2068            assert z.denominator() == 1, "bug in global_integral_model: %s" % ai
2069        return constructor.EllipticCurve(ai)
2070
2071    def integral_model(self):
2072        r"""
2073        Return a weierstrass model, $F$, of self with integral coefficients,
2074        along with a morphism $\phi$ of points on self to points on $F$.
2075       
2076        EXAMPLES:
2077            sage: E = EllipticCurve([1/2,0,0,5,1/3])
2078            sage: F, phi = E.integral_model()
2079            sage: F
2080            Elliptic Curve defined by y^2 + 3*x*y  = x^3 + 6480*x + 15552 over Rational Field
2081            sage: phi
2082            Generic morphism:
2083              From: Abelian group of points on Elliptic Curve defined by y^2 + 1/2*x*y  = x^3 + 5*x + 1/3 over Rational Field
2084              To:   Abelian group of points on Elliptic Curve defined by y^2 + 3*x*y  = x^3 + 6480*x + 15552 over Rational Field
2085              Via:  (u,r,s,t) = (1/6, 0, 0, 0)
2086            sage: P = E([4/9,41/27])
2087            sage: phi(P)
2088            (16 : 328 : 1)
2089            sage: phi(P) in F
2090            True
2091        """
2092        F = self.global_integral_model()
2093        return F, self.isomorphism_to(F)
2094
2095    def integral_weierstrass_model(self):
2096        r"""
2097        Return a model of the form $y^2 = x^3 + a*x + b$ for this curve with $a,b\in\Z$.
2098       
2099        EXAMPLES:
2100            sage: E = EllipticCurve('17a1')
2101            sage: E.integral_weierstrass_model()
2102            Elliptic Curve defined by y^2  = x^3 - 11*x - 890 over Rational Field
2103        """
2104        F = self.minimal_model().short_weierstrass_model()
2105        _,_,_,A,B = F.ainvs()
2106        for p in [2,3]:
2107            e=min(A.valuation(p)/4,B.valuation(p)/6).floor()
2108            A /= Integer(p**(4*e))
2109            B /= Integer(p**(6*e))
2110        return constructor.EllipticCurve([A,B])
2111
2112    def modular_degree(self, algorithm='sympow'):
2113        r"""
2114        Return the modular degree of this elliptic curve.
2115
2116        The result is cached.  Subsequence calls, even with a
2117        different algorithm, just returned the cached result.
2118
2119        INPUT:
2120           algorithm -- string:
2121              'sympow' -- (default) use Mark Watkin's (newer) C program sympow
2122              'magma' -- requires that MAGMA be installed (also implemented
2123                         by Mark Watkins)
2124
2125        \note{On 64-bit computers ec does not work, so \sage uses
2126        sympow even if ec is selected on a 64-bit computer.}
2127
2128        The correctness of this function when called with algorithm "ec"
2129        is subject to the following three hypothesis:
2130
2131        \begin{itemize}
2132
2133            \item Manin's conjecture: the Manin constant is 1
2134           
2135            \item Steven's conjecture: the $X_1(N)$-optimal quotient
2136            is the curve with minimal Faltings height.  (This is proved
2137            in most cases.)
2138
2139            \item The modular degree fits in a machine double, so it
2140            better be less than about 50-some bits.  (If you use sympow
2141            this contraint does not apply.)
2142
2143        \end{itemize}
2144
2145        Moreover for all algorithms, computing a certain value of an
2146        $L$-function ``uses a heuristic method that discerns when the
2147        real-number approximation to the modular degree is within epsilon
2148        [=0.01 for algorithm="sympow"] of the same integer for 3
2149        consecutive trials (which occur maybe every 25000 coefficients
2150        or so). Probably it could just round at some point. For
2151        rigour, you would need to bound the tail by assuming
2152        (essentially) that all the $a_n$ are as large as possible, but
2153        in practise they exhibit significant (square root)
2154        cancellation. One difficulty is that it doesn't do the sum in
2155        1-2-3-4 order; it uses 1-2-4-8---3-6-12-24--9-18-- (Euler
2156        product style) instead, and so you have to guess ahead of time
2157        at what point to curtail this expansion.''  (Quote from an
2158        email of Mark Watkins.)
2159
2160        \note{If the curve is loaded from the large Cremona database,
2161        then the modular degree is taken from the database.}
2162
2163        EXAMPLES:
2164            sage: E = EllipticCurve([0, -1, 1, -10, -20])
2165            sage: E
2166            Elliptic Curve defined by y^2 + y = x^3 - x^2 - 10*x - 20 over Rational Field
2167            sage: E.modular_degree()
2168            1                                     
2169            sage: E = EllipticCurve('5077a')
2170            sage: E.modular_degree()
2171            1984                                   
2172            sage: factor(1984)
2173            2^6 * 31
2174
2175            sage: EllipticCurve([0, 0, 1, -7, 6]).modular_degree()
2176            1984
2177            sage: EllipticCurve([0, 0, 1, -7, 6]).modular_degree(algorithm='sympow')
2178            1984
2179            sage: EllipticCurve([0, 0, 1, -7, 6]).modular_degree(algorithm='magma')  # optional
2180            1984
2181
2182        We compute the modular degree of the curve with rank 4 having smallest
2183        (known) conductor:
2184       
2185            sage: E = EllipticCurve([1, -1, 0, -79, 289])
2186            sage: factor(E.conductor())  # conductor is 234446
2187            2 * 117223
2188            sage: factor(E.modular_degree())
2189            2^7 * 2617
2190        """
2191        try:
2192            return self.__modular_degree
2193       
2194        except AttributeError:
2195            if algorithm == 'sympow':
2196                from sage.lfunctions.all import sympow
2197                m = sympow.modular_degree(self)
2198            elif algorithm == 'magma':
2199                m = rings.Integer(self._magma_().ModularDegree())
2200            else:
2201                raise ValueError, "unknown algorithm %s"%algorithm
2202            self.__modular_degree = m
2203            return m
2204
2205    def modular_parametrization(self):
2206        r"""
2207        Computes and returns the modular parametrization of this elliptic curve.
2208
2209        The curve is converted to a minimal model.
2210
2211        OUTPUT: A list of two Larent series [X(x),Y(x)] of degrees -2,
2212        -3 respectively, which satisfy the equation of the (minimal
2213        model of the) elliptic curve.  The are modular functions on
2214        $\Gamma_0(N)$ where $N$ is the conductor.
2215
2216        X.deriv()/(2*Y+a1*X+a3) should equal f(q)dq/q where f is
2217        self.q_expansion().
2218
2219        EXAMPLES:
2220            sage: E=EllipticCurve('389a1')
2221            sage: X,Y=E.modular_parametrization()
2222            sage: X
2223            q^-2 + 2*q^-1 + 4 + 7*q + 13*q^2 + 18*q^3 + 31*q^4 + 49*q^5 + 74*q^6 + 111*q^7 + 173*q^8 + 251*q^9 + 379*q^10 + 560*q^11 + 824*q^12 + 1199*q^13 + 1773*q^14 + 2365*q^15 + 3463*q^16 + 4508*q^17 + O(q^18)
2224            sage: Y
2225            -q^-3 - 3*q^-2 - 8*q^-1 - 17 - 33*q - 61*q^2 - 110*q^3 - 186*q^4 - 320*q^5 - 528*q^6 - 861*q^7 - 1383*q^8 - 2218*q^9 - 3472*q^10 - 5451*q^11 - 8447*q^12 - 13020*q^13 - 20083*q^14 - 29512*q^15 - 39682*q^16 + O(q^17)
2226
2227            The following should give 0, but only approximately:
2228            sage: q = X.parent().gen()
2229            sage: E.defining_polynomial()(X,Y,1) + O(q^11) == 0
2230            True
2231
2232            Note that below we have to change variable from x to q
2233            sage: a1,_,a3,_,_=E.a_invariants()
2234            sage: f=E.q_expansion(17)
2235            sage: q=f.parent().gen()
2236            sage: f/q == (X.derivative()/(2*Y+a1*X+a3))
2237            True
2238
2239        """
2240        R = LaurentSeriesRing(RationalField(),'q')
2241        XY = self.pari_mincurve().elltaniyama()
2242        return [1/R(1/XY[0]),1/R(1/XY[1])]
2243
2244    def cremona_label(self, space=False):
2245        """
2246        Return the Cremona label associated to (the minimal model) of
2247        this curve, if it is known.  If not, raise a RuntimeError
2248        exception.
2249
2250        EXAMPLES:
2251            sage: E=EllipticCurve('389a1')
2252            sage: E.cremona_label()
2253            '389a1'
2254
2255        The default database only contains conductors up to 10000, so
2256        any curve with conductor greater than that will cause an error
2257        to be raised:
2258
2259            sage: E=EllipticCurve([1,2,3,4,5]);
2260            sage: E.conductor()
2261            10351
2262            sage: E.cremona_label()
2263            Traceback (most recent call last):
2264            ...
2265            RuntimeError: Cremona label not known for Elliptic Curve defined by y^2 + x*y + 3*y = x^3 + 2*x^2 + 4*x + 5 over Rational Field.
2266        """
2267        try:
2268            if not space:
2269                return self.__cremona_label.replace(' ','')
2270            return self.__cremona_label
2271        except AttributeError:
2272            try:
2273                X = self.database_curve()
2274            except RuntimeError:
2275                raise RuntimeError, "Cremona label not known for %s."%self
2276            self.__cremona_label = X.__cremona_label
2277            return self.cremona_label(space)
2278
2279    label = cremona_label
2280
2281    def torsion_order(self):
2282        """
2283        Return the order of the torsion subgroup.
2284
2285        EXAMPLES:
2286            sage: e = EllipticCurve('11a')
2287            sage: e.torsion_order()
2288            5
2289            sage: type(e.torsion_order())
2290            <type 'sage.rings.integer.Integer'>
2291            sage: e = EllipticCurve([1,2,3,4,5])
2292            sage: e.torsion_order()
2293            1
2294            sage: type(e.torsion_order())
2295            <type 'sage.rings.integer.Integer'>
2296        """
2297        try:
2298            return self.__torsion_order
2299        except AttributeError:
2300            self.__torsion_order = self.torsion_subgroup().order()
2301            return self.__torsion_order
2302
2303    def torsion_subgroup(self, flag=0):
2304        """
2305        Returns the torsion subgroup of this elliptic curve.
2306
2307        INPUT:
2308            flag -- (default: 0)  chooses PARI algorithm:
2309              flag = 0: uses Doud algorithm
2310              flag = 1: uses Lutz-Nagell algorithm
2311
2312        OUTPUT:
2313            The EllipticCurveTorsionSubgroup instance associated to this elliptic curve.
2314
2315        EXAMPLES:
2316            sage: EllipticCurve('11a').torsion_subgroup()
2317            Torsion Subgroup isomorphic to Multiplicative Abelian Group isomorphic to C5 associated to the Elliptic Curve defined by y^2 + y = x^3 - x^2 - 10*x - 20 over Rational Field
2318            sage: EllipticCurve('37b').torsion_subgroup()
2319            Torsion Subgroup isomorphic to Multiplicative Abelian Group isomorphic to C3 associated to the Elliptic Curve defined by y^2 + y = x^3 + x^2 - 23*x - 50 over Rational Field
2320
2321            sage: e = EllipticCurve([-1386747,368636886]);e
2322            Elliptic Curve defined by y^2  = x^3 - 1386747*x + 368636886 over Rational Field
2323            sage: G = e.torsion_subgroup(); G
2324            Torsion Subgroup isomorphic to Multiplicative Abelian
2325            Group isomorphic to C8 x C2 associated to the Elliptic
2326            Curve defined by y^2 = x^3 - 1386747*x + 368636886 over
2327            Rational Field
2328            sage: G.0
2329            (1227 : 22680 : 1)
2330            sage: G.1
2331            (282 : 0 : 1)
2332            sage: list(G)
2333            [1, P1, P0, P0*P1, P0^2, P0^2*P1, P0^3, P0^3*P1, P0^4, P0^4*P1, P0^5, P0^5*P1, P0^6, P0^6*P1, P0^7, P0^7*P1]       
2334        """
2335        try:
2336            return self.__torsion_subgroup
2337        except AttributeError:
2338            self.__torsion_subgroup = rational_torsion.EllipticCurveTorsionSubgroup(self, flag)
2339            self.__torsion_order = self.__torsion_subgroup.order()
2340            return self.__torsion_subgroup
2341
2342    ## def newform_eval(self, z, prec):
2343##         """
2344##         The value of the newform attached to this elliptic curve at
2345##         the point z in the complex upper half plane, computed using
2346##         prec terms of the power series expansion.  Note that the power
2347##         series need not converge well near the real axis.
2348##         """
2349##         raise NotImplementedError
2350
2351    def root_number(self):
2352        """
2353        Returns the root number of this elliptic curve.
2354
2355        This is 1 if the order of vanishing of the L-function L(E,s)
2356        at 1 is even, and -1 if it is odd.
2357
2358        EXAMPLES:
2359            sage: EllipticCurve('11a1').root_number()
2360            1
2361            sage: EllipticCurve('37a1').root_number()
2362            -1
2363            sage: EllipticCurve('389a1').root_number()
2364            1
2365        """
2366        try:
2367            return self.__root_number
2368        except AttributeError:
2369            self.__root_number = int(self.pari_mincurve().ellrootno())
2370        return self.__root_number
2371
2372
2373    def has_cm(self):
2374        """
2375        Returns True iff this elliptic curve has Complex Multiplication
2376
2377        EXAMPLES:
2378            sage: E=EllipticCurve('37a1')
2379            sage: E.has_cm()
2380            False
2381            sage: E=EllipticCurve('32a1')
2382            sage: E.has_cm()
2383            True
2384            sage: E.j_invariant()
2385            1728
2386        """
2387
2388        return CMJ.has_key(self.j_invariant())
2389
2390    def cm_discriminant(self):
2391        """
2392        Returns the associated quadratic discriminant if this elliptic
2393        curve has Complex Multiplication.
2394       
2395        A ValueError is raised if the curve does not have CM (see the
2396        function has_cm())
2397
2398        EXAMPLES:
2399            sage: E=EllipticCurve('32a1')
2400            sage: E.cm_discriminant()
2401            -4
2402            sage: E=EllipticCurve('121b1')
2403            sage: E.cm_discriminant()
2404            -11
2405            sage: E=EllipticCurve('37a1')
2406           sage: E.cm_discriminant()
2407           Traceback (most recent call last):
2408           ...
2409           ValueError: Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field does not have CM
2410        """
2411
2412        try:
2413            return CMJ[self.j_invariant()]
2414        except KeyError:
2415            raise ValueError, "%s does not have CM"%self
2416           
2417
2418    def quadratic_twist(self, D):
2419        """
2420        Return the global minimal model of the quadratic twist of this
2421        curve by D.
2422
2423        EXAMPLES:
2424            sage: E=EllipticCurve('37a1')
2425            sage: E7=E.quadratic_twist(7); E7
2426            Elliptic Curve defined by y^2  = x^3 - 784*x + 5488 over Rational Field
2427            sage: E7.conductor()
2428            29008
2429            sage: E7.quadratic_twist(7) == E
2430            True
2431        """
2432        return EllipticCurve_number_field.quadratic_twist(self, D).minimal_model()
2433       
2434
2435    ##########################################################
2436    # Isogeny class (currently just uses Cremona database.)
2437    ##########################################################
2438    def isogeny_class(self, algorithm="mwrank", verbose=False):
2439        r"""
2440        Return all curves over $\Q$ in the isogeny class of this
2441        elliptic curve.
2442
2443        INPUT:
2444            algorithm -- string:
2445                 "mwrank"   -- (default) use the mwrank C++ library
2446                 "database" -- use the Cremona database (only works if
2447                               curve is isomorphic to a curve in the database)
2448
2449        OUTPUT:
2450            Returns the sorted list of the curves isogenous to self.
2451            If algorithm is "mwrank", also returns the isogeny matrix (otherwise
2452            returns None as second return value). 
2453
2454        \note{The result is \emph{not} provably correct, in the sense
2455            that when the numbers are huge isogenies could be missed
2456            because of precision issues.}
2457
2458        \note{The ordering depends on which algorithm is used.}
2459
2460        EXAMPLES:
2461            sage: I, A = EllipticCurve('37b').isogeny_class('mwrank') 
2462            sage: I   # randomly ordered
2463            [Elliptic Curve defined by y^2 + y = x^3 + x^2 - 23*x - 50 over Rational Field,
2464             Elliptic Curve defined by y^2 + y = x^3 + x^2 - 1873*x - 31833 over Rational Field,
2465             Elliptic Curve defined by y^2 + y = x^3 + x^2 - 3*x +1 over Rational Field]
2466            sage: A
2467            [0 3 3]
2468            [3 0 0]
2469            [3 0 0]
2470
2471            sage: I, _ = EllipticCurve('37b').isogeny_class('database'); I
2472            [Elliptic Curve defined by y^2 + y = x^3 + x^2 - 1873*x - 31833 over Rational Field,
2473             Elliptic Curve defined by y^2 + y = x^3 + x^2 - 23*x - 50 over Rational Field,
2474             Elliptic Curve defined by y^2 + y = x^3 + x^2 - 3*x +1 over Rational Field]
2475
2476        This is an example of a curve with a $37$-isogeny:
2477            sage: E = EllipticCurve([1,1,1,-8,6])
2478            sage: E.isogeny_class ()
2479            ([Elliptic Curve defined by y^2 + x*y + y = x^3 + x^2 - 8*x + 6 over Rational Field,
2480              Elliptic Curve defined by y^2 + x*y + y = x^3 + x^2 - 208083*x - 36621194 over Rational Field],
2481             [ 0 37]
2482             [37  0])
2483
2484        This curve had numerous $2$-isogenies:
2485        sage: e=EllipticCurve([1,0,0,-39,90])
2486            sage: e.isogeny_class ()
2487            ([Elliptic Curve defined by y^2 + x*y  = x^3 - 39*x + 90 over Rational Field,
2488              Elliptic Curve defined by y^2 + x*y  = x^3 - 4*x -1 over Rational Field,
2489              Elliptic Curve defined by y^2 + x*y  = x^3 + x over Rational Field,
2490              Elliptic Curve defined by y^2 + x*y  = x^3 - 49*x - 136 over Rational Field,
2491              Elliptic Curve defined by y^2 + x*y  = x^3 - 34*x - 217 over Rational Field,
2492              Elliptic Curve defined by y^2 + x*y  = x^3 - 784*x - 8515 over Rational Field],
2493             [0 2 0 0 0 0]
2494             [2 0 2 2 0 0]
2495             [0 2 0 0 0 0]
2496             [0 2 0 0 2 2]
2497             [0 0 0 2 0 0]
2498             [0 0 0 2 0 0])
2499
2500        See \url{http://modular.ucsd.edu/Tables/nature/} for more interesting
2501        examples of isogeny structures.
2502        """
2503        #if algorithm == "gp":
2504           
2505       #     return sum([L for _, L in self.isogenous_curves(algorithm="gp")], [self])
2506       
2507        if algorithm == "mwrank":
2508            try:
2509                E = self.mwrank_curve()
2510            except ValueError:
2511                E = self.minimal_model().mwrank_curve()
2512            I, A = E.isogeny_class(verbose=verbose)
2513            mat = matrix.MatrixSpace(rings.IntegerRing(), len(A))(A)
2514            I = [constructor.EllipticCurve(ainvs) for ainvs in I]
2515            return I, mat
2516       
2517        elif algorithm == "database":
2518           
2519            try:
2520                label = self.cremona_label(space=False)
2521            except RuntimeError:
2522                raise RuntimeError, "unable to to find %s in the database"%self
2523            db = sage.databases.cremona.CremonaDatabase()
2524            I = db.isogeny_class(label)
2525            I.sort()
2526            return I, None
2527       
2528        else:
2529
2530            raise ValueError, "unknown algorithm '%s'%"%algorithm
2531   
2532    def isogeny_graph(self):
2533        r"""
2534        Returns a graph representing the isogeny class of this elliptic curve,
2535        where the vertices are isogenous curves over $\Q$ and the edges are
2536        prime degree isogenies labeled by their degree.
2537
2538        EXAMPLES:
2539            sage: LL = []
2540            sage: for e in cremona_optimal_curves(range(1, 38)):
2541            ...    G = e.isogeny_graph()
2542            ...    already = False
2543            ...    for H in LL:
2544            ...        if G.is_isomorphic(H):
2545            ...            already = True
2546            ...            break
2547            ...    if not already:
2548            ...        LL.append(G)
2549            ...
2550            sage: graphs_list.show_graphs(LL)
2551           
2552            sage: E = EllipticCurve('195a')
2553            sage: G = E.isogeny_graph()
2554            sage: for v in G: print v, G.get_vertex(v)
2555            ...
2556            0 Elliptic Curve defined by y^2 + x*y  = x^3 - 110*x + 435 over Rational Field
2557            1 Elliptic Curve defined by y^2 + x*y  = x^3 - 115*x + 392 over Rational Field
2558            2 Elliptic Curve defined by y^2 + x*y  = x^3 + 210*x + 2277 over Rational Field
2559            3 Elliptic Curve defined by y^2 + x*y  = x^3 - 520*x - 4225 over Rational Field
2560            4 Elliptic Curve defined by y^2 + x*y  = x^3 + 605*x - 19750 over Rational Field
2561            5 Elliptic Curve defined by y^2 + x*y  = x^3 - 8125*x - 282568 over Rational Field
2562            6 Elliptic Curve defined by y^2 + x*y  = x^3 - 7930*x - 296725 over Rational Field
2563            7 Elliptic Curve defined by y^2 + x*y  = x^3 - 130000*x - 18051943 over Rational Field
2564            sage: G.plot(edge_labels=True).save('isogeny_graph.png')
2565        """
2566        from sage.graphs.graph import Graph
2567        L, M = self.isogeny_class()
2568        G = Graph(M, format='weighted_adjacency_matrix')
2569        d = {}
2570        for v in G.vertices():
2571            d[v] = L[v]
2572        G.set_vertices(d)
2573        return G
2574
2575    ##########################################################
2576    # Galois Representations
2577    ##########################################################
2578
2579    def is_reducible(self, p):
2580        """
2581        Return True if the mod-p representation attached
2582        to E is reducible.
2583
2584        INPUT:
2585            p -- a prime number
2586
2587        NOTE: The answer is cached.
2588
2589        EXAMPLES:
2590            sage: E = EllipticCurve('121a'); E
2591            Elliptic Curve defined by y^2 + x*y + y = x^3 + x^2 - 30*x - 76 over Rational Field
2592            sage: E.is_reducible(7)
2593            False
2594            sage: E.is_reducible(11)
2595            True
2596            sage: EllipticCurve('11a').is_reducible(5)
2597            True
2598            sage: e = EllipticCurve('11a2')
2599            sage: e.is_reducible(5)
2600            True
2601            sage: e.torsion_order()
2602            1       
2603        """
2604        try:
2605            return self.__is_reducible[p]
2606        except AttributeError:
2607            self.__is_reducible = {}
2608        except KeyError:
2609            pass
2610
2611        if not arith.is_prime(p):
2612            raise ValueError, 'p (=%s) must be prime'%p
2613        # we do is_surjective first, since this is
2614        # much easier than computing isogeny_class
2615        t, why = self.is_surjective(p) 
2616        if t == True:
2617            self.__is_reducible[p] = False
2618            return False  # definitely not reducible
2619        isogeny_matrix = self.isogeny_class()[ 1 ]
2620        v = isogeny_matrix.row(0) # first row
2621        for a in v:
2622            if a != 0 and a % p == 0:
2623                self.__is_reducible[p] = True
2624                return True
2625        self.__is_reducible[p] = False
2626        return False
2627
2628    def is_irreducible(self, p):
2629        """
2630        Return True if the mod p represenation is irreducible.
2631
2632        EXAMPLES:
2633            sage: e = EllipticCurve('37b')
2634            sage: e.is_irreducible(2)
2635            True
2636            sage: e.is_irreducible(3)
2637            False
2638            sage: e.is_reducible(2)
2639            False
2640            sage: e.is_reducible(3)
2641            True       
2642        """
2643        return not self.is_reducible(p)
2644       
2645    def is_surjective(self, p, A=1000):
2646        """
2647        Return True if the mod-p representation attached to E
2648        is surjective, False if it is not, or None if we were
2649        unable to determine whether it is or not.
2650
2651        NOTE: The answer is cached.
2652
2653        INPUT:
2654            p -- int (a prime number)
2655            A -- int (a bound on the number of a_p to use)
2656
2657        OUTPUT:
2658            a 2-tuple:
2659            -- surjective or (probably) not
2660            -- information about what it is if not surjective
2661
2662        EXAMPLES:
2663            sage: e = EllipticCurve('37b')
2664            sage: e.is_surjective(2)
2665            (True, None)
2666            sage: e.is_surjective(3)
2667            (False, '3-torsion')
2668           
2669
2670        REMARKS:
2671       
2672            1.  If p >= 5 then the mod-p representation is surjective
2673                if and only if the p-adic representation is
2674                surjective.  When p = 2, 3 there are counterexamples.
2675                See a very recent paper of Elkies for more details
2676                when p=3.
2677
2678            2.  When p <= 3 this function always gives the correct
2679                result irregardless of A, since it explicitly
2680                determines the p-division polynomial.
2681           
2682        """
2683        if not arith.is_prime(p):
2684            raise TypeError, "p (=%s) must be prime."%p
2685        A = int(A)
2686        key = (p, A)
2687        try:
2688            return self.__is_surjective[key]
2689        except KeyError:
2690            pass
2691        except AttributeError:
2692            self.__is_surjective = {}
2693
2694        ans = self._is_surjective(p, A)
2695        self.__is_surjective[key] = ans
2696        return ans
2697
2698    def _is_surjective(self, p, A):
2699        T = self.torsion_subgroup().order()
2700        if T % p == 0:
2701            return False, "%s-torsion"%p
2702
2703        if p == 2:
2704            # E is isomorphic to  [0,b2,0,8*b4,16*b6]
2705            b2,b4,b6,b8=self.b_invariants()
2706            R = rings.PolynomialRing(self.base_ring(), 'x')
2707            x = R.gen()
2708            f = x**3 + b2*x**2 + 8*b4*x + 16*b6
2709            if not f.is_irreducible():
2710                return False, '2-torsion'
2711            if arith.is_square(f.discriminant()):
2712                return False, "A3"
2713            return True, None
2714       
2715        if p == 3:
2716            # Algorithm: Let f be the 3-division polynomial, which is
2717            # a polynomial of degree 4.  Then I claim that this
2718            # polynomial has Galois group S_4 if and only if the
2719            # representation rhobar_{E,3} is surjective.  If the group
2720            # is S_4, then S_4 is a quotient of the image of
2721            # rhobar_{E,3}.  Since S_4 has order 24 and GL_2(F_3)
2722            # has order 48, the only possibility we have to consider
2723            # is that the image of rhobar is isomorphic to S_4.
2724            # But this is not the case because S_4 is not a subgroup
2725            # of GL_2(F_3).    If it were, it would be normal, since
2726            # it would have index 2.  But there is a *unique* normal
2727            # subgroup of GL_2(F_3) of index 2, namely SL_2(F_3),
2728            # and SL_2(F_3) is not isomorphic to S_4 (S_4 has a normal
2729            # subgroup of index 2 and SL_2(F_3) does not.)
2730            # (What's a simple way to see that SL_2(F_3) is the
2731            # unique index-2 normal subgroup?  I didn't see an obvious
2732            # reason, so just used the NormalSubgroups command in MAGMA
2733            # and it output exactly one of index 2.)
2734
2735            # Here's Noam Elkies proof for the other direction:
2736           
2737            #> Let E be an elliptic curve over Q.  Is the mod-3
2738            #> representation E[3]  surjective if and only if the
2739            #> (degree 4) division polynomial has Galois group S_4?  I
2740            #> can see why the group being S_4 implies the
2741            #> representation is surjective, but the converse is not
2742            #> clear to me.
2743            # I would have thought that this is the easier part: to
2744            # say that E[3] is surjective is to say the 3-torsion
2745            # field Q(E[3]) has Galois group GL_2(Z/3) over Q.  Let
2746            # E[3]+ be the subfield fixed by the element -1 of
2747            # GL_2(Z/3).  Then E[3] has Galois group PGL_2(Z/3), which
2748            # is identified with S_4 by its action on the four
2749            # 3-element subgroups of E[3].  Each such subgroup is in
2750            # turn determined by the x-coordinate shared by its two
2751            # nonzero points.  So, if E[3] is surjective then any
2752            # permutation of those x-coordinates is realized by some
2753            # element of Gal(E[3]+/Q).  Thus the Galois group of the
2754            # division polynomial (whose roots are those
2755            # x-coordinates) maps surjectively to S_4, which means it
2756            # equals S_4.
2757           
2758           
2759            f = self.division_polynomial(3)
2760            if not f.is_irreducible():
2761                return False, "reducible_3-divpoly"
2762            n = pari(f).polgalois()[0]
2763            if n == 24:
2764                return True, None
2765            else:
2766                return False, "3-divpoly_galgroup_order_%s"%n
2767       
2768        if self.has_cm():
2769            return False, "CM"
2770        an = self.anlist(A)
2771        ell = 0
2772        Np = self.conductor() * p
2773        signs = []
2774        while True:
2775            ell = arith.next_prime(ell)
2776            if ell >= A: break
2777            if Np % ell != 0:
2778                a_ell = an[int(ell)]
2779                if a_ell % p != 0:
2780                    s = arith.kronecker(a_ell**2 - 4*ell, p)
2781                    #print ell, s
2782                    if s == 0: continue
2783                    if not (s in signs):
2784                        signs.append(s)
2785                        if len(signs) == 2:
2786                            return True, None
2787           
2788        # could do something further here...
2789        return False, signs
2790
2791    def is_semistable(self):
2792        """
2793        Return True iff this elliptic curve is semi-stable at all primes.
2794
2795        EXAMPLES:
2796            sage: E=EllipticCurve('37a1')
2797            sage: E.is_semistable()
2798            True
2799            sage: E=EllipticCurve('90a1')
2800            sage: E.is_semistable()
2801            False
2802        """
2803        if self.base_ring() != Q:
2804            raise NotImplementedError, "is_semistable only implemented for curves over the rational numbers."
2805        return self.conductor().is_squarefree()
2806
2807    def reducible_primes(self):
2808        r"""
2809        Returns a list of the primes $p$ such that the mod $p$
2810        representation $\rho_{E,p}$ is reducible.  For all other
2811        primes the representation is irreducible.
2812
2813        NOTE -- this is \emph{not} provably correct in general.
2814        See the documentation for \code{self.isogeny_class}.
2815
2816        EXAMPLES:
2817            sage: E = EllipticCurve('225a')
2818            sage: E.reducible_primes()
2819            [3]
2820        """
2821        try:
2822            return self.__reducible_primes
2823        except AttributeError:
2824            pass
2825        C, I = self.isogeny_class(algorithm='mwrank')
2826        X = set(I.list())
2827        R = [p for p in X if arith.is_prime(p)]
2828        self.__reducible_primes = R
2829        return R
2830
2831    def non_surjective(self, A=1000):
2832        r"""
2833        Returns a list of primes p such that the mod-p representation
2834        $\rho_{E,p}$ *might* not be surjective (this list usually
2835        contains 2, because of shortcomings of the algorithm).  If p
2836        is not in the returned list, then rho_{E,p} is provably
2837        surjective (see A. Cojocaru's paper).  If the curve has CM
2838        then infinitely many representations are not surjective, so we
2839        simply return the sequence [(0,"cm")] and do no further computation.
2840
2841        INPUT:
2842            A -- an integer
2843        OUTPUT:
2844            list -- if curve has CM, returns [(0,"cm")].  Otherwise, returns a
2845                    list of primes where mod-p representation very likely
2846                    not surjective.   At any prime not in this list,
2847                    the representation is definitely surjective.
2848        EXAMPLES:
2849            sage: E = EllipticCurve([0, 0, 1, -38, 90])  # 361A
2850            sage: E.non_surjective()   # CM curve
2851            [(0, 'cm')]
2852
2853            sage: E = EllipticCurve([0, -1, 1, 0, 0]) # X_1(11)
2854            sage: E.non_surjective()
2855            [(5, '5-torsion')]
2856
2857            sage: E = EllipticCurve([0, 0, 1, -1, 0]) # 37A
2858            sage: E.non_surjective()
2859            []
2860
2861            sage: E = EllipticCurve([0,-1,1,-2,-1])   # 141C
2862            sage: E.non_surjective()
2863            [(13, [1])]
2864
2865        ALGORITHM:
2866            When p<=3 use division polynomials.  For 5 <= p <= B,
2867            where B is Cojocaru's bound, use the results in Section 2
2868            of Serre's inventiones paper"Sur Les Representations Modulaires Deg
2869            Degre 2 de Galqbar Over Q."
2870        """
2871        if self.has_cm():
2872            misc.verbose("cm curve")
2873            return [(0,"cm")]
2874        N = self.conductor()
2875        if self.is_semistable():
2876            C = 11
2877            misc.verbose("semistable -- so bound is 11")
2878        else:
2879            C = 1 + 4*sqrt(6)*int(N)/3 * sqrt(mul([1+1.0/int(p) for p,_ in factor(N)]))
2880            misc.verbose("conductor = %s, and bound is %s"%(N,C))
2881        C = 1 + 4*sqrt(6)*int(N)/3 * sqrt(mul([1+1.0/int(p) for p,_ in factor(N)]))
2882        misc.verbose("conductor = %s, and bound is %s"%(N,C))
2883        B = []
2884        p = 2
2885        while p <= C:
2886            t, v = self.is_surjective(p, A=A)
2887            misc.verbose("(%s,%s,%s)"%(p,t,v))
2888            if not t:
2889                B.append((p,v))
2890            p = next_prime(p)
2891        return B
2892
2893    def is_ordinary(self, p, ell=None):
2894        """
2895        Return True precisely when the mod-p representation attached
2896        to this elliptic curve is ordinary at ell.
2897
2898        INPUT:
2899            p -- a prime
2900            ell - a prime (default: p)
2901
2902        OUTPUT:
2903            bool
2904
2905        EXAMPLES:
2906            sage: E=EllipticCurve('37a1')
2907            sage: E.is_ordinary(37)
2908            True
2909            sage: E=EllipticCurve('32a1')
2910            sage: E.is_ordinary(2)
2911            False
2912            sage: [p for p in prime_range(50) if E.is_ordinary(p)]
2913            [5, 13, 17, 29, 37, 41]
2914        """
2915        if ell is None:
2916            ell = p
2917        return self.ap(ell) % p != 0
2918
2919    def is_good(self, p, check=True):
2920        """
2921        Return True if $p$ is a prime of good reduction for $E$.
2922
2923        INPUT:
2924            p -- a prime
2925
2926        OUTPUT:
2927            bool
2928
2929        EXAMPLES:
2930            sage: e = EllipticCurve('11a')
2931            sage: e.is_good(-8)
2932            Traceback (most recent call last):
2933            ...
2934            ValueError: p must be prime
2935            sage: e.is_good(-8, check=False)
2936            True
2937        """
2938        if check:
2939            if not arith.is_prime(p):
2940                raise ValueError, "p must be prime"
2941        return self.conductor() % p != 0
2942       
2943
2944    def is_supersingular(self, p, ell=None):
2945        """
2946        Return True precisely when p is a prime of good reduction and
2947        the mod-p representation attached to this elliptic curve is
2948        supersingular at ell.
2949
2950        INPUT:
2951            p -- a prime
2952            ell - a prime (default: p)
2953
2954        OUTPUT:
2955            bool
2956
2957        EXAMPLES:
2958            sage: E=EllipticCurve('37a1')
2959            sage: E.is_supersingular(37)
2960            False
2961            sage: E=EllipticCurve('32a1')
2962            sage: E.is_supersingular(2)
2963            False
2964            sage: E.is_supersingular(7)
2965            True
2966            sage: [p for p in prime_range(50) if E.is_supersingular(p)]
2967            [3, 7, 11, 19, 23, 31, 43, 47]
2968        """
2969        if ell is None:
2970            ell = p
2971        return self.is_good(p) and not self.is_ordinary(p, ell)
2972
2973    def supersingular_primes(self, B):
2974        """
2975        Return a list of all supersingular primes for this elliptic curve
2976        up to and possibly including B.
2977
2978        EXAMPLES:
2979            sage: e = EllipticCurve('11a')
2980            sage: e.aplist(20)
2981            [-2, -1, 1, -2, 1, 4, -2, 0]
2982            sage: e.supersingular_primes(1000)
2983            [2, 19, 29, 199, 569, 809]
2984
2985            sage: e = EllipticCurve('27a')
2986            sage: e.aplist(20)
2987            [0, 0, 0, -1, 0, 5, 0, -7]
2988            sage: e.supersingular_primes(97)
2989            [2, 5, 11, 17, 23, 29, 41, 47, 53, 59, 71, 83, 89]
2990            sage: e.ordinary_primes(97)
2991            [7, 13, 19, 31, 37, 43, 61, 67, 73, 79, 97]
2992            sage: e.supersingular_primes(3)
2993            [2]
2994            sage: e.supersingular_primes(2)
2995            [2]
2996            sage: e.supersingular_primes(1)
2997            []           
2998        """
2999        v = self.aplist(max(B, 3))
3000        P = arith.prime_range(max(B,3)+1)
3001        N = self.conductor()
3002        return [P[i] for i in [0,1] if P[i] <= B and v[i]%P[i]==0 and N%P[i] != 0] + \
3003                      [P[i] for i in range(2,len(v)) if v[i] == 0 and N%P[i] != 0]
3004
3005    def ordinary_primes(self, B):
3006        """
3007        Return a list of all ordinary primes for this elliptic curve
3008        up to and possibly including B.
3009
3010        EXAMPLES:
3011            sage: e = EllipticCurve('11a')
3012            sage: e.aplist(20)
3013            [-2, -1, 1, -2, 1, 4, -2, 0]
3014            sage: e.ordinary_primes(97)
3015            [3, 5, 7, 11, 13, 17, 23, 31, 37, 41, 43, 47, 53, 59, 61, 67, 71, 73, 79, 83, 89, 97]
3016            sage: e = EllipticCurve('49a')
3017            sage: e.aplist(20)
3018            [1, 0, 0, 0, 4, 0, 0, 0]
3019            sage: e.supersingular_primes(97)
3020            [3, 5, 13, 17, 19, 31, 41, 47, 59, 61, 73, 83, 89, 97]
3021            sage: e.ordinary_primes(97)
3022            [2, 11, 23, 29, 37, 43, 53, 67, 71, 79]
3023            sage: e.ordinary_primes(3)
3024            [2]
3025            sage: e.ordinary_primes(2)
3026            [2]
3027            sage: e.ordinary_primes(1)
3028            []
3029        """
3030        v = self.aplist(max(B, 3) )
3031        P = arith.prime_range(max(B,3) +1)
3032        return [P[i] for i in [0,1] if P[i] <= B and v[i]%P[i]!=0] +\
3033               [P[i] for i in range(2,len(v)) if v[i] != 0]
3034
3035    def eval_modular_form(self, points, prec):
3036        """
3037        Evaluate the L-series of this elliptic curve at points in CC
3038
3039        INPUT:
3040              points--  a list of points in the half-plane of convergence
3041              prec--    precision
3042
3043        OUTPUT:
3044              A list of values L(E,s) for s in points
3045
3046        NOTE:  Better examples are welcome.
3047
3048        EXAMPLES:
3049            sage: E=EllipticCurve('37a1')
3050            sage: E.eval_modular_form([1.5+I,2.0+I,2.5+I],0.000001)
3051            [0, 0, 0]
3052        """
3053        if not isinstance(points, (list,xrange)):
3054            try:
3055                points = list(points)
3056            except TypeError:
3057                return self.eval_modular_form([points],prec)
3058        an = self.pari_mincurve().ellan(prec)
3059        s = 0
3060        I = pari("I")
3061        pi = pari("Pi")
3062        c = pari(2)*pi*I
3063        ans = []
3064        for z in points:
3065            s = pari(0)
3066            r0 = (c*z).exp()
3067            r = r0
3068            for n in xrange(1,prec):
3069                s += an[n-1]*r     
3070                r *= r0
3071            ans.append(s.python())
3072        return ans
3073
3074    def _multiple_of_degree_of_isogeny_to_optimal_curve(self):
3075        """
3076        Internal function returning an integer m such that the degree of
3077        the isogeny between this curve and the optimal curve in its
3078        isogeny class is a divisor of m.
3079
3080        EXAMPLES:
3081            sage: E=EllipticCurve('11a1')
3082            sage: E._multiple_of_degree_of_isogeny_to_optimal_curve()
3083            25
3084            sage: E=EllipticCurve('11a2')
3085            sage: E._multiple_of_degree_of_isogeny_to_optimal_curve()
3086            5
3087            sage: E=EllipticCurve('11a3')
3088            sage: E._multiple_of_degree_of_isogeny_to_optimal_curve()
3089            5
3090        """
3091        M = self.isogeny_class()[1]
3092        return Integer(misc.prod([x for x in M.row(0) if x], 1))
3093               
3094    ########################################################################
3095    # Functions related to bounding the order of Sha (provably correctly!)
3096    # Heegner points and Kolyvagin's theorem
3097    ########################################################################
3098
3099    def sha(self):
3100        """
3101        Return an object of class
3102        'sage.schemes.elliptic_curves.sha.Sha' attached to this
3103        elliptic curve.
3104
3105        This can be used in functions related to bounding the order of
3106        Sha (The Tate-Shafarevich group of the curve).
3107
3108        EXAMPLES:
3109            sage: E=EllipticCurve('37a1')
3110            sage: S=E.sha()
3111            sage: S
3112            <class 'sage.schemes.elliptic_curves.sha.Sha'>
3113            sage: S.bound_kolyvagin()
3114            ([2], 1)
3115        """
3116        try:
3117            return self.__sha
3118        except AttributeError:
3119            from sha import Sha
3120            self.__sha = Sha(self)
3121            return self.__sha
3122   
3123    def satisfies_heegner_hypothesis(self, D):
3124        """
3125        Returns True precisely when D is a fundamental discriminant
3126        that satisfies the Heegner hypothesis for this elliptic curve.
3127
3128        EXAMPLES:
3129            sage: E = EllipticCurve('11a1')
3130            sage: E.satisfies_heegner_hypothesis(-7)
3131            True
3132            sage: E.satisfies_heegner_hypothesis(-11)
3133            False
3134        """
3135        if not number_field.is_fundamental_discriminant(D):
3136            return False
3137        if arith.GCD(D, self.conductor()) != 1:
3138            return False
3139        for p, _ in factor(self.conductor()):
3140            if kronecker_symbol(D,p) != 1:
3141                return False
3142        return True
3143
3144    def heegner_discriminants(self, bound):
3145        """
3146        Return the list of self's Heegner discriminants between -1 and -bound.
3147
3148        INPUT:
3149            bound (int) -- upper bound for -discriminant
3150
3151        OUTPUT:
3152            The list of Heegner discriminants between -1 and -bound for the given elliptic curve.
3153
3154        EXAMPLE:
3155            sage: E=EllipticCurve('11a')
3156            sage: E.heegner_discriminants(30)
3157            [-7, -8, -19, -24]
3158        """
3159        return [-D for D in xrange(1,bound) if self.satisfies_heegner_hypothesis(-D)]
3160
3161    def heegner_discriminants_list(self, n):
3162        """
3163        Return the list of self's first n Heegner discriminants smaller than -5.
3164
3165        INPUT:
3166            n (int) -- the number of discriminants to compute
3167
3168        OUTPUT:
3169            The list of the first n Heegner discriminants smaller than -5 for the given elliptic curve.
3170
3171        EXAMPLE:
3172            sage: E=EllipticCurve('11a')
3173            sage: E.heegner_discriminants_list(4)
3174            [-7, -8, -19, -24]
3175        """
3176        v = []
3177        D = -5
3178        while len(v) < n:
3179            while not self.satisfies_heegner_hypothesis(D):
3180                D -= 1
3181            v.append(D)
3182            D -= 1
3183        return v
3184
3185    def heegner_point_height(self, D, prec=2):
3186        """
3187        Use the Gross-Zagier formula to compute the Neron-Tate
3188        canonical height over K of the Heegner point corresponding to
3189        D, as an Interval (since it's computed to some precision using
3190        L-functions).
3191
3192        INPUT:
3193            D (int) -- fundamental discriminant (=/= -3, -4)
3194            prec (int) -- (default: 2), use prec*sqrt(N) + 20 terms
3195                          of L-series in computations, where N is the
3196                          conductor.
3197
3198        OUTPUT:
3199            Interval that contains the height of the Heegner point.
3200
3201        EXAMPLE:
3202            sage: E = EllipticCurve('11a')
3203            sage: E.heegner_point_height(-7)
3204            [0.22226977 ... 0.22227479]
3205        """
3206
3207        if not self.satisfies_heegner_hypothesis(D):
3208            raise ArithmeticError, "Discriminant (=%s) must be a fundamental discriminant that satisfies the Heegner hypothesis."%D
3209        if D == -3 or D == -4:
3210            raise ArithmeticError, "Discriminant (=%s) must not be -3 or -4."%D
3211        eps = self.root_number()
3212        L1_vanishes = self.lseries().L1_vanishes()
3213        if eps == 1 and L1_vanishes:
3214            return IR(0) # rank even hence >= 2, so Heegner point is torsion.
3215        alpha = R(sqrt(abs(D)))/(2*self.period_lattice().complex_area())
3216        F = self.quadratic_twist(D)
3217        E = self
3218        k_E = prec*sqrt(E.conductor()) + 20
3219        k_F = prec*sqrt(F.conductor()) + 20
3220
3221        MIN_ERR = R('1e-6')  # we assume that regulator and
3222                             # discriminant, etc., computed to this accuracy (which is easily the case).
3223                             # this should be made more intelligent / rigorous relative
3224                             # to the rest of the system.
3225        if eps == 1:   # E has even rank
3226            LF1, err_F = F.lseries().deriv_at1(k_F)
3227            LE1, err_E = E.lseries().at1(k_E)
3228            err_F = max(err_F, MIN_ERR)
3229            err_E = max(err_E, MIN_ERR)
3230            return IR(alpha-MIN_ERR,alpha+MIN_ERR) * IR(LE1-err_E,LE1+err_E) * IR(LF1-err_F,LF1+err_F)
3231
3232        else:          # E has odd rank
3233            LE1, err_E = E.lseries().deriv_at1(k_E)
3234            LF1, err_F = F.lseries().at1(k_F)
3235            err_F = max(err_F, MIN_ERR)
3236            err_E = max(err_E, MIN_ERR)
3237            return IR(alpha-MIN_ERR,alpha+MIN_ERR) * IR(LE1-err_E,LE1+err_E) * IR(LF1-err_F,LF1+err_F)
3238       
3239
3240    def heegner_index(self, D,  min_p=2, prec=5, verbose=False):
3241        r"""
3242        Return an interval that contains the index of the Heegner
3243        point $y_K$ in the group of K-rational points modulo torsion
3244        on this elliptic curve, computed using the Gross-Zagier
3245        formula and/or a point search, or the index divided by $2$.
3246
3247        NOTES: If \code{min_p} is bigger than 2 then the index can be
3248        off by any prime less than \code{min_p}.   This function
3249        returns the index divided by $2$ exactly when $E(\Q)_{/tor}$
3250        has index $2$ in $E(K)_{/tor}$.
3251
3252        INPUT:
3253            D (int) -- Heegner discriminant
3254            min_p (int) -- (default: 2) only rule out primes >= min_p
3255                           dividing the index. 
3256            verbose (bool) -- (default: False); print lots of mwrank
3257                              search status information when computing
3258                              regulator
3259            prec (int) -- (default: 5), use prec*sqrt(N) + 20 terms
3260                          of L-series in computations, where N is the
3261                          conductor.
3262                         
3263        OUTPUT:
3264            an interval that contains the index
3265
3266        EXAMPLES:
3267            sage: E = EllipticCurve('11a')
3268            sage: E.heegner_discriminants(50)
3269            [-7, -8, -19, -24, -35, -39, -40, -43]
3270            sage: E.heegner_index(-7)
3271            [0.99999332 .. 1.0000077]
3272
3273            sage: E = EllipticCurve('37b')
3274            sage: E.heegner_discriminants(100)
3275            [-3, -4, -7, -11, -40, -47, -67, -71, -83, -84, -95]
3276            sage: E.heegner_index(-95)          # long time (1 second)
3277            [1.9999923 .. 2.0000077]
3278
3279        Current discriminants -3 and -4 are not supported:
3280            sage: E.heegner_index(-3)
3281            Traceback (most recent call last):
3282            ...
3283            ArithmeticError: Discriminant (=-3) must not be -3 or -4.
3284
3285        The curve 681b returns an interval that contains $3/2$.
3286        This is because $E(\Q)$ is not saturated in $E(K)$.  The
3287        true index is $3$:
3288            sage: E = EllipticCurve('681b')
3289            sage: I = E.heegner_index(-8); I
3290            [1.4999942 .. 1.5000058]
3291            sage: 2*I
3292            [2.9999885 .. 3.0000115]
3293
3294        In fact, whenever the returned index has a denominator
3295        of $2$, the true index is got by multiplying the returned
3296        index by $2$.  Unfortunately, this is not an if and only
3297        if condition, i.e., sometimes the index must be multiplied
3298        by $2$ even though the denominator is not $2$.
3299        """
3300        # First compute upper bound on height of Heegner point.
3301        tm = misc.verbose("computing heegner point height...")
3302        h0 = self.heegner_point_height(D, prec=prec)
3303       
3304        # We divide by 2 to get the height **over Q** of the
3305        # Heegner point on the twist.
3306
3307        ht = h0/2
3308        misc.verbose('Height of heegner point = %s'%ht, tm)
3309       
3310        if self.root_number() == 1:
3311            F = self.quadratic_twist(D)
3312        else:
3313            F = self
3314        h  = ht.upper()
3315        misc.verbose("Heegner height bound = %s"%h)
3316        B = F.CPS_height_bound()
3317        misc.verbose("CPS bound = %s"%B)
3318        c = h/(min_p**2) + B
3319        misc.verbose("Search would have to be up to height = %s"%c)
3320       
3321        if c > _MAX_HEIGHT or F is self:
3322            misc.verbose("Doing direct computation of MW group.")
3323            reg = F.regulator(verbose=verbose)
3324            return self.__adjust_heegner_index(ht/IR(reg))
3325       
3326        # Do naive search to eliminate possibility that Heegner point
3327        # is divisible by p<min_p, without finding Heegner point.
3328        misc.verbose("doing point search")
3329        P = F.point_search(c, verbose=verbose)
3330        misc.verbose("done with point search")
3331        P = [x for x in P if x.order() == oo]
3332        if len(P) == 0:
3333            return IR(1)
3334
3335        misc.verbose("saturating")
3336        S, I, reg = F.saturation(P, verbose=verbose)
3337        misc.verbose("done saturating")
3338        return self.__adjust_heegner_index(ht/IR(reg))
3339
3340    def __adjust_heegner_index(self, a):
3341        r"""
3342        Take the square root of the interval that contains the Heegner
3343        index.
3344
3345        EXAMPLES:
3346            sage: E = EllipticCurve('11a1')
3347            sage: a = RIF(sqrt(2))-1.4142135623730951
3348            sage: E._EllipticCurve_rational_field__adjust_heegner_index(a)
3349            [0.0000000... .. 1.490116...e-8]
3350        """
3351        if a.lower() < 0:
3352            a = IR((0, a.upper()))
3353        return a.sqrt()
3354       
3355
3356    def heegner_index_bound(self, D=0,  prec=5, verbose=True, max_height=_MAX_HEIGHT):
3357        """
3358        Assume self has rank 0.
3359       
3360        Return a list v of primes such that if an odd prime p divides
3361        the index of the Heegner point in the group of rational
3362        points *modulo torsion*, then p is in v.
3363
3364        If 0 is in the interval of the height of the Heegner point
3365        computed to the given prec, then this function returns v = 0.
3366        This does not mean that the Heegner point is torsion, just
3367        that it is very likely torsion.
3368
3369        If we obtain no information from a search up to max_height, e.g.,
3370        if the Siksek et al. bound is bigger than max_height, then
3371        we return v = -1.
3372
3373        INPUT:
3374            D (int) -- (deault: 0) Heegner discriminant; if 0, use the
3375                       first discriminant < -4 that satisfies the Heegner
3376                       hypothesis
3377            verbose (bool) -- (default: True)
3378            prec (int) -- (default: 5), use prec*sqrt(N) + 20 terms
3379                          of L-series in computations, where N is the
3380                          conductor.
3381            max_height (float) -- should be <= 21; bound on logarithmic
3382                                  naive height used in point searches.
3383                                  Make smaller to make this function
3384                                  faster, at the expense of possibly
3385                                  obtaining a worse answer.  A good
3386                                  range is between 13 and 21.
3387                                 
3388        OUTPUT:
3389            v -- list or int (bad primes or 0 or -1)
3390            D -- the discriminant that was used (this is useful if D was
3391                 automatically selected).
3392
3393        EXAMPLES:
3394            sage: E = EllipticCurve('11a1')
3395            sage: E.heegner_index_bound(verbose=False)
3396            ([2], -7)
3397        """
3398        max_height = min(float(max_height), _MAX_HEIGHT)
3399        if self.root_number() != 1:
3400            raise RuntimeError, "The rank must be 0."
3401       
3402        if D == 0:
3403            D = -5
3404            while not self.satisfies_heegner_hypothesis(D):
3405                D -= 1
3406       
3407        # First compute upper bound on Height of Heegner point.
3408        ht = self.heegner_point_height(D, prec=prec)
3409        if 0 in ht:
3410            return 0, D
3411        F = self.quadratic_twist(D)
3412        h  = ht.upper()
3413        misc.verbose("Heegner height bound = %s"%h)
3414        B = F.CPS_height_bound()
3415        misc.verbose("CPS bound = %s"%B)
3416        H = h
3417        p = 3
3418        while True:
3419            c = h/(2*p**2) + B
3420            if c < max_height:
3421                break
3422            if p > 100:
3423                break
3424            p = next_prime(p)
3425        misc.verbose("Using p = %s"%p)
3426
3427        if c > max_height:
3428            misc.verbose("No information by searching only up to max_height (=%s)."%c)
3429            return -1, D
3430           
3431        misc.verbose("Searching up to height = %s"%c)
3432        eps = 10e-5
3433
3434        def _bound(P):
3435            """
3436            We will use this function below in two places.  It
3437            bounds the index using a nontrivial point.
3438            """
3439            assert len(P) == 1
3440           
3441            S, I, reg = F.saturation(P, verbose=verbose)
3442            h = IR(reg-eps,reg+eps)           
3443            ind2 = ht/(h/2)
3444            misc.verbose("index squared = %s"%ind2)
3445            ind = ind2.sqrt()
3446            misc.verbose("index = %s"%ind)
3447            # Compute upper bound on square root of index.
3448            if ind.absolute_diameter() < 1:
3449                t, i = ind.is_int()
3450                if t:   # unique integer in interval, so we've found exact index squared.
3451                    return arith.prime_divisors(i), D
3452            raise RuntimeError, "Unable to compute bound for e=%s, D=%s (try increasing precision)"%(self,D)
3453       
3454        # First try a quick search, in case we get lucky and find
3455        # a generator.
3456        P = F.point_search(13, verbose=verbose)
3457        P = [x for x in P if x.order() == oo]
3458        if len(P) > 0:
3459            return _bound(P)
3460
3461        # Do search to eliminate possibility that Heegner point is
3462        # divisible by primes up to p, without finding Heegner point.
3463        P = F.point_search(c, verbose=verbose)
3464        P = [x for x in P if x.order() == oo]
3465        if len(P) == 0:
3466            # We've eliminated the possibility of a divisor up to p.
3467            return arith.prime_range(3,p), D
3468        else:
3469            return _bound(P)
3470
3471    padic_regulator = padics.padic_regulator
3472
3473    padic_height_pairing_matrix = padics.padic_height_pairing_matrix
3474
3475    padic_height = padics.padic_height
3476    padic_height_via_multiply = padics.padic_height_via_multiply
3477
3478    padic_sigma = padics.padic_sigma
3479    padic_sigma_truncated = padics.padic_sigma_truncated
3480       
3481    padic_E2 = padics.padic_E2
3482       
3483    matrix_of_frobenius = padics.matrix_of_frobenius
3484   
3485    # def weierstrass_p(self):
3486    #         # TODO: add allowing negative valuations for power series
3487    #         return 1/t**2 + a1/t + rings.frac(1,12)*(a1-8*a2) -a3*t \
3488    #                - (a4+a1*a3)*t**2  + O(t**3)
3489
3490
3491    def mod5family(self):
3492        """
3493        Return the family of all elliptic curves with the same mod-5
3494        representation as self.
3495
3496        EXAMPLES:
3497        sage: E=EllipticCurve('32a1')
3498        sage: E.mod5family()
3499        Elliptic Curve defined by y^2  = x^3 + 4*x over Fraction Field of Univariate Polynomial Ring in t over Rational Field
3500        """
3501        E = self.short_weierstrass_model()
3502        a = E.a4()
3503        b = E.a6()
3504        return mod5family.mod5family(a,b)
3505
3506    def tate_curve(self, p):
3507        r"""
3508        Creates the Tate Curve over the $p$-adics associated to this elliptic curves.
3509
3510        This Tate curve a $p$-adic curve with split multiplicative
3511        reduction of the form $y^2+xy=x^3+s_4 x+s_6$ which is
3512        isomorphic to the given curve over the algebraic closure of
3513        $\QQ_p$.  Its points over $\QQ_p$ are isomorphic to
3514        $\QQ_p^{\times}/q^{\Z}$ for a certain parameter $q\in\Z_p$.
3515       
3516        INPUT:
3517
3518            p -- a prime where the curve has multiplicative reduction.
3519           
3520
3521        EXAMPLES:
3522            sage: e = EllipticCurve('130a1')
3523            sage: e.tate_curve(2)
3524            2-adic Tate curve associated to the Elliptic Curve defined by y^2 + x*y + y = x^3 - 33*x + 68 over Rational Field
3525
3526        The input curve must have multiplicative reduction at the prime.
3527            sage: e.tate_curve(3)
3528            Traceback (most recent call last):
3529            ...
3530            ValueError: The elliptic curve must have multiplicative reduction at 3
3531
3532        We compute with $p=5$:           
3533            sage: T = e.tate_curve(5); T
3534            5-adic Tate curve associated to the Elliptic Curve defined by y^2 + x*y + y = x^3 - 33*x + 68 over Rational Field
3535
3536        We find the Tate parameter $q$:
3537            sage: T.parameter(prec=5)
3538            3*5^3 + 3*5^4 + 2*5^5 + 2*5^6 + 3*5^7 + O(5^8)
3539
3540        We compute the $L$-invariant of the curve:
3541            sage: T.L_invariant(prec=10)
3542            5^3 + 4*5^4 + 2*5^5 + 2*5^6 + 2*5^7 + 3*5^8 + 5^9 + O(5^10)
3543        """
3544        try:
3545            return self._tate_curve[p]
3546        except AttributeError:
3547            self._tate_curve = {}
3548        except KeyError:
3549            pass
3550       
3551        Eq = ell_tate_curve.TateCurve(self,p)
3552        self._tate_curve[p] = Eq
3553        return Eq
3554
3555
3556    def integral_points(self, mw_base='auto',tors_points='auto'):
3557        """
3558        Computes all integral points on the elliptic curve E which has Mordell-Weil
3559        base mw_base and torsions points specified by tors_points.
3560        Default 'auto' relies on self.gens() and self.torsion_subgroup().gens()
3561       
3562        INPUT:
3563            self -- EllipticCurve_Rational_Field
3564            mw_base -- list of EllipticCurvePoint representing the Mordell-Weil base
3565                    (default: 'auto' - calls self.gens())
3566            tors_points -- list of EllipticCurvePoint giving all torsion points of E
3567                        (default: 'auto' - calls self.torsion_subgroup().gens())
3568           
3569        OUTPUT:
3570            set of all integral points on E
3571           
3572        HINTS:
3573            - The complexity increases exponentially in the rank of curve E.
3574            - It can help if you try another Mordell-Weil base, because the
3575            computation time depends on this, too.   
3576           
3577        EXAMPLES:
3578        A curve of rank 3 with no torsion points
3579       
3580            sage: E=EllipticCurve([0,0,1,-7,6])
3581            sage: P1=E.point((2,0)); P2=E.point((-1,3)); P3=E.point((4,6))
3582            sage: a=E.integral_points([P1,P2,P3], []); a 
3583            set([(3 : -4 : 1), (21 : 95 : 1), (-3 : 0 : 1), (21 : -96 : 1), (93 : -897 : 1), (3 : 3 : 1), (-1 : 3 : 1), (8 : 21 : 1), (-3 : -1 : 1), (342 : 6324 : 1), (37 : -225 : 1), (37 : 224 : 1), (11 : -36 : 1), (406 : -8181 : 1), (4 : -7 : 1), (2 : 0 : 1), (52 : 374 : 1), (93 : 896 : 1), (14 : 51 : 1), (14 : -52 : 1), (8 : -22 : 1), (4 : 6 : 1), (-2 : -4 : 1), (816 : -23310 : 1), (11 : 35 : 1), (2 : -1 : 1), (0 : 2 : 1), (1 : -1 : 1), (-2 : 3 : 1), (0 : -3 : 1), (406 : 8180 : 1), (52 : -375 : 1), (816 : 23309 : 1), (-1 : -4 : 1), (1 : 0 : 1), (342 : -6325 : 1)])
3584           
3585        It is also possible not to specify mw_base and tors_points but because the
3586        Mordell-Weil base found by built-in functions differs this takes longer
3587       
3588            sage: E=EllipticCurve([0,0,1,-7,6])
3589            sage: a=E.integral_points(); a
3590            set([(3 : -4 : 1), (21 : 95 : 1), (-3 : 0 : 1), (21 : -96 : 1), (93 : -897 : 1), (3 : 3 : 1), (-1 : 3 : 1), (4 : 6 : 1), (-3 : -1 : 1), (342 : 6324 : 1), (37 : -225 : 1), (37 : 224 : 1), (11 : -36 : 1), (406 : -8181 : 1), (4 : -7 : 1), (2 : 0 : 1), (52 : 374 : 1), (93 : 896 : 1), (14 : 51 : 1), (14 : -52 : 1), (8 : -22 : 1), (8 : 21 : 1), (-2 : -4 : 1), (816 : -23310 : 1), (11 : 35 : 1), (2 : -1 : 1), (0 : 2 : 1), (1 : -1 : 1), (-2 : 3 : 1), (0 : -3 : 1), (406 : 8180 : 1), (52 : -375 : 1), (816 : 23309 : 1), (-1 : -4 : 1), (1 : 0 : 1), (342 : -6325 : 1)])
3591           
3592       
3593        Another example with rank 5 and no torsion points
3594           
3595            sage: E=EllipticCurve([-879984,319138704])
3596            sage: P1=E.point((540,1188)); P2=E.point((576,1836))
3597            sage: P3=E.point((468,3132)); P4=E.point((612,3132))
3598            sage: P5=E.point((432,4428))
3599            sage: a=E.integral_points([P1,P2,P3,P4,P5],[]); a        # long time (!)
3600            set([(155356 : -61232804 : 1), (2717460 : -4479656508 : 1), (468 : -3132 : 1), (513 : -1647 : 1), (-764 : 23356 : 1), (155356 : 61232804 : 1), (65052 : -16589988 : 1), (792 : 10908 : 1), (2204663452 : -103517423460188 : 1), (36828 : -7065252 : 1), (1732 : 63172 : 1), (-792 : -22788 : 1), (673 : -5633 : 1), (-792 : 22788 : 1), (2592 : 124308 : 1), (-620 : 25028 : 1), (720 : -7668 : 1), (1368 : 40932 : 1), (-620 : -25028 : 1), (540 : -1188 : 1), (1072 : 24652 : 1), (108 : 15012 : 1), (1260648 : -1415437308 : 1), (612 : -3132 : 1), (576 : -1836 : 1), (1072 : -24652 : 1), (1260648 : 1415437308 : 1), (-72 : -19548 : 1), (432 : 4428 : 1), (396 : -5724 : 1), (520 : -1468 : 1), (193320 : 84998268 : 1), (-684 : 24516 : 1), (972 : 19548 : 1), (-900 : 19548 : 1), (-423 : -24813 : 1), (43200 : 8976852 : 1), (9972 : -991548 : 1), (14265 : 1700163 : 1), (78696 : 22074876 : 1), (376 : 6436 : 1), (196992 : 87431508 : 1), (-684 : -24516 : 1), (585 : 2133 : 1), (792 : -10908 : 1), (4680 : 314172 : 1), (-1080 : 3132 : 1), (2448 : -113292 : 1), (1732 : -63172 : 1), (-279 : 23301 : 1), (193320 : -84998268 : 1), (3060 : -162108 : 1), (-72 : 19548 : 1), (172 : -13148 : 1), (376 : -6436 : 1), (2385 : 108567 : 1), (36828 : 7065252 : 1), (114280 : -38631428 : 1), (612 : 3132 : 1), (540 : 1188 : 1), (114280 : 38631428 : 1), (4680 : -314172 : 1), (5940 : -452412 : 1), (2448 : 113292 : 1), (-423 : 24813 : 1), (-864 : 20844 : 1), (-279 : -23301 : 1), (8668 : 802468 : 1), (720 : 7668 : 1), (17856 : 2382804 : 1), (2385 : -108567 : 1), (-1080 : -3132 : 1), (432 : -4428 : 1), (3060 : 162108 : 1), (513 : 1647 : 1), (972 : -19548 : 1), (3537 : 203607 : 1), (368892 : 224051292 : 1), (468 : 3132 : 1), (396 : 5724 : 1), (17856 : -2382804 : 1), (2113 : -88847 : 1), (9972 : 991548 : 1), (14265 : -1700163 : 1), (43200 : -8976852 : 1), (36 : 16956 : 1), (108 : -15012 : 1), (-864 : -20844 : 1), (-900 : -19548 : 1), (2717460 : 4479656508 : 1), (585 : -2133 : 1), (8668 : -802468 : 1), (2204663452 : 103517423460188 : 1), (196992 : -87431508 : 1), (1368 : -40932 : 1), (78696 : -22074876 : 1), (172 : 13148 : 1), (576 : 1836 : 1), (3537 : -203607 : 1), (2592 : -124308 : 1), (520 : 1468 : 1), (2113 : 88847 : 1), (368892 : -224051292 : 1), (65052 : 16589988 : 1), (36 : -16956 : 1), (673 : 5633 : 1), (-764 : -23356 : 1), (5940 : 452412 : 1)])
3601           
3602           
3603        NOTES:
3604            - This function uses the algorithm given in [Co1] and the computation of
3605            the complex elliptic logarithm as presented in [Co2].
3606            REFERENCES:
3607                -[Co1] Cohen H., Number Theory Vol I: Tools and Diophantine Equations
3608                        GTM 239, Springer 2007
3609                -[Co2] Cohen H., A Course in Computational Algebraic Number Theory
3610                        GTM 138, Springer 1996
3611        AUTHORS:
3612            - Michael Mardaus (2008-07)
3613            - Tobias Nagel (2008-07)
3614        """
3615    ###############################################################################
3616    # INPUT CHECK #################################################################
3617        if mw_base=='auto':
3618            mw_base = self.gens()
3619        if tors_points=='auto':
3620            tors_points = self.torsion_subgroup().gens()
3621        try:
3622            len_tors = len(tors_points)
3623        except TypeError:
3624            raise TypeError, "tors_points must be a list"
3625        try:
3626            r = len(mw_base)
3627        except TypeError:
3628            raise TypeError, 'mw_base must be a list'
3629        if (r == 0) and (len_tors == 0):
3630            raise RuntimeError, 'Both base points and torsions points are not specified'
3631        try: # test if all given points are on curve (try not necessary because of no new error)
3632            if r < len_tors:
3633                for i in range(r):
3634                    trash = self.point(mw_base[i])
3635                    trash = self.point(tors_points[i])
3636                for i in range(r,len_tors):
3637                    trash = self.point(tors_points[i])
3638            else:
3639                for i in range(len_tors):
3640                    trash = self.point(mw_base[i])
3641                    trash = self.point(tors_points[i])
3642                for i in range(len_tors,r):
3643                    trash = self.point(mw_base[i])
3644        except TypeError:
3645            raise
3646    #END INPUT-CHECK###############################################################
3647    ###############################################################################
3648               
3649    ###############################################################################
3650    # INTERNAL FUNCTIONS ##########################################################
3651   
3652    ############################## begin ################################
3653        def is_int(a):
3654            if floor(a)==ceil(a):
3655                return 1
3656            else:
3657                return 0
3658    ##############################  end  ################################           
3659    ############################## begin ################################
3660        def height_of_curve():
3661            #Computes height of curve 'self' according to Co1
3662            if j == 0:
3663                h_j = 1
3664            else:
3665                h_j = max(log(abs(j.numerator()), e), log(abs(j.denominator()), e))
3666            if (g2 != 0) and (g3 != 0):
3667                h_gs = max(1, log(abs(g2), e), log(abs(g3), e))
3668            elif g2 == 0:
3669                if g3 == 0:
3670                    return max(1,h_j)
3671                h_gs = max(1, log(abs(g3), e))
3672            else:
3673                h_gs = max(1, log(abs(g2), e))
3674            return max(1,h_j, h_gs)
3675    ##############################  end  ################################
3676    ############################## begin ################################
3677        def log_plus(x):
3678            if x==0:
3679                return 1
3680            else:
3681                return max(1, log(abs(x), e))
3682    ##############################  end  ################################
3683    ############################## begin ################################
3684        def height_paaring_matrix(list):
3685            #computes the height paaring matrix associated with the bilinear form of the
3686            #canonical (Nero-Tate) height function
3687            point_height = []
3688            dim=len(list)
3689            for i in range(dim):
3690                point_height.append(list[i].height())
3691            M = matrix.MatrixSpace(R, dim)
3692            mat = M()
3693            for j in range(dim):
3694                for k in range(j,dim):             
3695                    mat[j,k]=((list[j]+list[k]).height() - point_height[j] - point_height[k])/2
3696                    if j!=k:
3697                        mat[k,j]=mat[j,k]
3698            return mat
3699    ##############################  end  ################################
3700    ############################## begin ################################
3701        def min_eigenvalue(M):
3702            #computes the minimal eigenvalue of M
3703            eigenvalues = []
3704            eigenvalues =M.charpoly().real_roots();
3705            return min(eigenvalues)
3706    ##############################  end  ################################
3707    ############################## begin ################################
3708        def prod(list):
3709            #multiplies all elements in list
3710            tmp = 1
3711            for i in range(len(list)):
3712                tmp = tmp * list[i]
3713            return tmp
3714    ##############################  end  ################################
3715    ############################## begin ################################
3716        def extract_realroots(list, eps = 1e-8):
3717            #function is used to extract real roots from result of
3718            #function roots() listed in list
3719            #at the same time cutting off impreciseness of roots()
3720            #RETURN: sorted list of roots in RR
3721            tmp = []
3722            for i in range(0, len(list)):
3723                im = list[i][0].imag()
3724                re = list[i][0].real()
3725                if ((im < eps) and (im > -eps)): #imaginary part in eps
3726                    if ((re < eps) and (re > -eps)): #both parts are imprecise (also 0 included)
3727                        tmp.append(0) 
3728                    else:
3729                        tmp.append(re)
3730            tmp.sort()
3731            return tmp
3732    ##############################  end  ################################
3733    ############################## begin ################################
3734        def extract_roots(list, eps = 1e-8):
3735            #function is used to extract roots from result of
3736            #function roots() listed in list
3737            #cutting off impreciseness of roots()
3738            #RETURN: sorted list of roots in C
3739            tmp = []
3740            for i in range(0, len(list)):
3741                im = list[i][0].imag()
3742                re = list[i][0].real()
3743                if ((im < eps) and (im > -eps)):
3744                    if ((re < eps) and (re > -eps)): #both parts are imprecise (also 0 included)
3745                        tmp.append(0 + 0*I)
3746                    else:
3747                        tmp.append(re + 0*I)
3748                elif ((re < eps) and (re > -eps)):
3749                    tmp.append(im*I)
3750                else: #no impreciseness
3751                    tmp.append(list[i][0])
3752            return tmp
3753    ##############################  end  ################################
3754    ############################## begin ################################
3755        def in_egg(P,e1,e2):
3756            #test if P is on compact component of E
3757            #only possible for disc(E)>0, otherwise there is no compact component
3758            if (P.xy()[0] >= e1) and (P.xy()[0] <= e2):
3759                return True
3760            else:
3761                return False
3762    ##############################  end  ################################
3763    ############################## begin ################################
3764        def point_preprocessing(list,e1,e2):
3765            #Transforms mw_base so that at most one point is on the
3766            #compact component of the curve
3767            Q = []
3768            mem = -1
3769            for i in range(0,len(list)):
3770                if in_egg(list[i],e1,e2):
3771                    if mem == -1:
3772                        mem = i
3773                    else:
3774                        Q.append(list[i]+list[mem])
3775                        mem = i
3776                else:
3777                    Q.append(list[i])
3778            if mem != -1: #add last point, which is not in egg, to Q
3779                Q.append(2*list[mem])
3780            return Q
3781    ##############################  end  ################################
3782    ############################## begin ################################
3783        def modified_height(i):#computes modified height if base point i
3784            return max(mw_base[i].height(),h_E,c7*mw_base_log[i]**2)
3785    ##############################  end  ################################
3786    ############################## begin ################################
3787        def complex_elliptic_log(E,P):
3788                #Algorithm according to Co2
3789                x = P.xy()[0]
3790                y = P.xy()[1]
3791                #Initialize
3792                a1 = E.a1()
3793                a3 = E.a3()
3794                b2 = E.b2()
3795                b4 = E.b4()
3796                b6 = E.b6()
3797                if disc < 0:
3798                #Connected Case
3799                    PZ = rings.PolynomialRing(C, 'var')
3800                    var = PZ.gen()
3801                    tmp = (4*var**3+b2*var**2+2*b4*var+b6).roots()
3802                    real_roots = extract_realroots(tmp) 
3803                    if len(real_roots) != 1:
3804                        raise ValueError, ' no or more than one real root despite disc < 0'
3805                    else:
3806                        e1 = real_roots[0]
3807                        beta = sqrt(3*e1**2 + (b2/2)*e1 + b4/2)
3808                        alpha = 3*e1 + b2/4
3809                        a = R(2*sqrt(beta))
3810                        b = R(sqrt(alpha + 2*beta))
3811                        c = R((x - e1 + beta)/sqrt(x - e1))
3812                        while (a - b) > 1e-12:
3813                            a = R((a + b)/2)
3814                            b = R(sqrt(a*b))
3815                            c = R((c + sqrt(c**2 + b**2 - a**2))/2)
3816                        if (2*y + a1*x + a3)*((x - e1)**2 - beta**2) < 0:
3817                            z = R(arcsin(a/c)/a)
3818                        else:
3819                            z = R((pi - arcsin(a/c))/a)
3820                        if (2*y + a1*x + a3) > 0:
3821                            z = z + pi/a
3822                        return z
3823                else:
3824                    a2 = E.a2()
3825                    #Disconnected Case
3826                    PZ = rings.PolynomialRing(C, 'var')
3827                    var = PZ.gen()
3828                    tmp = (4*var**3+b2*var**2+2*b4*var+b6).roots()
3829                    real_roots = extract_realroots(tmp)
3830                    if len(real_roots)!=3:
3831                        raise ValueError, 'Not excatly 3 real roots found despite disc > 0'
3832                    else:
3833                        e1 = real_roots[2]
3834                        e2 = real_roots[1]
3835                        e3 = real_roots[0]
3836                        a = R(sqrt(e1 - e3))
3837                        b = R(sqrt(e1 - e2))
3838                        if x < e1:
3839                            f = 1
3840                            lam = y/(x - e3)
3841                            x = lam**2 + a1*lam - a2 - x - e3
3842                        else:
3843                            f = 0
3844                        c = sqrt(x - e3)
3845                        while (a-b) > 1e-12:
3846                            a = R((a + b)/2)
3847                            b = R(sqrt(a*b))
3848                            c = R((c + sqrt(c**2 + b**2 - a**2))/2)
3849                        #Connected component
3850                        if (((f == 0) and (2*y + a1*x + a3) < 0) or ((f == 1) and (2*y + a1*x + a3) >= 0)):
3851                            z = R(arcsin(a/c)/a)
3852                        else:
3853                            z = R((pi - arcsin(a/c))/a)
3854                        if f == 0:
3855                            return z
3856                        else:
3857                            #other component
3858                            w2=E.period_lattice().basis()[1]
3859                            return (z + w2/2)
3860    ##############################  end  ################################
3861    ############################## begin ################################
3862        def search_points(lowerb, upperb):
3863        # Search of integral points P on E with x(P) between lowerbd
3864        # and upperbd, using the explicit defintion of E
3865            a = self.a_invariants()
3866            for x in range(lowerb,upperb):
3867                rs = x**3 + a[1]*x**2 + a[3]*x + a[4] # right hand side of E
3868                y_coef = a[0]*x + a[2] #coefficient at y
3869                disc = y_coef**2 + 4*rs
3870                if disc > 0:
3871                    sol1 = R((-y_coef + sqrt(disc))/2)
3872                    sol2 = R((-y_coef - sqrt(disc))/2)
3873                    if is_int(sol1): #tests if sol1 is integral.
3874                        int_points.add(self.point((x,sol1)))
3875                    if is_int(sol2):
3876                        int_points.add(self.point((x,sol2))) 
3877                elif disc == 0:
3878                    sol1 = R((-y_coef + sqrt(disc))/2)
3879                    if is_int(sol1): #tests if sol1 is integral.
3880                        int_points.add(self.point((x,sol1)))
3881    ##############################  end  ################################
3882    ############################## begin ################################
3883        def search_remaining_points():
3884        #search integral points on curve E written as linear combination
3885        #of n times the mordell-weil base points and torsion points
3886        #(n is bounded by H_q, which will be computed at runtime)
3887            ###################### begin ########################
3888            #Computes b-adic representation of num
3889            #OUTPUT: list with b-adic representation list[0]=base^0 ... list[i]=base^i
3890            def badic_repres(num,base):
3891                new_num=[]
3892                for i in range(r):
3893                    new_num.append(num%base)
3894                    num=num//base
3895                return new_num
3896            ######################  end  ########################
3897            OP=self.point((0,1,0))
3898            #qw as tmp memory, for coefficients in linear combination
3899            #qw = []
3900            #for go in range(r):
3901            #    qw.append('init')
3902            #end
3903            for i in range(1,((2*H_q+1)**r)/2):
3904                koeffs = badic_repres(i,2*H_q+1)
3905                P = self.point((0,1,0))
3906                for index in range(r):
3907                    P = P + (koeffs[index]-H_q)*mw_base[index]
3908                    #qw[index]=koeffs[index]-H_q #s.a.
3909                if P!=OP:#P itself could be an integral point (torsion point would be (0:1:0) )
3910                    if is_int(P.xy()[0]) and is_int(P.xy()[1]):
3911                        int_points.add(P)
3912                        int_points.add(-P)
3913                        #print P, qw
3914                for j in range(len_tors): # P + torsion points
3915                    P = P + tors_points[j]
3916                    if P!=OP:
3917                        if is_int(P.xy()[0]) and is_int(P.xy()[1]):
3918                            int_points.add(P)
3919                            int_points.add(-P)
3920    ##############################  end  ################################
3921   
3922    # END Internal functions ######################################################
3923    ###############################################################################
3924   
3925        e = exp(1)
3926        I = constants.I
3927        pi = R(constants.pi)
3928   
3929        g2 = R((self.c4()/12))
3930        g3 = R((self.c6()/216))
3931        disc = self.discriminant()
3932        j = self.j_invariant()
3933        b2 = self.b2() 
3934       
3935        PR = rings.PolynomialRing(C, 'var')
3936        var = PR.gen()
3937        tmp = (4*var**3-g2*var-g3).roots()
3938        if disc > 0: #compact and noncompact component -> 3 roots in RR
3939            roots = extract_realroots(tmp)
3940            e1 = roots[0]
3941            e2 = roots[1]
3942            e3 = roots[2]
3943            mw_base = point_preprocessing(mw_base,e1,e2) #at most one point in E^{egg}
3944   
3945        elif disc < 0: #only noncompact component => 1 root in RR (=: e3), 2 roots in C
3946            roots = extract_roots(tmp)
3947            if roots[0].imag() == 0:
3948                e3 = roots[0].real()
3949                e2 = roots[1]
3950                e1 = roots[2]
3951            elif roots[1].imag() == 0:
3952                e3 = roots[1].real()
3953                e2 = roots[0]
3954                e1 = roots[2]
3955            else:
3956                e3 = roots[2].real()
3957                e2 = roots[1]
3958                e1 = roots[0]
3959               
3960        #curve of rank 0 doesn't need the whole algorithm
3961        if r == 0:
3962            int_points = set([])
3963            if disc > 0:
3964                #Points in egg have e2>=x>=e1
3965                search_points(ceil(e1), floor(e2)+1)
3966            #Points in noncompact component with x(P)< 2*max(|e1|,|e2|,|e3|) , espec. x(P)>=e3
3967            search_points(ceil(e3), ceil(2*max(abs(e1),abs(e2),abs(e3))))
3968            OP=self.point((0,1,0))
3969            for i in range(len_tors):
3970                P_tmp=tors_points[i]
3971                while P_tmp != OP:
3972                    if is_int(P_tmp.xy()[0]) and is_int(P_tmp.xy()[1]):
3973                            int_points.add(P_tmp)
3974                    P_tmp = P_tmp + tors_points[i]     
3975            return list(int_points)
3976        # Algorithm presented in [Co1]
3977        h_E = height_of_curve()
3978        w1 = self.period_lattice().basis()[0]
3979        w2 = self.period_lattice().basis()[1]
3980        mu = (log((abs(disc)),e)  + log_plus(j))/6 + log_plus(b2/12)
3981        if b2!=0:
3982            mu=mu+log(2,e)
3983       
3984        c1 = R(exp(mu + 2.14))
3985        c2 = min_eigenvalue(height_paaring_matrix(mw_base))
3986        c3 = (w1**2)*abs(b2)/48 + 8
3987        c5 = R(sqrt(c1*c3))
3988        c7 = abs((3*pi)/((w1**2) * (w1/w2).imag()))
3989       
3990        mw_base_log = [] #contains \Phi(Q_i)
3991        mod_h_list = []  #contains h_m(Q_i)
3992        c9_help_list = []
3993        for i in range(0,r):
3994            mw_base_log.append(abs(complex_elliptic_log(self,mw_base[i])))
3995            mod_h_list.append(modified_height(i))
3996            c9_help_list.append(sqrt(mod_h_list[i])/mw_base_log[i])
3997   
3998        c8 = max(e*h_E,max(mod_h_list))
3999        c9 = e/sqrt(c7) * min(c9_help_list)
4000        n=r+1
4001        c10 = R(2 * 10**(8+7*n) * R((2/e)**(2 * n**2)) * (n+1)**(4 * n**2 + 10 * n) * log(c9,e)**(-2*n - 1) * prod(mod_h_list))
4002   
4003        top = 128 #arbitrary first upper bound
4004        bottom = 0
4005        log_c9=log(c9,e); log_c5=log(c5,e)
4006        log_r_top = log(r*(1e1)**top, e)
4007   
4008        while R(c10*(log_r_top+log_c9)*(log(log_r_top,e)+h_E+log_c9)**(n+1)) > R(c2/2 * ((1e1)**top)**2 - log_c5):
4009            #initial bound 'top' too small, upshift of search interval
4010            bottom = top
4011            top = 2*top
4012        while top >= bottom: #binary-search like search for fitting exponent (bound)
4013            bound = floor(bottom + (top - bottom)/2)
4014            log_r_bound = log(r*(1e1)**bound, e)
4015            if R(c10*(log_r_bound+log_c9)*(log(log_r_bound, e)+h_E+log_c9)**(n+1)) > R(c2/2 * ((1e1)**bound)**2 - log_c5):
4016                bottom = bound + 1
4017            else:
4018                top = bound - 1
4019   
4020        H_q = R((1e1)**bound)
4021        break_cond = 0 #at least one reduction step
4022        #reduction via LLL
4023        while break_cond < 0.9: #as long as the improvement of the new bound in comparison to the old is greater than 10%
4024            c = R((H_q**n)*10)  #c has to be greater than H_q^n
4025            M = matrix.MatrixSpace(Z,n)
4026            m = M.identity_matrix()
4027            for i in range(r):
4028                m[i, r] = R(c*mw_base_log[i]).round()
4029            m[r,r] = R(c*w1).round()
4030   
4031            #LLL - implemented in sage - operates on rows not on columns
4032            m_LLL = m.LLL()
4033            m_gram = m_LLL.gram_schmidt()[0]
4034            b1_norm = R(m_LLL.row(0).norm())
4035   
4036            #compute constant c1 ~ c1_LLL of Corollary 2.3.17 and hence d(L,0)^2 ~ d_L_0
4037            c1_LLL = -1
4038            for i in range(n):
4039                tmp = R(b1_norm/(m_gram.row(i).norm()))
4040                if tmp > c1_LLL:
4041                    c1_LLL = tmp
4042   
4043            if c1_LLL < 0:
4044                raise RuntimeError, 'Unexpected intermediate result. Please try another Mordell-Weil base'
4045           
4046            d_L_0 = R(b1_norm**2 / c1_LLL)
4047   
4048            #Reducing of upper bound
4049            Q = r * H_q**2
4050            T = (1 + (3/2*r*H_q))/2 
4051            if d_L_0 < R(T**2+Q):
4052                d_L_0 = 10*(T**2*Q)
4053            low_bound = R((sqrt(d_L_0 - Q) - T)/c)
4054   
4055            #new bound according to low_bound and upper bound [c_5 exp((-c_2*H_q^2)/2)] provided by Corollary 8.7.3
4056            H_q_new = R(sqrt(log(low_bound/c5, e)/(-c2/2)))
4057            H_q_new = ceil(H_q_new)
4058            break_cond = R(H_q_new/H_q)
4059            H_q = H_q_new
4060            #END LLL-Reduction loop
4061   
4062   
4063        int_points = set([])
4064        ##for a more detailed view how the function works uncomment
4065        ##all print - statements below
4066        if disc > 0:
4067            ##Points in egg have e2>=x>=e1
4068            search_points(ceil(e1), floor(e2)+1)
4069            #print 'Points on compact component \n',list(int_points)
4070           
4071        ##Points in noncompact component with x(P)< 2*max(|e1|,|e2|,|e3|) , espec. x(P)>=e3
4072        search_points(ceil(e3), ceil(2*max(abs(e1),abs(e2),abs(e3))))
4073        #print 'Points on compact component and non-compact component part 1 \n', list(int_points)
4074       
4075        #print 'starting search of remainig points using bound ',H_q
4076        search_remaining_points()
4077        #print 'Integral points:\n',list(int_points)
4078        #print 'Total amount of points:',len(int_points)
4079   
4080        return int_points
4081
4082
4083def cremona_curves(conductors):
4084    """
4085    Return iterator over all known curves (in database) with conductor
4086    in the list of conductors.
4087
4088    EXAMPLES:
4089        sage: [(E.label(), E.rank()) for E in cremona_curves(srange(35,40))]
4090        [('35a1', 0),
4091        ('35a2', 0),
4092        ('35a3', 0),
4093        ('36a1', 0),
4094        ('36a2', 0),
4095        ('36a3', 0),
4096        ('36a4', 0),
4097        ('37a1', 1),
4098        ('37b1', 0),
4099        ('37b2', 0),
4100        ('37b3', 0),
4101        ('38a1', 0),
4102        ('38a2', 0),
4103        ('38a3', 0),
4104        ('38b1', 0),
4105        ('38b2', 0),
4106        ('39a1', 0),
4107        ('39a2', 0),
4108        ('39a3', 0),
4109        ('39a4', 0)]
4110    """
4111    if isinstance(conductors, (int,long, rings.RingElement)):
4112        conductors = [conductors]
4113    return sage.databases.cremona.CremonaDatabase().iter(conductors)
4114
4115def cremona_optimal_curves(conductors):
4116    """
4117    Return iterator over all known optimal curves (in database) with
4118    conductor in the list of conductors.
4119
4120    EXAMPLES:
4121        sage: [(E.label(), E.rank()) for E in cremona_optimal_curves(srange(35,40))]
4122        [('35a1', 0),
4123        ('36a1', 0),
4124        ('37a1', 1),
4125        ('37b1', 0),
4126        ('38a1', 0),
4127        ('38b1', 0),
4128        ('39a1', 0)]
4129    """
4130    if isinstance(conductors, (int,long,rings.RingElement)):
4131        conductors = [conductors]
4132    return sage.databases.cremona.CremonaDatabase().iter_optimal(conductors)
4133