Ticket #3674: ell_rational_field.py

"""
Elliptic curves over the rational numbers

AUTHORS:
   -- William Stein (2005): first version
   -- William Stein (2006-02-26): fixed Lseries_extended which didn't work
            because of changes elsewhere in SAGE.
   -- David Harvey (2006-09): Added padic_E2, padic_sigma, padic_height,
            padic_regulator methods.
   -- David Harvey (2007-02): reworked padic-height related code
   -- Christian Wuthrich (2007): added padic sha computation
   -- David Roe (2007-9): moved sha, l-series and p-adic functionality to separate files.
   -- John Cremona (2008-01)
   -- Tobias Nagel & Michael Mardaus (2008-07): added integral_points
"""

#*****************************************************************************
#       Copyright (C) 2005,2006,2007 William Stein <wstein@gmail.com>
#
#  Distributed under the terms of the GNU General Public License (GPL)
#
#    This code is distributed in the hope that it will be useful,
#    but WITHOUT ANY WARRANTY; without even the implied warranty of
#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
#    General Public License for more details.
#
#  The full text of the GPL is available at:
#
#                  http://www.gnu.org/licenses/
#*****************************************************************************

import ell_point
import formal_group
import rational_torsion
from ell_generic import EllipticCurve_generic
from ell_number_field import EllipticCurve_number_field

import sage.groups.all
import sage.rings.arith as arith
import sage.rings.all as rings
import sage.rings.number_field.number_field as number_field
import sage.misc.misc as misc
from sage.misc.all import verbose 
import sage.functions.constants as constants
import sage.modular.modform.constructor
import sage.modular.modform.element
from sage.misc.functional import log
from sage.rings.padics.factory import Zp, Qp

# Use some interval arithmetic to guarantee correctness.  We assume
# that alpha is computed to the precision of a float.
# IR = rings.RIF
#from sage.rings.interval import IntervalRing; IR = IntervalRing()

import sage.matrix.all as matrix
import sage.databases.cremona
from   sage.libs.pari.all import pari
import sage.functions.transcendental as transcendental
import math
from sage.calculus.calculus import sqrt, arcsin, floor, ceil
import sage.libs.mwrank.all as mwrank
import constructor
from sage.interfaces.all import gp

import ell_modular_symbols
import padic_lseries
import padics

from lseries_ell import Lseries_ell

import mod5family

from sage.rings.all import (
    PowerSeriesRing, LaurentSeriesRing, O, 
    infinity as oo,
    Integer,
    Integers,
    IntegerRing, RealField,
    ComplexField, RationalField)

import gp_cremona
import sea

from gp_simon import simon_two_descent

import ell_tate_curve

factor = arith.factor
exp = math.exp
mul = misc.mul
next_prime = arith.next_prime
kronecker_symbol = arith.kronecker_symbol

Q = RationalField()         
C = ComplexField()
R = RealField()
Z = IntegerRing()
IR = rings.RealIntervalField(20)

_MAX_HEIGHT=21

# CMJ is a dict of pairs (j,D) where j is a rational CM j-invariant
# and D is the corresponding quadratic discriminant

CMJ={ 0: -3, 54000: -12, -12288000: -27, 1728: -4, 287496: -16, 
      -3375: -7, 16581375: -28, 8000: -8, -32768: -11, -884736: -19, 
      -884736000: -43, -147197952000: -67, -262537412640768000: -163}



class EllipticCurve_rational_field(EllipticCurve_number_field):
    """
    Elliptic curve over the Rational Field.
    """
    def __init__(self, ainvs, extra=None):
        if extra != None:   # possibility of two arguments (the first would be the field)
            ainvs = extra
        if isinstance(ainvs, str):
            label = ainvs
            X = sage.databases.cremona.CremonaDatabase()[label]
            EllipticCurve_number_field.__init__(self, Q, X.a_invariants())
            for attr in ['rank', 'torsion_order', 'cremona_label', 'conductor',
                         'modular_degree', 'gens', 'regulator']:
                s = "_EllipticCurve_rational_field__"+attr
                if hasattr(X,s):
                    setattr(self, s, getattr(X, s))
            return
        EllipticCurve_number_field.__init__(self, Q, ainvs)
        self.__np = {}
        self.__gens = {}
        self.__rank = {}
        self.__regulator = {}
        if self.base_ring() != Q:
            raise TypeError, "Base field (=%s) must be the Rational Field."%self.base_ring()
        
    def _set_rank(self, r):
        """
        Internal function to set the cached rank of this elliptic curve to r.

        WARNING:  No checking is done!  Not intended for use by users.

        EXAMPLES:
            sage: E=EllipticCurve('37a1')
            sage: E._set_rank(99)  # bogus value -- not checked
            sage: E.rank()         # returns bogus cached value
            99
            sage: E.gens()         # causes actual rank to be computed
            [(0 : -1 : 1)]
            sage: E.rank()         # the correct rank
            1
        """
        self.__rank = {}
        self.__rank[True] = Integer(r)

    def _set_torsion_order(self, t):
        """
        Internal function to set the cached torsion order of this elliptic curve to t.

        WARNING:  No checking is done!  Not intended for use by users.

        EXAMPLES:
            sage: E=EllipticCurve('37a1')
            sage: E._set_torsion_order(99)  # bogus value -- not checked
            sage: E.torsion_order()         # returns bogus cached value
            99
            sage: T = E.torsion_subgroup()  # causes actual torsion to be computed
            sage: E.torsion_order()         # the correct value
            1
        """
        self.__torsion_order = Integer(t)

    def _set_cremona_label(self, L):
        """
        Internal function to set the cached label of this elliptic curve to L.

        WARNING:  No checking is done!  Not intended for use by users.

        EXAMPLES:
            sage: E=EllipticCurve('37a1')
            sage: E._set_cremona_label('bogus')
            sage: E.label()
            'bogus'
            sage: E.database_curve().label()
            '37a1'
            sage: E.label() # no change
            'bogus'
            sage: E._set_cremona_label(E.database_curve().label())
            sage: E.label() # now it is correct
            '37a1'
        """
        self.__cremona_label = L

    def _set_conductor(self, N):
        """
        Internal function to set the cached conductor of this elliptic curve to N.

        WARNING: No checking is done!  Not intended for use by users.
        Setting to the wrong value will cause strange problems (see
        examples).
        
        EXAMPLES:
            sage: E=EllipticCurve('37a1')
            sage: E._set_conductor(99)      # bogus value -- not checked
            sage: E.conductor()             # returns bogus cached value
            99

        This will not work since the conductor is used when searching
        the database:
            sage: E._set_conductor(E.database_curve().conductor()) 
            Traceback (most recent call last): 
            ...  
            RuntimeError: Elliptic curve ... not in the database.
            sage: E._set_conductor(EllipticCurve(E.a_invariants()).database_curve().conductor()) 
            sage: E.conductor()             # returns correct value
            37
        """
        self.__conductor_pari = Integer(N)

    def _set_modular_degree(self, deg):
        """
        Internal function to set the cached modular degree of this elliptic curve to deg.

        WARNING:  No checking is done!

        EXAMPLES:
            sage: E=EllipticCurve('5077a1')
            sage: E.modular_degree()
            1984
            sage: E._set_modular_degree(123456789)
            sage: E.modular_degree()
            123456789
            sage: E._set_modular_degree(E.database_curve().modular_degree())
            sage: E.modular_degree()
            1984
        """
        self.__modular_degree = Integer(deg)
        
    def _set_gens(self, gens):
        """
        Internal function to set the cached generators of this elliptic curve to gens.

        WARNING:  No checking is done!

        EXAMPLES:
            sage: E=EllipticCurve('5077a1')
            sage: E.rank()
            3
            sage: E.gens() # random
            [(-2 : 3 : 1), (-7/4 : 25/8 : 1), (1 : -1 : 1)]
            sage: E._set_gens([]) # bogus list
            sage: E.rank()        # unchanged 
            3
            sage: E._set_gens(E.database_curve().gens())
            sage: E.gens()
            [(-2 : 3 : 1), (-7/4 : 25/8 : 1), (1 : -1 : 1)]
        """
        self.__gens = {}
        self.__gens[True] = [self.point(x, check=True) for x in gens]
        self.__gens[True].sort()

    def is_integral(self):
        """
        Returns True if this elliptic curve has integral coefficients (in Z)

        EXAMPLES:
            sage: E=EllipticCurve(QQ,[1,1]); E
            Elliptic Curve defined by y^2  = x^3 + x +1 over Rational Field
            sage: E.is_integral()
            True
            sage: E2=E.change_weierstrass_model(2,0,0,0); E2
            Elliptic Curve defined by y^2  = x^3 + 1/16*x + 1/64 over Rational Field
            sage: E2.is_integral()
            False
         """
        try:
            return self.__is_integral
        except AttributeError:
            one = Integer(1)
            self.__is_integral = bool(misc.mul([x.denominator() == 1 for x in self.ainvs()]))
            return self.__is_integral
            

    def mwrank(self, options=''):
        """
        Run Cremona's mwrank program on this elliptic curve and
        return the result as a string.

        INPUT:
            options -- string; passed when starting mwrank.  The format is
        q p<precision> v<verbosity> b<hlim_q> x<naux>  c<hlim_c> l t o s d>]

        OUTPUT:
            string -- output of mwrank on this curve

        NOTE: The output is a raw string and completely illegible
            using automatic display, so it is recommended to use print
            for legible outout

        EXAMPLES:
            sage: E = EllipticCurve('37a1')
            sage: E.mwrank() #random
            ...
            sage: print E.mwrank()
            Curve [0,0,1,-1,0] :        Basic pair: I=48, J=-432
            disc=255744
            ...
            Generator 1 is [0:-1:1]; height 0.05111...
            <BLANKLINE>
            Regulator = 0.05111...
            <BLANKLINE>
            The rank and full Mordell-Weil basis have been determined unconditionally.
            ...

        Options to mwrank can be passed:
            sage: E = EllipticCurve([0,0,0,877,0])

        Run mwrank with 'verbose' flag set to 0 but list generators if found
            sage: print E.mwrank('-v0 -l')
            Curve [0,0,0,877,0] :   0 <= rank <= 1
            Regulator = 1
            
        Run mwrank again, this time with a higher bound for point
        searching on homogeneous spaces:

            sage: print E.mwrank('-v0 -l -b11')
            Curve [0,0,0,877,0] :   Rank = 1
            Generator 1 is [29604565304828237474403861024284371796799791624792913256602210:-256256267988926809388776834045513089648669153204356603464786949:490078023219787588959802933995928925096061616470779979261000]; height 95.980371987964
            Regulator = 95.980371987964
        """
        if options == "":
            from sage.interfaces.all import mwrank
        else:
            from sage.interfaces.all import Mwrank
            mwrank = Mwrank(options=options)
        return mwrank(self.a_invariants())

    def conductor(self, algorithm="pari"):
        """
        Returns the conductor of the elliptic curve.

        INPUT:
            algorithm -- str, (default: "pari")
                   "pari"   -- use the PARI C-library ellglobalred
                               implementation of Tate's algorithm
                   "mwrank" -- use Cremona's mwrank implementation of
                               Tate's algorithm; can be faster if the
                               curve has integer coefficients (TODO:
                               limited to small conductor until mwrank
                               gets integer factorization)
                   "gp" -- use the GP interpreter.
                   "all" -- use both implementations, verify that the
                            results are the same (or raise an error),
                            and output the common value.
                                     
        EXAMPLE:
            sage: E = EllipticCurve([1, -1, 1, -29372, -1932937])
            sage: E.conductor(algorithm="pari")
            3006
            sage: E.conductor(algorithm="mwrank")
            3006
            sage: E.conductor(algorithm="gp")
            3006
            sage: E.conductor(algorithm="all")
            3006

        NOTE: The conductor computed using each algorithm is cached separately.
        Thus calling E.conductor("pari"), then E.conductor("mwrank") and
        getting the same result checks that both systems compute the same answer.
        """

        if algorithm == "pari":
            try:
                return self.__conductor_pari
            except AttributeError:
                self.__conductor_pari = Integer(self.pari_mincurve().ellglobalred()[0])
            return self.__conductor_pari

        elif algorithm == "gp":
            try:
                return self.__conductor_gp
            except AttributeError:
                self.__conductor_gp = Integer(gp.eval('ellglobalred(ellinit(%s,0))[1]'%self.a_invariants()))
                return self.__conductor_gp

        elif algorithm == "mwrank":
            try:
                return self.__conductor_mwrank
            except AttributeError:
                if self.is_integral():
                    self.__conductor_mwrank = Integer(self.mwrank_curve().conductor())
                else:
                    self.__conductor_mwrank = Integer(self.minimal_model().mwrank_curve().conductor())
            return self.__conductor_mwrank

        elif algorithm == "all":
            N1 = self.conductor("pari")
            N2 = self.conductor("mwrank")
            N3 = self.conductor("gp")
            if N1 != N2 or N2 != N3:
                raise ArithmeticError, "Pari, mwrank and gp compute different conductors (%s,%s,%s) for %s"%(
                    N1, N2, N3, self)
            return N1
        else:
            raise RuntimeError, "algorithm '%s' is not known."%algorithm

    ####################################################################
    #  Access to PARI curves related to this curve.
    ####################################################################

    def pari_curve(self, prec = None, factor = 1):
        """
        Return the PARI curve corresponding to this elliptic curve.

        INPUT:
            prec -- The precision of quantities calculated for the
                    returned curve (in decimal digits).  if None, defaults
                    to factor * the precision of the largest cached curve
                    (or 10 if none yet computed)
            factor -- the factor to increase the precision over the
                      maximum previously computed precision.  Only used if
                      prec (which gives an explicit precision) is None.
                
        EXAMPLES:
            sage: E = EllipticCurve([0, 0,1,-1,0])
            sage: e = E.pari_curve()
            sage: type(e)
            <type 'sage.libs.pari.gen.gen'>
            sage: e.type()
            't_VEC'
            sage: e.ellan(10)
            [1, -2, -3, 2, -2, 6, -1, 0, 6, 4]

            sage: E = EllipticCurve(RationalField(), ['1/3', '2/3'])
            sage: e = E.pari_curve()
            sage: e.type()
            't_VEC'
            sage: e[:5]
            [0, 0, 0, 1/3, 2/3]
        """
        if prec is None:
            try:
                L = self._pari_curve.keys()
                L.sort()
                if factor == 1:
                    return self._pari_curve[L[-1]]
                else:
                    prec = int(factor * L[-1])
                    self._pari_curve[prec] = pari(self.a_invariants()).ellinit(precision=prec)
                    return self._pari_curve[prec]
            except AttributeError:
                pass
        try:
            return self._pari_curve[prec]
        except AttributeError:
            prec = 10
            self._pari_curve = {}
        except KeyError:
            pass
        self._pari_curve[prec] = pari(self.a_invariants()).ellinit(precision=prec)
        return self._pari_curve[prec]

    def pari_mincurve(self, prec = None):
        """
        Return the PARI curve corresponding to a minimal model
        for this elliptic curve.

        INPUT:
        prec -- The precision of quantities calculated for the
                returned curve (in decimal digits).  if None, defaults
                to the precision of the largest cached curve (or 10 if
                none yet computed)

        EXAMPLES:
            sage: E = EllipticCurve(RationalField(), ['1/3', '2/3'])
            sage: e = E.pari_mincurve()
            sage: e[:5]
            [0, 0, 0, 27, 486]
            sage: E.conductor()
            47232
            sage: e.ellglobalred()
            [47232, [1, 0, 0, 0], 2]
        """
        if prec is None:
            try:
                L = self._pari_mincurve.keys()
                L.sort()
                return self._pari_mincurve[L[len(L) - 1]]
            except AttributeError:
                pass
        try:
            return self._pari_mincurve[prec]
        except AttributeError:
            prec = 10
            self._pari_mincurve = {}
        except KeyError:
            pass
        e = self.pari_curve(prec)
        mc, change = e.ellminimalmodel()
        self._pari_mincurve[prec] = mc
        # self.__min_transform = change
        return mc

    def database_curve(self):
        """
        Return the curve in the elliptic curve database isomorphic to
        this curve, if possible.  Otherwise raise a RuntimeError
        exception.  

        EXAMPLES:
            sage: E = EllipticCurve([0,1,2,3,4])
            sage: E.database_curve()
            Elliptic Curve defined by y^2  = x^3 + x^2 + 3*x + 5 over Rational Field

        NOTES: The model of the curve in the database can be different
               than the Weierstrass model for this curve, e.g.,
               database models are always minimal.
        """
        try:
            return self.__database_curve
        except AttributeError:
            misc.verbose("Looking up %s in the database."%self)
            D = sage.databases.cremona.CremonaDatabase()
            ainvs = self.minimal_model().ainvs()
            try:
                self.__database_curve = D.elliptic_curve_from_ainvs(self.conductor(), ainvs)
            except RuntimeError:
                raise RuntimeError, "Elliptic curve %s not in the database."%self
            return self.__database_curve
            

    def Np(self, p):
        """
        The number of points on E modulo p, where p is a prime, not
        necessarily of good reduction.  (When p is a bad prime, also
        counts the singular point.)

        EXAMPLES:
            sage: E = EllipticCurve([0, -1, 1, -10, -20])
            sage: E.Np(2)
            5
            sage: E.Np(3)
            5
            sage: E.conductor()
            11
            sage: E.Np(11)
            11
        """
        if self.conductor() % p == 0:
            return p + 1 - self.ap(p)
        #raise ArithmeticError, "p (=%s) must be a prime of good reduction"%p
        if p < 1125899906842624:   # TODO: choose more wisely?
            return p+1 - self.ap(p)
        else:
            return self.sea(p)

    def sea(self, p, early_abort=False):
        r"""
        Return the number of points on $E$ over $\F_p$ computed using
        the SEA algorithm, as implemented in PARI by Christophe Doche
        and Sylvain Duquesne.

        INPUT:
            p -- a prime number
            early_abort -- bool (default: Falst); if True an early abort technique
                       is used and the computation is interrupted as soon
                       as a small divisor of the order is detected.

        \note{As of 2006-02-02 this function does not work on
        Microsoft Windows under Cygwin (though it works under
        vmware of course).}

        EXAMPLES:
            sage: E = EllipticCurve('37a')
            sage: E.sea(next_prime(10^30))
            1000000000000001426441464441649
        """
        p = rings.Integer(p)
        return sea.ellsea(self.minimal_model().a_invariants(), p, early_abort=early_abort)

    #def __pari_double_prec(self):
    #    EllipticCurve_number_field._EllipticCurve__pari_double_prec(self)
    #    try:
    #        del self._pari_mincurve
    #    except AttributeError:
    #        pass
        
    ####################################################################
    #  Access to mwrank
    ####################################################################
    def mwrank_curve(self, verbose=False):
        """
        Construct an mwrank_EllipticCurve from this elliptic curve

        The resulting mwrank_EllipticCurve has available methods from
        John Cremona's eclib library.

        EXAMPLES:
            sage: E=EllipticCurve('11a1')
            sage: EE=E.mwrank_curve()
            sage: EE
            y^2+ y = x^3 - x^2 - 10*x - 20
            sage: type(EE)
            <class 'sage.libs.mwrank.interface.mwrank_EllipticCurve'>
            sage: EE.isogeny_class()
            ([[0, -1, 1, -10, -20], [0, -1, 1, -7820, -263580], [0, -1, 1, 0, 0]],
            [[0, 5, 5], [5, 0, 0], [5, 0, 0]])
        """
        try:
            return self.__mwrank_curve
        except AttributeError:
            pass
        self.__mwrank_curve = mwrank.mwrank_EllipticCurve(
            self.ainvs(), verbose=verbose)
        return self.__mwrank_curve

    def two_descent(self, verbose=True,
                    selmer_only = False,
                    first_limit = 20,
                    second_limit = 8,
                    n_aux = -1,
                    second_descent = 1):
        """
        Compute 2-descent data for this curve.

        INPUT:
            verbose     -- (default: True) print what mwrank is doing
                           If False, *no output* is 
            selmer_only -- (default: False) selmer_only switch
            first_limit -- (default: 20) firstlim is bound on |x|+|z|
            second_limit-- (default: 8)  secondlim is bound on log max {|x|,|z| },
                                         i.e. logarithmic
            n_aux       -- (default: -1) n_aux only relevant for general
                           2-descent when 2-torsion trivial; n_aux=-1 causes default
                           to be used (depends on method)
            second_descent -- (default: True) second_descent only relevant for
                           descent via 2-isogeny
        OUTPUT:

            Nothing -- nothing is returned (though much is printed
                                            unless verbose=False)

        EXAMPLES:
            sage: E=EllipticCurve('37a1')
            sage: E.two_descent(verbose=False) # no output
        """
        self.mwrank_curve().two_descent(verbose, selmer_only,
                                        first_limit, second_limit,
                                        n_aux, second_descent)


    ####################################################################
    #  Etc.
    ####################################################################

    def aplist(self, n, python_ints=False):
        r"""
        The Fourier coefficients $a_p$ of the modular form attached to
        this elliptic curve, for all primes $p\leq n$.

        INPUT:
            n -- integer
            python_ints -- bool (default: False); if True return a list of
                      Python ints instead of SAGE integers.

        OUTPUT:
            -- list of integers

        EXAMPLES:
            sage: e = EllipticCurve('37a')
            sage: e.aplist(1)
            []
            sage: e.aplist(2)
            [-2]
            sage: e.aplist(10)
            [-2, -3, -2, -1]
            sage: v = e.aplist(13); v
            [-2, -3, -2, -1, -5, -2]
            sage: type(v[0])
            <type 'sage.rings.integer.Integer'>
            sage: type(e.aplist(13, python_ints=True)[0])
            <type 'int'>
        """
        # How is this result dependant on the real precision in pari?  At all?
        e = self.pari_mincurve()
        v = e.ellaplist(n, python_ints=True)
        if python_ints:
            return v
        else:
            return [Integer(a) for a in v]
        
        

    def anlist(self, n, python_ints=False):
        """
        The Fourier coefficients up to and including $a_n$ of the
        modular form attached to this elliptic curve.  The ith element
        of the return list is a[i].

        INPUT:
            n -- integer
            python_ints -- bool (default: False); if True return a list of
                      Python ints instead of SAGE integers.

        OUTPUT:
            -- list of integers

        EXAMPLES:
            sage: E = EllipticCurve([0, -1, 1, -10, -20])
            sage: E.anlist(3)
            [0, 1, -2, -1]
            
            sage: E = EllipticCurve([0,1])
            sage: E.anlist(20)
            [0, 1, 0, 0, 0, 0, 0, -4, 0, 0, 0, 0, 0, 2, 0, 0, 0, 0, 0, 8, 0]
        """
        n = int(n)
        # How is this result dependent on the real precision in Pari?  At all?
        e = self.pari_mincurve()
        if n >= 2147483648:
            raise RuntimeError, "anlist: n (=%s) must be < 2147483648."%n

        v = [0] + e.ellan(n, python_ints=True)
        if not python_ints:
            v = [Integer(x) for x in v]
        return v
        
        
        # There is some overheard associated with coercing the PARI
        # list back to Python, but it's not bad.  It's better to do it
        # this way instead of trying to eval the whole list, since the
        # int conversion is done very sensibly.  NOTE: This would fail
        # if a_n won't fit in a C int, i.e., is bigger than
        # 2147483648; however, we wouldn't realistically compute
        # anlist for n that large anyways.
        #
        # Some relevant timings:
        #
        # E <--> [0, 1, 1, -2, 0]   389A
        #  E = EllipticCurve([0, 1, 1, -2, 0]);   // SAGE or MAGMA
        #  e = E.pari_mincurve()
        #  f = ellinit([0,1,1,-2,0]);
        #
        #  Computation                                              Time (1.6Ghz Pentium-4m laptop)
        #  time v:=TracesOfFrobenius(E,10000);  // MAGMA            0.120
        #  gettime;v=ellan(f,10000);gettime/1000                    0.046
        #  time v=e.ellan (10000)                                   0.04
        #  time v=E.anlist(10000)                                   0.07
        
        #  time v:=TracesOfFrobenius(E,100000);  // MAGMA           1.620
        #  gettime;v=ellan(f,100000);gettime/1000                   0.676
        #  time v=e.ellan (100000)                                  0.7
        #  time v=E.anlist(100000)                                  0.83

        #  time v:=TracesOfFrobenius(E,1000000);  // MAGMA          20.850
        #  gettime;v=ellan(f,1000000);gettime/1000                  9.238
        #  time v=e.ellan (1000000)                                 9.61
        #  time v=E.anlist(1000000)                                 10.95  (13.171 in cygwin vmware)
        
        #  time v:=TracesOfFrobenius(E,10000000);  //MAGMA          257.850
        #  gettime;v=ellan(f,10000000);gettime/1000      FAILS no matter how many allocatemem()'s!!
        #  time v=e.ellan (10000000)                                139.37  
        #  time v=E.anlist(10000000)                                136.32
        #
        #  The last SAGE comp retries with stack size 40MB,
        #  80MB, 160MB, and succeeds last time.  It's very interesting that this
        #  last computation is *not* possible in GP, but works in py_pari!
        #

    def q_expansion(self, prec):
        r""" 
        Return the $q$-expansion to precision prec of the newform