Ticket #12015: trac_12015-part1.patch

File trac_12015-part1.patch, 101.7 KB (added by William Stein, 11 years ago)

first patch -- this implements all functionality for first version and passes all its tests. Unfortunately, the doctest coverage score is 36%. To finish this ticket, I have to get it to 100% and deal with the bugs that getting the coverage up reveals.

  • module_list.py

    # HG changeset patch
    # User William Stein <wstein@gmail.com>
    # Date 1321129751 28800
    # Node ID 6522f45f8c8b377f87d751c3d1e5a9538cca9f5e
    # Parent  76a3c98cb70a82b14e3f5ba5aac899efcb77465a
    trac 12015: L-series attached to general Euler products
    
    diff --git a/module_list.py b/module_list.py
    a b  
    532532
    533533    ################################
    534534    ##
     535    ## sage.lfunctions
     536    ##
     537    ################################
     538
     539    Extension("sage.lfunctions.eulerprod_fast",
     540              ["sage/lfunctions/eulerprod_fast.pyx"]),
     541
     542    ################################
     543    ##
    535544    ## sage.libs
    536545    ##
    537546    ################################
  • new file sage/lfunctions/eulerprod.py

    diff --git a/sage/lfunctions/eulerprod.py b/sage/lfunctions/eulerprod.py
    new file mode 100644
    - +  
     1"""
     2General L-series
     3
     4AUTHOR:
     5    - William Stein
     6
     7TODO:
     8    - Symmetric powers (and modular degree -- see trac 9758)
     9    - Triple product L-functions: Gross-Kudla, Zhang, etc -- see the code in triple_prod/triple.py
     10    - Support L-calc L-function
     11    - Make it so we use exactly one GP session for *all* of the Dokchitser L-functions
     12    - Tensor products
     13    - Genus 2 curves, via smalljac and genus2reduction
     14    - Fast L-series of elliptic curves over number fields (not just sqrt(5)), via smalljac
     15    - Inverse of number_of_coefficients function.
     16"""
     17
     18#################################################################################
     19#
     20# (c) Copyright 2011 William Stein
     21#
     22#  This file is part of SAGE
     23#
     24#  SAGE is free software: you can redistribute it and/or modify
     25#  it under the terms of the GNU General Public License as published by
     26#  the Free Software Foundation, either version 2 of the License, or
     27#  (at your option) any later version.
     28#
     29#  SAGE is distributed in the hope that it will be useful,
     30#  but WITHOUT ANY WARRANTY; without even the implied warranty of
     31#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
     32#  GNU General Public License for more details.
     33#
     34#  You should have received a copy of the GNU General Public License
     35#  along with this program.  If not, see <http://www.gnu.org/licenses/>.
     36#
     37#################################################################################
     38
     39import copy, math, types
     40
     41from sage.all import prime_range, cached_method, sqrt, SR, vector
     42from sage.rings.all import is_RationalField, ZZ, Integer, QQ, O, ComplexField, CDF, primes, infinity as oo
     43from sage.schemes.elliptic_curves.ell_generic import is_EllipticCurve
     44from sage.misc.all import prod
     45from sage.modular.abvar.abvar import is_ModularAbelianVariety
     46from sage.modular.dirichlet import is_DirichletCharacter
     47from sage.lfunctions.dokchitser import Dokchitser
     48from sage.modular.modsym.space import is_ModularSymbolsSpace
     49from sage.modular.abvar.abvar import is_ModularAbelianVariety
     50from sage.rings.number_field.number_field_base import is_NumberField
     51import sage.modular.modform.element
     52from sage.modular.all import Newform
     53from sage.structure.factorization import Factorization
     54from sage.misc.mrange import cartesian_product_iterator
     55
     56from sage.lfunctions.eulerprod_fast import extend_multiplicatively_generic
     57
     58I = sqrt(-1)
     59
     60def prec(s):
     61    """
     62    Return precision of s, if it has a precision attribute.  Otherwise
     63    return 53.  This is useful internally in this module.
     64
     65    EXAMPLES::
     66
     67        sage: from sage.lfunctions.eulerprod import prec
     68        sage: prec(ComplexField(100)(1))
     69        100
     70        sage: prec(RealField(125)(1))
     71        125
     72        sage: prec(1/3)
     73        53   
     74    """
     75    if hasattr(s, 'prec'):
     76        return s.prec()
     77    return 53
     78
     79def norm(a):
     80    """
     81    Return the norm of a, for a in either a number field or QQ.
     82
     83    This is a function used internally in this module, mainly because
     84    elements of QQ and ZZ have no norm method.
     85
     86    EXAMPLES::
     87
     88        sage: from sage.lfunctions.eulerprod import norm
     89        sage: K.<a> = NumberField(x^2-x-1)
     90        sage: norm(a+5)
     91        29
     92        sage: (a+5).norm()
     93        29
     94        sage: norm(17)
     95        17
     96    """
     97    try:
     98        return a.norm()
     99    except AttributeError:
     100        return a
     101
     102def tiny(prec):
     103    """
     104    Return a number that we consider tiny to the given precision prec
     105    in bits.   This is used in various places as "zero" to the given
     106    precision for various checks, e.g., of correctness of the
     107    functional equation.
     108    """
     109    return max(1e-8, 1.0/2**(prec-1))
     110
     111def prime_below(P):
     112    """
     113    Return the prime in ZZ below the prime P (or element of QQ).
     114   
     115    EXAMPLES::
     116
     117        sage: from sage.lfunctions.eulerprod import prime_below
     118        sage: K.<a> = NumberField(x^2-x-1)
     119        sage: prime_below(K.prime_above(11))
     120        11
     121        sage: prime_below(K.prime_above(5))
     122        5
     123        sage: prime_below(K.prime_above(3))
     124        3
     125        sage: prime_below(7)
     126        7
     127    """
     128    try:
     129        return P.smallest_integer()
     130    except AttributeError:
     131        return ZZ(P)
     132
     133
     134class LSeriesDerivative(object):
     135    """
     136    The formal derivative of an L-series.
     137
     138    EXAMPLES::
     139
     140        sage: from sage.lfunctions.eulerprod import LSeries
     141        sage: L = LSeries('delta')
     142        sage: L.derivative()
     143        First derivative of L-function associated to Ramanujan's Delta (a weight 12 cusp form)
     144        sage: L.derivative()(11/2)
     145        0.125386233743526
     146
     147    We directly create an instance of the class (users shouldn't need to do this)::
     148
     149        sage: from sage.lfunctions.eulerprod import LSeriesDerivative
     150        sage: Ld = LSeriesDerivative(L, 2); Ld
     151        Second derivative of L-function associated to Ramanujan's Delta (a weight 12 cusp form)
     152        sage: type(Ld)
     153        <class 'sage.lfunctions.eulerprod.LSeriesDerivative'>
     154    """
     155    def __init__(self, lseries, k):
     156        """
     157        INPUT:
     158            - lseries -- any LSeries object (derives from LseriesAbstract)
     159            - k -- positive integer
     160        """
     161        k = ZZ(k)
     162        if k <= 0:
     163            raise ValueError, "k must be a positive integer"
     164        self._lseries = lseries
     165        self._k = k
     166
     167    def __cmp__(self, right):
     168        return cmp((self._lseries,self._k), (right._lseries, right._k))
     169
     170    def __call__(self, s):
     171        """
     172        Return the value of this derivative at s, which must coerce to a
     173        complex number.  The computation is to the same precision as s,
     174        or to 53 bits of precision if s is exact.
     175
     176        As usual, if s has large imaginary part, then there could be
     177        substantial precision loss (a warning is printed in that case).
     178
     179        EXAMPLES::
     180
     181            sage: from sage.lfunctions.eulerprod import LSeries
     182            sage: L = LSeries(EllipticCurve('389a'))
     183            sage: L1 = L.derivative(1); L1(1)
     184            -1.94715429754927e-20
     185            sage: L1(2)
     186            0.436337613850735
     187            sage: L1(I)
     188            -19.8890471908356 + 31.2633280771869*I
     189            sage: L2 = L.derivative(2); L2(1)
     190            1.51863300057685
     191            sage: L2(I)
     192            134.536162459604 - 62.6542402272310*I
     193        """
     194        return self._lseries._function(prec(s)).derivative(s, self._k)
     195
     196    def __repr__(self):
     197        """
     198        String representation of this derivative of an L-series.
     199
     200        EXAMPLES::
     201
     202            sage: from sage.lfunctions.eulerprod import LSeries
     203            sage: L = LSeries('delta')
     204            sage: L.derivative(1).__repr__()
     205            "First derivative of L-function associated to Ramanujan's Delta (a weight 12 cusp form)"
     206            sage: L.derivative(2).__repr__()
     207            "Second derivative of L-function associated to Ramanujan's Delta (a weight 12 cusp form)"
     208            sage: L.derivative(3).__repr__()
     209            "Third derivative of L-function associated to Ramanujan's Delta (a weight 12 cusp form)"
     210            sage: L.derivative(4).__repr__()
     211            "4-th derivative of L-function associated to Ramanujan's Delta (a weight 12 cusp form)"
     212            sage: L.derivative(2011).__repr__()
     213            "2011-th derivative of L-function associated to Ramanujan's Delta (a weight 12 cusp form)"       
     214        """
     215        k = self._k
     216        if k == 1:
     217            kth = 'First'
     218        elif k == 2:
     219            kth = 'Second'
     220        elif k == 3:
     221            kth = 'Third'
     222        else:
     223            kth = '%s-th'%k
     224        return "%s derivative of %s"%(kth, self._lseries)
     225
     226    def derivative(self, k=1):
     227        """
     228        Return the k-th derivative of this derivative object.
     229       
     230        EXAMPLES::
     231
     232            sage: from sage.lfunctions.eulerprod import LSeries
     233            sage: f = Newforms(43,2,names='a')[1]; f
     234            q + a1*q^2 - a1*q^3 + (-a1 + 2)*q^5 + O(q^6)
     235            sage: L = LSeries(f); L1 = L.derivative()
     236            sage: L1(1)
     237            0.331674007376949
     238            sage: L(1)
     239            0.620539857407845
     240            sage: L = LSeries(f); L1 = L.derivative(); L1
     241            First derivative of L-series of a degree 2 newform of level 43 and weight 2
     242            sage: L1.derivative()
     243            Second derivative of L-series of a degree 2 newform of level 43 and weight 2
     244            sage: L1.derivative(3)
     245            4-th derivative of L-series of a degree 2 newform of level 43 and weight 2
     246       
     247        """
     248        if k == 0:
     249            return self
     250        return LSeriesDerivative(self._lseries, self._k + k)
     251
     252class LSeriesParentClass(object):
     253    def __contains__(self, x):
     254        return isinstance(x, (LSeriesAbstract, LSeriesProduct))
     255       
     256    def __call__(self, x):
     257        """
     258        EXAMPLES::
     259
     260            sage: from sage.lfunctions.eulerprod import LSeries
     261            sage: L = LSeries('zeta')
     262            sage: P = L.parent(); P
     263            All L-series objects and their products
     264            sage: P(L) is L
     265            True
     266
     267            sage: P(L^3)
     268            (Riemann Zeta function viewed as an L-series)^3
     269            sage: (L^3).parent()
     270            All L-series objects and their products
     271
     272        You can make the L-series attached to an object by coercing
     273        that object into the parent::
     274           
     275            sage: P(EllipticCurve('11a'))
     276            L-series of Elliptic Curve defined by y^2 + y = x^3 - x^2 - 10*x - 20 over Rational Field
     277
     278        We also allow coercing in 1, because this is very useful for
     279        formal factorizations and products::
     280
     281            sage: P(1)
     282            1
     283
     284        Other numbers do not coerce in, of course::
     285
     286            sage: P(2)
     287            Traceback (most recent call last):
     288            ...
     289            TypeError
     290        """
     291        if isinstance(x, LSeriesAbstract):
     292            return x
     293        elif isinstance(x, LSeriesProduct):
     294            return x
     295        elif x == 1:
     296            return x
     297        else:
     298            try:
     299                return LSeries(x)
     300            except NotImplementedError:
     301                raise TypeError
     302
     303    def __repr__(self):
     304        """
     305        Return string representation of this parent object.
     306
     307            sage: from sage.lfunctions.eulerprod import LSeriesParent
     308            sage: LSeriesParent.__repr__()
     309            'All L-series objects and their products'
     310        """
     311        return "All L-series objects and their products"
     312   
     313    def __cmp__(self, right):
     314        """
     315        Returns equality if right is an instance of
     316        LSeriesParentClass; otherwise compare memory locations.
     317
     318        EXAMPLES::
     319
     320            sage: from sage.lfunctions.eulerprod import LSeriesParent, LSeriesParentClass
     321            sage: cmp(LSeriesParent, LSeriesParentClass())
     322            0
     323            sage: cmp(LSeriesParent, 7) != 0
     324            True
     325            sage: cmp(7, LSeriesParent) == -cmp(LSeriesParent, 7)
     326            True
     327        """
     328        if isinstance(right, LSeriesParentClass):
     329            return 0
     330        return cmp(id(self), id(right))
     331
     332
     333LSeriesParent = LSeriesParentClass()
     334
     335class LSeriesAbstract(object):
     336    r"""
     337    L-series defined by an Euler product.
     338
     339    The parameters that define the 'shape' of the L-series are:
     340   
     341             conductor, hodge_numbers, weight, epsilon, poles, base_field
     342
     343    Let gamma(s) = prod( Gamma((s+h)/2) for h in hodge_numbers ).
     344    Denote this L-series by L(s), and let `L^*(s) = A^s \gamma(s) L(s)`, where
     345    `A = \sqrt(N)/\pi^{d/2}`, where d = len(hodge_numbers) and N = conductor.
     346    Then the functional equation is
     347
     348                  Lstar(s) = epsilon * Lstar(weight - s).
     349
     350    To actually use this class we create a derived class that in
     351    addition implements a method _local_factor(P), that takes as input
     352    a prime ideal P of K=base_field, and returns a polynomial, which
     353    is typically the reversed characteristic polynomial of Frobenius
     354    at P of Gal(Kbar/K) acting on the maximal unramified quotient of
     355    some Galois representation.  This class automatically computes the
     356    Dirichlet series coefficients `a_n` from the local factors of the
     357    `L`-function.
     358   
     359    The derived class may optionally -- and purely as an optimization
     360    -- define a method self._precompute_local_factors(bound,
     361    prec=None), which is typically called before
     362   
     363        [_local_factor(P) for P with norm(P) < bound]
     364
     365    is called in the course of various computations.
     366    """
     367    def __init__(self,
     368                 conductor,
     369                 hodge_numbers,
     370                 weight,
     371                 epsilon,
     372                 poles,
     373                 residues,
     374                 base_field,
     375                 is_selfdual=True,
     376                 prec=53):
     377        """
     378        INPUT:
     379            - ``conductor`` -- list or number (in a subset of the
     380              positive real numbers); if the conductor is a list, then
     381              each conductor is tried in order (to the precision prec
     382              below) until we find one that works.
     383            - ``hodge_numbers`` -- list of numbers (in a subring of the complex numbers)
     384            - ``weight`` -- number (in a subset of the positive real numbers)
     385            - ``epsilon`` -- number (in a subring of the complex numbers)
     386            - ``poles`` -- list of numbers (in subring of complex numbers); poles of the *completed* L-function
     387            - ``residues`` -- list of residues at each pole given in poles or string "automatic"
     388            - ``base_field`` -- QQ or a number field; local L-factors
     389              correspond to nonzero prime ideals of this field.
     390            - ``is_selfdual`` -- bool (default: True)
     391            - ``prec`` -- integer (default: 53); precision to use when trying to figure
     392              out parameters using the functional equation
     393           
     394       
     395        EXAMPLES::
     396
     397            sage: from sage.lfunctions.eulerprod import LSeriesAbstract
     398            sage: L = LSeriesAbstract(conductor=1, hodge_numbers=[0], weight=1, epsilon=1, poles=[1], residues=[-1], base_field=QQ)
     399            sage: type(L)
     400            <class 'sage.lfunctions.eulerprod.LSeriesAbstract'>
     401            sage: L._conductor
     402            1
     403            sage: L._hodge_numbers
     404            [0]
     405            sage: L._weight
     406            1
     407            sage: L._epsilon
     408            1
     409            sage: L._poles
     410            [1]
     411            sage: L._residues
     412            [-1]
     413            sage: L._base_field
     414            Rational Field
     415            sage: L
     416            Euler Product L-series with conductor 1, Hodge numbers [0], weight 1, epsilon 1, poles [1], residues [-1] over Rational Field
     417        """
     418        self._anlist = {None:[]}
     419
     420        (self._conductor, self._hodge_numbers, self._weight, self._epsilon,
     421         self._poles, self._residues, self._base_field, self._is_selfdual) = (
     422            conductor, hodge_numbers, weight, epsilon, poles, residues,
     423            base_field, is_selfdual)
     424
     425        # the following parameters allow for specifying a list of
     426        # possibilities:
     427        #
     428        #    conductor -- list of possible conductors
     429        #    hodge_numbers -- give a list of lists
     430        #    weight -- list of possible weights
     431        #    poles (residues must be the default 'automatic')
     432        #    epsilon -- list of possible epsilon's.
     433
     434        # 1. Figure out for which parameters we have multiple options
     435        # 2. Run through them until checking each until one is found
     436        #    that works.
     437        v = []
     438        if isinstance(conductor, list):
     439            v.append('_conductor')
     440        if len(hodge_numbers)>0 and isinstance(hodge_numbers[0], list):
     441            v.append('_hodge_numbers')
     442        if isinstance(weight, list):
     443            v.append('_weight')
     444        if len(poles) > 0 and isinstance(poles[0], list):
     445            v.append('_poles')
     446        if isinstance(epsilon, list):
     447            v.append('_epsilon')
     448           
     449        if len(v) > 0:
     450            found_params = False
     451            for X in cartesian_product_iterator([getattr(self, attr) for attr in v]):
     452                kwds = dict((v[i],X[i]) for i in range(len(v)))
     453                if self._is_valid_parameters(prec=prec, save=True, **kwds):
     454                    found_params = True
     455                    break
     456            if not found_params:
     457                raise RuntimeError, "no choice of values for %s works"%(', '.join(v))
     458
     459    def _is_valid_parameters(self, prec=53, save=True, **kwds):
     460        valid = False
     461        try:
     462            old = [(k, getattr(self, k)) for k in kwds.keys()]
     463            for k,v in kwds.iteritems():
     464                setattr(self, k, v)
     465            self._function.clear_cache()
     466            self._function(prec=prec)
     467            try:
     468                self._function(prec=prec)
     469                valid = True
     470            except RuntimeError:
     471                pass
     472        finally:
     473            if not save:
     474                for k, v in old:
     475                    setattr(self, k, v)
     476            return valid
     477
     478    def __cmp__(self, right):
     479        if self is right:
     480            return 0
     481        for a in ['degree', 'weight', 'conductor', 'epsilon', 'base_field']:
     482            c = cmp(getattr(self, a)(), getattr(right, a)())
     483            if c: return c
     484        c = cmp(type(self), type(right))
     485        return self._cmp(right)
     486
     487    def _cmp(self, right):
     488        raise TypeError
     489
     490    def parent(self):
     491        """
     492        Return parent of this L-series, which is the collection of all
     493        L-series and their products.
     494
     495        EXAMPLES::
     496
     497            sage: from sage.lfunctions.eulerprod import LSeries
     498            sage: L = LSeries('delta');  P = L.parent(); P
     499            All L-series objects and their products
     500            sage: L in P
     501            True
     502        """
     503        return LSeriesParent
     504
     505    def __pow__(self, n):
     506        """
     507        Return the n-th power of this L-series, where n can be any integer.
     508
     509        EXAMPLES::
     510
     511            sage: from sage.lfunctions.eulerprod import LSeries
     512            sage: L = LSeries('delta'); 
     513            sage: L3 = L^3; L3
     514            (L-function associated to Ramanujan's Delta (a weight 12 cusp form))^3
     515            sage: L3(1)
     516            0.0000524870430366548
     517            sage: L(1)^3
     518            0.0000524870430366548
     519            sage: M = L^(-3); M(1)
     520            19052.3211471761
     521            sage: L(1)^(-3)
     522            19052.3211471761
     523
     524        Higher precision::
     525
     526            sage: M = L^(-3); M(RealField(100)(1))
     527            19052.321147176093380952680193
     528            sage: L(RealField(100)(1))^(-3)
     529            19052.321147176093380952680193
     530       
     531        Special case -- 0th power -- is not allowed::
     532
     533            sage: L^0
     534            Traceback (most recent call last):
     535            ...
     536            ValueError: product must be nonempty
     537
     538        """
     539        return LSeriesProduct([(self, ZZ(n))])
     540
     541    def __mul__(self, right):
     542        """
     543        Multiply two L-series, or an L-series times a formal product of L-series.
     544
     545        EXAMPLES::
     546
     547            sage: from sage.lfunctions.eulerprod import LSeries
     548            sage: d = LSeries('delta'); z = LSeries('zeta')
     549            sage: d * z
     550            (Riemann Zeta function viewed as an L-series) * (L-function associated to Ramanujan's Delta (a weight 12 cusp form))
     551            sage: d * (d * z)
     552            (Riemann Zeta function viewed as an L-series) * (L-function associated to Ramanujan's Delta (a weight 12 cusp form))^2
     553        """
     554        if isinstance(right, LSeriesAbstract):
     555            return LSeriesProduct([(self, 1), (right, 1)])
     556        elif isinstance(right, LSeriesProduct):
     557            return right * self
     558        raise TypeError
     559   
     560    def __div__(self, right):
     561        """
     562        Divide two L-series or formal L-series products.
     563
     564        EXAMPLES::
     565
     566            sage: from sage.lfunctions.eulerprod import LSeries
     567            sage: d = LSeries('delta'); z = LSeries('zeta')
     568            sage: d / z
     569            (Riemann Zeta function viewed as an L-series)^-1 * (L-function associated to Ramanujan's Delta (a weight 12 cusp form))
     570            sage: d / (z^3)
     571            (Riemann Zeta function viewed as an L-series)^-3 * (L-function associated to Ramanujan's Delta (a weight 12 cusp form))
     572        """
     573        if isinstance(right, LSeriesAbstract):
     574            return LSeriesProduct([(self, 1), (right, -1)])
     575        elif isinstance(right, LSeriesProduct):
     576            return LSeriesProduct(Factorization([(self, 1)]) / right._factorization)
     577        raise TypeError
     578
     579    def conductor(self):
     580        """
     581        Return the conductor of this L-series.
     582
     583        EXAMPLES::
     584
     585            sage: from sage.lfunctions.eulerprod import LSeries
     586            sage: LSeries('zeta').conductor()
     587            1
     588            sage: LSeries('delta').conductor()
     589            1
     590            sage: LSeries(EllipticCurve('11a')).conductor()
     591            11
     592            sage: LSeries(Newforms(33)[0]).conductor()
     593            33
     594            sage: LSeries(DirichletGroup(37).0).conductor()
     595            37
     596            sage: LSeries(kronecker_character(7)).conductor()
     597            28
     598            sage: kronecker_character(7).conductor()
     599            28
     600            sage: L = LSeries(EllipticCurve('11a3').base_extend(QQ[sqrt(2)]), prec=5)
     601            sage: L.conductor().factor()
     602            2^6 * 11^2
     603        """
     604        return self._conductor
     605
     606    def hodge_numbers(self):
     607        """
     608        Return the Hodge numbers of this L-series.  These define the local Gamma factors.
     609
     610        EXAMPLES::
     611
     612            sage: from sage.lfunctions.eulerprod import LSeries
     613            sage: LSeries('zeta').hodge_numbers()
     614            [0]
     615            sage: LSeries('delta').hodge_numbers()
     616            [0, 1]
     617            sage: LSeries(EllipticCurve('11a')).hodge_numbers()
     618            [0, 1]
     619            sage: LSeries(Newforms(43,names='a')[1]).hodge_numbers()
     620            [0, 1]
     621            sage: LSeries(DirichletGroup(37).0).hodge_numbers()
     622            [1]
     623            sage: LSeries(EllipticCurve(QQ[sqrt(-1)],[1,2]), prec=5).hodge_numbers() # long time
     624            [0, 0, 1, 1]
     625        """
     626        return list(self._hodge_numbers)
     627
     628    def weight(self):
     629        """
     630        Return the weight of this L-series.
     631       
     632        EXAMPLES::
     633       
     634            sage: from sage.lfunctions.eulerprod import LSeries
     635            sage: LSeries('zeta').weight()
     636            1
     637            sage: LSeries('delta').weight()
     638            12
     639            sage: LSeries(EllipticCurve('389a')).weight()
     640            2
     641            sage: LSeries(Newforms(43,names='a')[1]).weight()
     642            2
     643            sage: LSeries(Newforms(6,4)[0]).weight()
     644            4
     645            sage: LSeries(DirichletGroup(37).0).weight()
     646            1
     647            sage: L = LSeries(EllipticCurve('11a3').base_extend(QQ[sqrt(2)]),prec=5); L.weight()
     648            2
     649        """
     650        return self._weight
     651
     652    def poles(self):
     653        """
     654        Poles of the *completed* L-function with the extra Gamma
     655        factors included.
     656
     657        WARNING: These are not just the poles of self.
     658
     659        OUTPUT:
     660             - list of numbers
     661       
     662        EXAMPLES::
     663
     664            sage: from sage.lfunctions.eulerprod import LSeries
     665            sage: LSeries('zeta').poles()
     666            [0, 1]
     667            sage: LSeries('delta').poles()
     668            []
     669        """
     670        return list(self._poles)
     671   
     672    def residues(self, prec=None):
     673        """
     674        Residues of the *completed* L-function at each pole.
     675       
     676        EXAMPLES::
     677
     678            sage: from sage.lfunctions.eulerprod import LSeries
     679            sage: LSeries('zeta').residues()
     680            [9, 8]
     681            sage: LSeries('delta').residues()
     682            []
     683            sage: v = LSeries('zeta').residues()
     684            sage: v.append(10)
     685            sage: LSeries('zeta').residues()
     686            [9, 8]
     687
     688        The residues of the Dedekind Zeta function of a field are dynamically computed::
     689
     690            sage: K.<a> = NumberField(x^2 + 1)
     691            sage: L = LSeries(K); L
     692            Dedekind Zeta function of Number Field in a with defining polynomial x^2 + 1
     693
     694        If you just call residues you get back that they are automatically computed::
     695       
     696            sage: L.residues()
     697            'automatic'
     698
     699        But if you call with a specific precision, they are computed using that precision::
     700       
     701            sage: L.residues(prec=53)
     702            [-0.886226925452758]
     703            sage: L.residues(prec=200)
     704            [-0.88622692545275801364908374167057259139877472806119356410690]
     705        """
     706        if self._residues == 'automatic':
     707            if prec is None:
     708                return self._residues
     709            else:
     710                C = ComplexField(prec)
     711                return [C(a) for a in self._function(prec=prec).gp()('Lresidues')]
     712        else:
     713            return list(self._residues)
     714
     715    def base_field(self):
     716        """
     717        EXAMPLES::
     718
     719            sage: from sage.lfunctions.eulerprod import LSeries
     720            sage: LSeries('zeta').base_field()
     721            Rational Field
     722            sage: L = LSeries(EllipticCurve('11a3').base_extend(QQ[sqrt(2)]), prec=5); L.base_field()
     723            Number Field in sqrt2 with defining polynomial x^2 - 2
     724        """
     725        return self._base_field
     726
     727    def epsilon(self, prec=None):
     728        """
     729        EXAMPLES::
     730
     731            sage: from sage.lfunctions.eulerprod import LSeries
     732            sage: LSeries('zeta').epsilon()
     733            1
     734            sage: LSeries('delta').epsilon()
     735            1
     736            sage: LSeries(EllipticCurve('389a')).epsilon()
     737            1
     738            sage: LSeries(EllipticCurve('37a')).epsilon()
     739            -1
     740            sage: LSeries(Newforms(389,names='a')[1]).epsilon()
     741            -1
     742            sage: LSeries(Newforms(6,4)[0]).epsilon()
     743            1
     744            sage: LSeries(DirichletGroup(7).0).epsilon()
     745            1/7*I*((((((e^(2/21*I*pi) + 1)*e^(2/21*I*pi) + 1)*e^(1/21*I*pi) - 1)*e^(1/21*I*pi) - 1)*e^(2/21*I*pi) - 1)*e^(1/21*I*pi) - 1)*sqrt(7)*e^(1/21*I*pi)
     746            sage: LSeries(DirichletGroup(7).0).epsilon(prec=53)
     747            0.386513572759156 + 0.922283718859307*I
     748
     749        In this example, the epsilon factor is computed when the curve
     750        is created. The prec parameter determines the floating point precision
     751        used in computing the epsilon factor::
     752       
     753            sage: L = LSeries(EllipticCurve('11a3').base_extend(QQ[sqrt(2)]), prec=5); L.epsilon()
     754            -1
     755
     756        Here is extra confirmation that the rank is really odd over the quadratic field::
     757       
     758            sage: EllipticCurve('11a').quadratic_twist(2).rank()
     759            1
     760
     761        We can compute with the L-series too::
     762       
     763            sage: L(RealField(5)(2))
     764            0.53
     765
     766        For L-functions of newforms with nontrivial character, the
     767        epsilon factor is harder to find (we don't have a good
     768        algorithm implemented to find it) and might not even be 1 or
     769        -1, so it is set to 'solve'.  In this case, the functional
     770        equation is used to determine the solution.::
     771
     772            sage: f = Newforms(kronecker_character_upside_down(7),3)[0]; f
     773            q - 3*q^2 + 5*q^4 + O(q^6)
     774            sage: L = LSeries(f)
     775            sage: L.epsilon()
     776            'solve'
     777            sage: L(3/2)
     778            0.332981771482934
     779            sage: L.epsilon()
     780            1
     781
     782        Here is an example with nontrivial character::
     783
     784            sage: f = Newforms(DirichletGroup(7).0, 5, names='a')[0]; f
     785            q + a0*q^2 + ((zeta6 - 2)*a0 - zeta6 - 1)*q^3 + (-4*zeta6*a0 + 2*zeta6 - 2)*q^4 + ((4*zeta6 - 2)*a0 + 9*zeta6 - 18)*q^5 + O(q^6)
     786            sage: L = LSeries(f)
     787
     788        First trying to evaluate with the default (53 bits) of
     789        precision fails, since for some reason (that I do not
     790        understand) the program is unable to find a valid epsilon
     791        factor::
     792
     793            sage: L(0)
     794            Traceback (most recent call last):
     795            ...
     796            RuntimeError: unable to determine epsilon from functional equation working to precision 53, since we get epsilon=0.806362085925390 - 0.00491051026156292*I, which is not sufficiently close to 1
     797
     798        However, when we evaluate to 100 bits of precision it works::
     799       
     800            sage: L(RealField(100)(0))
     801            0
     802
     803        The epsilon factor is *not* known to infinite precision::
     804           
     805            sage: L.epsilon()
     806            'solve'
     807
     808        But it is now known to 100 bits of precision, and here it is::
     809       
     810            sage: L.epsilon(100)
     811            0.42563106101692403875896879406 - 0.90489678963824790765479396740*I
     812
     813        When we try to compute to higher precision, again Sage solves for the epsilon factor
     814        numerically::
     815           
     816            sage: L(RealField(150)(1))
     817            0.26128389551787271923496480408992971337929665 - 0.29870133769674001421149135036267324347896657*I
     818
     819        And now it is known to 150 bits of precision.  Notice that
     820        this is consistent with the value found above, and has
     821        absolute value (very close to) 1.
     822       
     823            sage: L.epsilon(150)
     824            0.42563106101692403875896879406038776338921622 - 0.90489678963824790765479396740501409301704122*I
     825            sage: abs(L.epsilon(150))
     826            1.0000000000000000000000000000000000000000000
     827        """
     828        if self._epsilon == 'solve' or (
     829                  hasattr(self._epsilon, 'prec') and (prec is None or self._epsilon.prec() < prec)):
     830            return 'solve'
     831        if prec is not None:
     832            C = ComplexField(prec)
     833            if isinstance(self._epsilon, list):
     834                return [C(x) for x in self._epsilon]
     835            return C(self._epsilon)
     836        return self._epsilon
     837
     838    def is_selfdual(self):
     839        """
     840        Return True if this L-series is self dual; otherwise, return False.
     841
     842        EXAMPLES::
     843
     844        Many L-series are self dual::
     845       
     846            sage: from sage.lfunctions.eulerprod import LSeries
     847            sage: LSeries('zeta').is_selfdual()
     848            True
     849            sage: LSeries('delta').is_selfdual()
     850            True
     851            sage: LSeries(Newforms(6,4)[0]).is_selfdual()
     852            True
     853
     854        Nonquadratic characters have non-self dual L-series::
     855       
     856            sage: LSeries(DirichletGroup(7).0).is_selfdual()
     857            False
     858            sage: LSeries(kronecker_character(7)).is_selfdual()
     859            True
     860
     861        Newforms with non-quadratic characters also have non-self dual L-seris::
     862           
     863            sage: L = LSeries(Newforms(DirichletGroup(7).0, 5, names='a')[0]); L.is_selfdual()
     864            False
     865        """
     866        return self._is_selfdual
     867
     868    @cached_method
     869    def degree(self):
     870        """
     871        Return the degree of this L-function, which is by definition the number of Gamma
     872        factors (e.g., the number of Hodge numbers) divided by the degree of the base field.
     873
     874        EXAMPLES::
     875
     876            sage: from sage.lfunctions.eulerprod import LSeries
     877            sage: LSeries('zeta').degree()
     878            1
     879            sage: LSeries(DirichletGroup(5).0).degree()
     880            1
     881            sage: LSeries(EllipticCurve('11a')).degree()
     882            2
     883
     884        The L-series attached to this modular symbols space of dimension 2 is a product
     885        of 2 degree 2 L-series, hence has degree 4::
     886
     887            sage: M = ModularSymbols(43,2,sign=1).cuspidal_subspace()[1]; M.dimension()
     888            2
     889            sage: L = LSeries(M); L
     890            L-series attached to Modular Symbols subspace of dimension 2 of Modular Symbols space of dimension 4 for Gamma_0(43) of weight 2 with sign 1 over Rational Field
     891            sage: L.factor()
     892            (L-series of a degree 2 newform of level 43 and weight 2) * (L-series of a degree 2 newform of level 43 and weight 2)
     893            sage: L.degree()
     894            4
     895           
     896            sage: x = var('x'); K.<a> = NumberField(x^2-x-1); LSeries(EllipticCurve([0,-a,a,0,0])).degree()
     897            2
     898        """
     899        n = len(self.hodge_numbers())
     900        d = self.base_field().degree()
     901        assert n % d == 0, "degree of base field must divide the number of Hodge numbers"
     902        return n//d
     903
     904    def twist(self, chi, conductor=None, epsilon=None, prec=53):
     905        r"""
     906        Return the quadratic twist of this L-series by the character chi, which
     907        must be a character of self.base_field().   Thus chi should take as input
     908        prime ideals (or primes) of the ring of integers of the base field, and
     909        output something that can be coerced to the complex numbers.
     910
     911        INPUT:
     912            - `\chi` -- 1-dimensional character
     913
     914        EXAMPLES::
     915
     916            sage: from sage.lfunctions.eulerprod import LSeries
     917            sage: E = EllipticCurve('11a')
     918            sage: L = LSeries(E); L
     919            L-series of Elliptic Curve defined by y^2 + y = x^3 - x^2 - 10*x - 20 over Rational Field
     920            sage: L3 = L.twist(DirichletGroup(3).0); L3
     921            Twist of L-series of Elliptic Curve defined by y^2 + y = x^3 - x^2 - 10*x - 20 over Rational Field by Dirichlet character modulo 3 of conductor 3 mapping 2 |--> -1
     922            sage: L3._chi
     923            Dirichlet character modulo 3 of conductor 3 mapping 2 |--> -1
     924            sage: L3._L
     925            L-series of Elliptic Curve defined by y^2 + y = x^3 - x^2 - 10*x - 20 over Rational Field
     926            sage: L3(1)
     927            1.68449633297548
     928            sage: F = E.quadratic_twist(-3)
     929            sage: L3.conductor()
     930            99
     931            sage: F.conductor()
     932            99
     933            sage: F.lseries()(1)
     934            1.68449633297548
     935            sage: L3.anlist(20)
     936            [0, 1, 2, 0, 2, -1, 0, -2, 0, 0, -2, -1, 0, 4, -4, 0, -4, 2, 0, 0, -2]
     937            sage: F.anlist(20)
     938            [0, 1, 2, 0, 2, -1, 0, -2, 0, 0, -2, -1, 0, 4, -4, 0, -4, 2, 0, 0, -2]
     939            sage: L3.anlist(1000) == F.anlist(1000)
     940            True
     941            sage: L3.local_factor(11)
     942            T + 1
     943
     944        A higher degree twist::
     945
     946            sage: L = LSeries(EllipticCurve('11a'))
     947            sage: L5 = L.twist(DirichletGroup(5).0); L5
     948            Twist of L-series of Elliptic Curve defined by y^2 + y = x^3 - x^2 - 10*x - 20 over Rational Field by Dirichlet character modulo 5 of conductor 5 mapping 2 |--> zeta4
     949            sage: L5(1)
     950            1.28009593569230 - 0.681843202124309*I
     951            sage: L5.epsilon()
     952            'solve'
     953            sage: L5.epsilon(53)
     954            0.989989082587826 + 0.141143956147310*I
     955            sage: L5.conductor()
     956            275
     957            sage: L5.taylor_series(center=1, degree=3)
     958            1.28009593569230 - 0.681843202124309*I + (-0.536450338806282 + 0.166075270978779*I)*z + (0.123743053129226 + 0.320802890011298*I)*z^2 + O(z^3)
     959
     960
     961        WARNING!! Twisting is not implemented in full generality when
     962        the conductors are not coprime.  One case where we run into
     963        trouble is when twisting lowers the level of a newform.  Below
     964        we take the form of level 11 and weight 2, twist it by the
     965        character chi of conductor 3 to get a form of level 99.  Then
     966        we take the L-series of the level 99 form, and twist that by
     967        chi, which should be the L-series attached to the form of
     968        level 11.  Unfortunately, our code for working out the local
     969        L-factors doesn't succeed in this case, hence the local factor
     970        is wrong, so the functional equation is not satisfied.
     971
     972            sage: f = Newform('11a'); f
     973            q - 2*q^2 - q^3 + 2*q^4 + q^5 + O(q^6)
     974            sage: L = LSeries(f)
     975            sage: chi = DirichletGroup(3).0
     976            sage: Lc = L.twist(chi); Lc
     977            Twist of L-series of a degree 1 newform of level 11 and weight 2 by Dirichlet character modulo 3 of conductor 3 mapping 2 |--> -1
     978            sage: Lc.anlist(20)
     979            [0, 1, 2, 0, 2, -1, 0, -2, 0, 0, -2, -1, 0, 4, -4, 0, -4, 2, 0, 0, -2]
     980            sage: g = Newform('99d'); g
     981            q + 2*q^2 + 2*q^4 - q^5 + O(q^6)
     982            sage: list(g.qexp(20))
     983            [0, 1, 2, 0, 2, -1, 0, -2, 0, 0, -2, -1, 0, 4, -4, 0, -4, 2]
     984            sage: Lt = Lc.twist(chi, conductor=11)
     985            Traceback (most recent call last):
     986            ...
     987            RuntimeError: no choice of values for _epsilon works
     988            sage: Lt = Lc.twist(chi,conductor=11, epsilon=1)
     989            sage: Lt(1)
     990            Traceback (most recent call last):
     991            ...
     992            RuntimeError: invalid L-series parameters: functional equation not satisfied
     993
     994        This is because the local factor is wrong::
     995       
     996            sage: Lt.local_factor(3)
     997            1
     998            sage: L.local_factor(3)
     999            3*T^2 + T + 1
     1000           
     1001       
     1002        """
     1003        return LSeriesTwist(self, chi=chi, conductor=conductor, epsilon=epsilon, prec=prec)
     1004
     1005    @cached_method
     1006    def local_factor(self, P, prec=None):
     1007        """
     1008        Return the local factor of the L-function at the prime P of
     1009        self._base_field.  The result is cached.
     1010       
     1011        INPUT:
     1012            - a prime P of the ring of integers of the base_field
     1013            - prec -- None or positive integer (bits of precision)
     1014           
     1015        OUTPUT:
     1016            - a polynomial, e.g., something like "1-a*T+p*T^2".
     1017
     1018        EXAMPLES:
     1019
     1020        You must overload this in the derived class::
     1021
     1022            sage: from sage.lfunctions.eulerprod import LSeriesAbstract
     1023            sage: L = LSeriesAbstract(conductor=1, hodge_numbers=[0], weight=1, epsilon=1, poles=[1], residues=[-1], base_field=QQ)
     1024            sage: L.local_factor(2)
     1025            Traceback (most recent call last):
     1026            ...
     1027            NotImplementedError: must be implemented in the derived class
     1028        """
     1029        return self._local_factor(P, prec=prec)
     1030
     1031    def _local_factor(self, P, prec=None):
     1032        """
     1033        Compute local factor at prime P.  This must be overwritten in the derived class.
     1034
     1035        EXAMPLES::
     1036
     1037            sage: from sage.lfunctions.eulerprod import LSeriesAbstract
     1038            sage: L = LSeriesAbstract(conductor=1, hodge_numbers=[0], weight=1, epsilon=1, poles=[1], residues=[-1], base_field=QQ)
     1039            sage: L.local_factor(2)
     1040            Traceback (most recent call last):
     1041            ...
     1042            NotImplementedError: must be implemented in the derived class
     1043        """
     1044        raise NotImplementedError, "must be implemented in the derived class"
     1045
     1046    def __repr__(self):
     1047        """
     1048        EXAMPLES::
     1049       
     1050            sage: from sage.lfunctions.eulerprod import LSeriesAbstract
     1051            sage: L = LSeriesAbstract(conductor=1, hodge_numbers=[0], weight=1, epsilon=1, poles=[1], residues=[-1], base_field=QQ)
     1052            sage: L.__repr__()
     1053            'Euler Product L-series with conductor 1, Hodge numbers [0], weight 1, epsilon 1, poles [1], residues [-1] over Rational Field'
     1054        """
     1055        return "Euler Product L-series with conductor %s, Hodge numbers %s, weight %s, epsilon %s, poles %s, residues %s over %s"%(
     1056            self._conductor, self._hodge_numbers, self._weight, self.epsilon(), self._poles, self._residues, self._base_field)
     1057
     1058    def _precompute_local_factors(self, bound, prec):
     1059        """
     1060        Derived classes may use this as a 'hint' that _local_factors
     1061        will soon get called for primes of norm less than the bound.
     1062
     1063        In the base class, this is a no-op, and it is not necessary to
     1064        overload this class.
     1065
     1066        INPUT:
     1067            - ``bound`` -- integer
     1068            - ``prec`` -- integer
     1069
     1070        EXAMPLES::
     1071
     1072            sage: from sage.lfunctions.eulerprod import LSeriesAbstract
     1073            sage: L = LSeriesAbstract(conductor=1, hodge_numbers=[0], weight=1, epsilon=1, poles=[1], residues=[-1], base_field=QQ)
     1074            sage: L._precompute_local_factors(100, 53)
     1075        """
     1076        pass
     1077
     1078    def _primes_above(self, p):
     1079        """
     1080        Return the primes of the ring of integers of the base field above the integer p.
     1081
     1082        INPUT:
     1083            - p -- prime integer (no type checking necessarily done)
     1084
     1085        EXAMPLES::
     1086
     1087            sage: from sage.lfunctions.eulerprod import LSeriesAbstract
     1088            sage: L = LSeriesAbstract(conductor=1, hodge_numbers=[0], weight=1, epsilon=1, poles=[1], residues=[-1], base_field=QQ)
     1089            sage: L._primes_above(3)
     1090            [3]
     1091            sage: L = LSeriesAbstract(conductor=1, hodge_numbers=[0], weight=1, epsilon=1, poles=[1], residues=[-1], base_field=QQ[sqrt(-1)])
     1092            sage: L._primes_above(5)
     1093            [Fractional ideal (I + 2), Fractional ideal (-I + 2)]
     1094            sage: L._primes_above(3)
     1095            [Fractional ideal (3)]
     1096        """
     1097        K = self._base_field
     1098        if is_RationalField(K):
     1099            return [p]
     1100        else:
     1101            return K.primes_above(p)
     1102
     1103    def zeros(self, n):
     1104        """
     1105        Return the imaginary parts of the first n nontrivial zeros of
     1106        this L-function on the critical line in the upper half plane,
     1107        as 32-bit real numbers.
     1108
     1109        INPUT:
     1110             - n -- nonnegative integer
     1111
     1112        EXAMPLES::
     1113       
     1114        """
     1115        return self._lcalc().zeros(n)
     1116
     1117    def _lcalc(self):
     1118        """
     1119        Return Rubinstein Lcalc object attached to this L-series. This
     1120        is useful both for evaluating the L-series, especially when
     1121        the imaginary part is large, and for computing the zeros in
     1122        the critical strip.
     1123
     1124        EXAMPLES::
     1125        """
     1126        # this will require using Rubinstein's L-calc
     1127        raise NotImplementedError
     1128
     1129    def anlist(self, bound, prec=None):
     1130        """
     1131        Return list `v` of Dirichlet series coefficients `a_n`for `n` up to and including bound,
     1132        where `v[n] = a_n`.  If at least one of the coefficients is ambiguous, e.g., the
     1133        local_factor method returns a list of possibilities, then this function instead
     1134        returns a generator over all possible `a_n` lists.
     1135
     1136        In particular, we include v[0]=0 as a convenient place holder
     1137        to avoid having `v[n-1] = a_n`.
     1138
     1139        INPUT:
     1140            - ``bound`` -- nonnegative integer
     1141       
     1142        EXAMPLES::
     1143
     1144            sage: from sage.lfunctions.eulerprod import LSeries; L = LSeries('zeta')
     1145            sage: L.anlist(30)
     1146            [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
     1147            sage: from sage.lfunctions.eulerprod import LSeries; L = LSeries(EllipticCurve('11a'))
     1148            sage: L.anlist(30)
     1149            [0, 1, -2, -1, 2, 1, 2, -2, 0, -2, -2, 1, -2, 4, 4, -1, -4, -2, 4, 0, 2, 2, -2, -1, 0, -4, -8, 5, -4, 0, 2]
     1150            sage: K.<a> = NumberField(x^2-x-1); L = LSeries(EllipticCurve([0,-a,a,0,0]))
     1151            sage: L.anlist(30)
     1152            [0, 1, 0, 0, -2, -1, 0, 0, 0, -4, 0, 3, 0, 0, 0, 0, 0, 0, 0, 5, 2, 0, 0, 0, 0, -4, 0, 0, 0, 11, 0]
     1153        """
     1154        # First check if we know anlist to infinite bit precision up to given bound:
     1155        if len(self._anlist[None]) > bound:
     1156            if prec is None:
     1157                # request it to infinite precision
     1158                return self._anlist[None][:bound+1]
     1159            else:
     1160                # finite precision request
     1161                C = ComplexField(prec)
     1162                return [C(a) for a in self._anlist[None]]
     1163
     1164        if prec is not None:
     1165            # check if numerically computed already to at least this precision
     1166            t = [z for z in self._anlist.iteritems() if z[0] >= prec and len(z[1]) > bound]
     1167            if len(t) > 0:
     1168                C = ComplexField(prec)
     1169                return [C(a) for a in t[0][1][:bound+1]]
     1170
     1171        self._precompute_local_factors(bound+1, prec=prec)
     1172        compute_anlist_multiple = False
     1173        LF = []
     1174        for p in prime_range(bound+1):
     1175            lf = []
     1176            for P in self._primes_above(p):
     1177                if norm(P) <= bound:
     1178                    F = self._local_factor(P, prec)
     1179                    if isinstance(F, list):
     1180                        compute_anlist_multiple = True
     1181                    lf.append(F)
     1182            LF.append((p, lf))
     1183                       
     1184        coefficients = self._compute_anlist(LF, bound, prec)
     1185        if not compute_anlist_multiple:
     1186            coefficients = list(coefficients)[0]
     1187            # save answer in cache
     1188            self._anlist[prec] = coefficients
     1189        return coefficients
     1190
     1191    def _compute_anlist(self, LF, bound, prec):
     1192        """
     1193        Iterator over possible anlists, given LF, bound, and prec.
     1194
     1195        INPUT:
     1196            - ``LF`` -- list of pairs (p, [local factors (or lists of them) at primes over p])
     1197            - ``bound`` -- positive integer
     1198            - ``prec`` -- positive integer (bits of precision)
     1199        """
     1200        K = self._base_field
     1201        coefficients = [0,1] + [0]*(bound-1)
     1202       
     1203        for i, (p, v) in enumerate(LF):
     1204            if len(v) > 0:
     1205                # Check for the possibility of several different
     1206                # choices of Euler factor at a given prime.   If this happens,
     1207                # we switch gears.
     1208                some_list = False
     1209                for j, lf in enumerate(v):
     1210                    if isinstance(lf, list):
     1211                        some_list = True
     1212                        for f in list(lf):
     1213                            LF0 = copy.deepcopy(LF)
     1214                            LF0[i][1][j] = f
     1215                            for z in self._compute_anlist(LF0, bound, prec):
     1216                                yield z
     1217                if some_list:
     1218                    return
     1219                # Not several factors -- just compute the a_{p^r} up to the required bound:
     1220                f = prod(v)
     1221                accuracy_p = int(math.floor(math.log(bound)/math.log(p))) + 1
     1222                T = f.parent().gen()
     1223                series_p = (f + O(T**accuracy_p))**(-1)
     1224                for j in range(1, accuracy_p):
     1225                    coefficients[p**j] = series_p[j]
     1226                   
     1227        # fill in non-prime power coefficients
     1228        extend_multiplicatively_generic(coefficients)
     1229        yield list(coefficients)
     1230
     1231    def _symbolic_(self, R, bound=10, prec=None):
     1232        """
     1233        Convert self into the symbolic ring as a truncated Dirichleter series, including
     1234        terms up to `n^s` where n=bound.
     1235
     1236        EXAMPLES::
     1237       
     1238            sage: from sage.lfunctions.eulerprod import LSeries
     1239            sage: L = LSeries(EllipticCurve('37a'))
     1240            sage: SR(L)
     1241            -2/2^s - 3/3^s + 2/4^s - 2/5^s + 6/6^s - 1/7^s + 6/9^s + 4/10^s + 1           
     1242            sage: L._symbolic_(SR, 20)
     1243            -2/2^s - 3/3^s + 2/4^s - 2/5^s + 6/6^s - 1/7^s + 6/9^s + 4/10^s - 5/11^s - 6/12^s - 2/13^s + 2/14^s + 6/15^s - 4/16^s - 12/18^s - 4/20^s + 1
     1244        """
     1245        s = R.var('s')
     1246        a = self.anlist(bound, prec)
     1247        return sum(a[n]/n**s for n in range(1,bound+1))
     1248
     1249    def __call__(self, s):
     1250        """
     1251        Return value of this L-function at s.  If s is a real or
     1252        complex number to prec bits of precision, then the result is
     1253        also computed to prec bits of precision.  If s has infinite or
     1254        unknown precision, then the L-value is computed to 53 bits of
     1255        precision.
     1256
     1257        EXAMPLES::
     1258
     1259            sage: from sage.lfunctions.eulerprod import LSeries
     1260            sage: L = LSeries(EllipticCurve('37a'))
     1261            sage: L(1)
     1262            0
     1263            sage: L(2)
     1264            0.381575408260711
     1265            sage: z = L(RealField(100)(2)); z
     1266            0.38157540826071121129371040958
     1267            sage: z.prec()
     1268            100
     1269            sage: L(RealField(150)(2))
     1270            0.38157540826071121129371040958008663667709753
     1271
     1272        WARNING: There will be precision loss (with a warning) if the imaginary
     1273        part of the input is large::
     1274
     1275            sage: L = LSeries('zeta')
     1276            sage: L(1/2 + 40*I)
     1277            verbose -1 (...: dokchitser.py, __call__) Warning: Loss of 14 decimal digits due to cancellation
     1278            0.793046013671137 - 1.04127377821427*I
     1279            sage: L(ComplexField(200)(1/2 + 40*I))
     1280            verbose -1 (...: dokchitser.py, __call__) Warning: Loss of 14 decimal digits due to cancellation
     1281            0.79304495256192867196489258889793696080572220439302833315881 - 1.0412746146510650200518905953910554313275550685861559488384*I
     1282
     1283        An example with a pole::
     1284
     1285            sage: L = LSeries('zeta')
     1286            sage: L(2)  # good
     1287            1.64493406684823
     1288            sage: L(1)  # a pole!
     1289            Traceback (most recent call last):
     1290            ...
     1291            ZeroDivisionError: pole at 1
     1292        """
     1293        if s in self._poles:
     1294            raise ZeroDivisionError, "pole at %s"%s
     1295        return self._function(prec(s))(s)
     1296
     1297    def derivative(self, k=1):
     1298        """
     1299        Return the k-th derivative of self.
     1300
     1301        INPUT:
     1302            - k -- (default: 1) nonnegative integer
     1303
     1304        EXAMPLES::
     1305       
     1306            sage: from sage.lfunctions.eulerprod import LSeries
     1307            sage: L = LSeries('zeta')
     1308            sage: L.derivative()
     1309            First derivative of Riemann Zeta function viewed as an L-series
     1310
     1311        We numerically approximate the derivative at two points and compare
     1312        with evaluating the derivative object::
     1313       
     1314            sage: eps=1e-10;  (L(2+eps) - L(2))/eps
     1315            -0.937547817159157
     1316            sage: Lp = L.derivative(); Lp(2)
     1317            -0.937548254315844
     1318            sage: eps=1e-10;  (L(2+I+eps) - L(2+I))/eps
     1319            0.0624900131640516 + 0.489033813444451*I
     1320            sage: Lp(2+I)
     1321            0.0624900021906470 + 0.489033591679899*I
     1322       
     1323        Higher derivatives::
     1324       
     1325            sage: L.derivative(2)
     1326            Second derivative of Riemann Zeta function viewed as an L-series
     1327            sage: L.derivative(2)(2)
     1328            1.98928023429890
     1329            sage: L.derivative(3)
     1330            Third derivative of Riemann Zeta function viewed as an L-series
     1331            sage: L.derivative(4)
     1332            4-th derivative of Riemann Zeta function viewed as an L-series
     1333            sage: L.derivative(5)
     1334            5-th derivative of Riemann Zeta function viewed as an L-series
     1335
     1336        Derivative of derivative::
     1337       
     1338            sage: L.derivative().derivative()
     1339            Second derivative of Riemann Zeta function viewed as an L-series
     1340
     1341        Using the derivative function in Sage works::
     1342       
     1343            sage: derivative(L)
     1344            First derivative of Riemann Zeta function viewed as an L-series
     1345            sage: derivative(L,2)
     1346            Second derivative of Riemann Zeta function viewed as an L-series
     1347        """
     1348        if k == 0:
     1349            return self
     1350        return LSeriesDerivative(self, k)
     1351
     1352    def taylor_series(self, center=None, degree=6, variable='z', prec=53):
     1353        """
     1354        Return the Taylor series expansion of self about the given
     1355        center point to the given degree in the specified variable
     1356        numerically computed to the precision prec in bits.  If the
     1357        center is not specified it defaults to weight / 2.
     1358
     1359        INPUT:
     1360            - ``center`` -- None or number that coerces to the complex numbers
     1361            - ``degree`` -- integer
     1362            - ``variable`` -- string or symbolic variable
     1363            - ``prec`` -- positive integer (floating point bits of precision)
     1364
     1365        EXAMPLES:::
     1366
     1367            sage: from sage.lfunctions.eulerprod import LSeries; L = LSeries('zeta')
     1368            sage: L.taylor_series()
     1369            -1.46035450880959 - 3.92264613920915*z - 8.00417850696433*z^2 - 16.0005515408865*z^3 - 31.9998883216853*z^4 - 64.0000050055172*z^5 + O(z^6)
     1370            sage: RealField(53)(zeta(1/2))
     1371            -1.46035450880959
     1372            sage: L.taylor_series(center=2, degree=4, variable='t', prec=30)
     1373            1.6449341 - 0.93754825*t + 0.99464012*t^2 - 1.0000243*t^3 + O(t^4)
     1374            sage: RealField(30)(zeta(2))
     1375            1.6449341
     1376             
     1377        """
     1378        if center is None:
     1379            center = ComplexField(prec)(self._weight)/2
     1380        return self._function(prec).taylor_series(center, degree, variable)
     1381
     1382    def analytic_rank(self, tiny=1e-8, prec=53):
     1383        center = ComplexField(prec)(self._weight) / 2
     1384        degree = 4
     1385        while True:
     1386            f = self.taylor_series(center, degree, prec=prec)
     1387            i = 0
     1388            while i < degree:
     1389                if abs(f[i]) > tiny:
     1390                    return i
     1391                i += 1
     1392            degree += 2
     1393
     1394    @cached_method
     1395    def _function(self, prec=53, T=1.2):
     1396        """
     1397        Return Dokchitser object that allows for computation of this
     1398        L-series computed to enough terms so that the functional
     1399        equation checks out with the given value of T and precision.
     1400
     1401        This is used behind the scenes for evaluation and computation
     1402        of Taylor series.
     1403        """
     1404        eps = self.epsilon(prec)
     1405        return self._dokchitser(prec, eps, T=T)
     1406
     1407    def _dokchitser_unitialized(self, prec, epsilon):
     1408        # Create the Dokchitser object
     1409        if epsilon == 'solve':
     1410            eps = 'X'
     1411        else:
     1412            eps = epsilon
     1413        return Dokchitser(conductor = self.conductor(), gammaV = self.hodge_numbers(), weight = self.weight(),
     1414                       eps = eps, poles = self.poles(), residues = self.residues(),
     1415                       prec = prec)
     1416
     1417    def number_of_coefficients(self, prec=53, T=1.2):
     1418        """
     1419        Return the number of Dirichlet series coefficients that will
     1420        be needed in order to evaluate this L-series (near the real
     1421        line) to prec bits of precision and have the functional
     1422        equation test pass with the given value of T.
     1423
     1424        INPUT:
     1425            - prec -- integer
     1426       
     1427        EXAMPLES::
     1428
     1429            sage: from sage.lfunctions.eulerprod import LSeries       
     1430            sage: L = LSeries(DirichletGroup(5).0)
     1431            sage: L.number_of_coefficients(20)
     1432            8
     1433            sage: L.number_of_coefficients()
     1434            11
     1435            sage: L.number_of_coefficients(1000)
     1436            43
     1437        """
     1438        return self._dokchitser_unitialized(prec, self.epsilon(prec)).num_coeffs(T)
     1439
     1440    def _dokchitser(self, prec, epsilon, T=1.2):
     1441        L = self._dokchitser_unitialized(prec, epsilon)
     1442       
     1443        # Find out how many coefficients of the Dirichlet series are needed
     1444        # to compute to the requested precision.
     1445        n = L.num_coeffs(T=T)
     1446        #if n >= 500:   # TODO: for debugging only -- remove later
     1447        #    print "num coeffs =", n
     1448
     1449        # Compute the Dirichlet series coefficients
     1450        X = self.anlist(n, prec)
     1451        if isinstance(X, types.GeneratorType):
     1452            # Several possible coefficients -- we try them until finding on that works.
     1453            coeff_lists = X
     1454        else:
     1455            # Only one to try
     1456            coeff_lists = [X]
     1457
     1458        tiny0 = tiny(prec)
     1459        for coeffs in coeff_lists:
     1460            # Define a string that when evaluated in PARI defines a function
     1461            # a(k), which returns the Dirichlet coefficient a_k.
     1462            s = 'v=[%s]; a(k)=v[k];'%','.join([str(z) if isinstance(z, (int,long,Integer)) else z._pari_init_() for z in coeffs[1:]])
     1463
     1464            # Tell the L-series / PARI about the coefficients.
     1465
     1466            if self.is_selfdual():
     1467                L.init_coeffs('a(k)', pari_precode = s)
     1468            else:
     1469                # Have to directly call gp_eval, since case of functional equation having two different
     1470                # (conjugate) L-functions isn't supported in Dokchitser class (yet).
     1471                L._Dokchitser__init = True
     1472                L._gp_eval(s)
     1473                L._gp_eval('initLdata("a(k)",1,"conj(a(k))")')
     1474
     1475            if epsilon == 'solve':
     1476                cmd = "sgneq = Vec(checkfeq()); sgn = -sgneq[2]/sgneq[1]; sgn"
     1477                epsilon = ComplexField(prec)(L._gp_eval(cmd))
     1478                if abs(abs(epsilon)-1) > tiny0:
     1479                    raise RuntimeError, "unable to determine epsilon from functional equation working to precision %s, since we get epsilon=%s, which is not sufficiently close to 1"%(prec, epsilon)
     1480                # 1, -1 are two common special cases, where it is clear what the
     1481                # infinite precision version is.
     1482                if epsilon == 1:
     1483                    self._epsilon = 1
     1484                elif epsilon == -1:
     1485                    self._epsilon = -1
     1486                else:
     1487                    self._epsilon = epsilon
     1488
     1489            fe = L.check_functional_equation()
     1490            if abs(fe) <= tiny0:
     1491                # one worked!
     1492                self._anlist[prec] = coeffs
     1493                return L
     1494            else:
     1495                pass
     1496           
     1497        # They all failed.
     1498        raise RuntimeError, "invalid L-series parameters: functional equation not satisfied"
     1499
     1500    def check_functional_equation(self, T, prec=53):
     1501        return self._function(prec=prec).check_functional_equation(T)
     1502
     1503class LSeriesProductEvaluator(object):
     1504    def __init__(self, factorization, prec):
     1505        self._factorization = factorization
     1506        self._prec = prec
     1507
     1508    def __call__(self, s):
     1509        try:
     1510            v = self._functions
     1511        except AttributeError:
     1512            self._functions = [(L._function(self._prec),e) for L,e in self._factorization]
     1513            v = self._functions
     1514        return prod(f(s)**e for f,e in v)
     1515
     1516class LSeriesProduct(object):
     1517    """
     1518    A formal product of L-series.
     1519    """
     1520    def __init__(self, F):
     1521        """
     1522        INPUT:
     1523            - `F` -- list of pairs (L,e) where L is an L-function and e is a nonzero integer.
     1524        """
     1525        if not isinstance(F, Factorization):
     1526            F = Factorization(F)
     1527        F.sort(cmp)
     1528        if len(F) == 0:
     1529            raise ValueError, "product must be nonempty"
     1530        self._factorization = F
     1531
     1532    def __cmp__(self, right):
     1533        # TODO: make work even if right not a product
     1534        return cmp(self._factorization, right._factorization)
     1535
     1536    def is_selfdual(self):
     1537        """
     1538        Return True if every factor of self is self dual; otherwise, return False.
     1539       
     1540        EXAMPLES::
     1541
     1542            sage: from sage.lfunctions.eulerprod import LSeries
     1543            sage: L1 = LSeries('zeta'); L1.is_selfdual()
     1544            True
     1545            sage: L2 = LSeries(DirichletGroup(9).0); L2.is_selfdual()
     1546            False
     1547            sage: (L1*L1).is_selfdual()
     1548            True
     1549            sage: (L1*L2).is_selfdual()
     1550            False
     1551            sage: (L2*L2).is_selfdual()
     1552            False
     1553        """
     1554        is_selfdual = set(L.is_selfdual() for L,e in self._factorization)
     1555        if len(is_selfdual) > 1:
     1556            return False
     1557        else:
     1558            return list(is_selfdual)[0]
     1559
     1560    def conductor(self):
     1561        """
     1562        Return the conductor of this product, which we define to be
     1563        the product with multiplicities of the conductors of the
     1564        factors of self.
     1565
     1566        EXAMPLES::
     1567
     1568            sage: from sage.lfunctions.eulerprod import LSeries
     1569            sage: J = J0(33); L = LSeries(J)
     1570            sage: L.conductor()
     1571            3993
     1572            sage: L.conductor().factor()
     1573            3 * 11^3
     1574            sage: J.decomposition()
     1575            [
     1576            Simple abelian subvariety 11a(1,33) of dimension 1 of J0(33),
     1577            Simple abelian subvariety 11a(3,33) of dimension 1 of J0(33),
     1578            Simple abelian subvariety 33a(1,33) of dimension 1 of J0(33)
     1579            ]
     1580
     1581        Of course, the conductor as we have defined it need not be an integer::
     1582
     1583            sage: L = LSeries(J[0])/LSeries(J[2])^2; L
     1584            (L-series of a degree 1 newform of level 11 and weight 2) * (L-series of a degree 1 newform of level 33 and weight 2)^-2
     1585            sage: L.conductor()
     1586            1/99
     1587        """
     1588        return prod(L.conductor()**e for L, e in self._factorization)
     1589
     1590    def __pow__(self, n):
     1591        """
     1592        Return the n-th power of this formal L-series product, where n can be any integer.
     1593       
     1594        EXAMPLES::
     1595
     1596            sage: from sage.lfunctions.eulerprod import LSeries
     1597            sage: L = LSeries('zeta') * LSeries(DirichletGroup(9).0)
     1598            sage: L(2)
     1599            1.80817853715812 + 0.369298119218816*I
     1600            sage: L = LSeries('zeta') * LSeries(DirichletGroup(9).0)
     1601            sage: L3 = L^3; L3
     1602            (Riemann Zeta function viewed as an L-series)^3 * (L-series attached to Dirichlet character modulo 9 of conductor 9 mapping 2 |--> zeta6)^3
     1603            sage: L3(2)
     1604            5.17205298762567 + 3.57190597873829*I
     1605            sage: CC((zeta(2)*LSeries(DirichletGroup(9).0)(2))^3)
     1606            5.17205298762567 + 3.57190597873829*I
     1607            sage: L^(-1999)
     1608            (Riemann Zeta function viewed as an L-series)^-1999 * (L-series attached to Dirichlet character modulo 9 of conductor 9 mapping 2 |--> zeta6)^-1999
     1609            sage: (L^(-1999)) (2)
     1610            8.90248311986228e-533 - 6.20123089437732e-533*I
     1611        """
     1612        return LSeriesProduct(self._factorization**ZZ(n))
     1613
     1614    def __mul__(self, right):
     1615        if isinstance(right, LSeriesAbstract):
     1616            return LSeriesProduct(self._factorization * Factorization([(right,1)]))
     1617        elif isinstance(right, LSeriesProduct):
     1618            return LSeriesProduct(self._factorization * right._factorization)
     1619        else:
     1620            raise TypeError
     1621
     1622    def __div__(self, right):
     1623        if isinstance(right, LSeriesAbstract):
     1624            return LSeriesProduct(self._factorization * Factorization([(right,-1)]))
     1625        elif isinstance(right, LSeriesProduct):
     1626            return LSeriesProduct(self._factorization / right._factorization)
     1627        raise TypeError
     1628
     1629    def parent(self):
     1630        return LSeriesParent
     1631   
     1632    def factor(self):
     1633        return self._factorization
     1634
     1635    def hodge_numbers(self):
     1636        """
     1637        Return the Hodge numbers of this product of L-series.
     1638        """
     1639
     1640    def degree(self):
     1641        return sum(e*L.degree() for L,e in self._factorization)
     1642
     1643    def local_factor(self, P, prec=None):
     1644        return prod(L.local_factor(P,prec)**e for L,e in self._factorization)
     1645
     1646    def __getitem__(self, *args, **kwds):
     1647        return self._factorization.__getitem__(*args, **kwds)
     1648
     1649    def __repr__(self):
     1650        return self.factor().__repr__()
     1651
     1652    def __call__(self, s):
     1653        return self._function(prec(s))(s)
     1654
     1655    def derivative(self, k=1):
     1656        raise NotImplementedError
     1657
     1658    def taylor_series(self, center=None, degree=6, variable='z', prec=53):
     1659        """
     1660        EXAMPLE::
     1661
     1662            sage: from sage.lfunctions.eulerprod import LSeries
     1663            sage: L1 = LSeries('zeta'); L2 = LSeries('delta')
     1664            sage: f = L1 * L2; f
     1665            (Riemann Zeta function viewed as an L-series) * (L-function associated to Ramanujan's Delta (a weight 12 cusp form))
     1666            sage: f.taylor_series(center=2, degree=4, variable='w', prec=30)
     1667            0.24077647 + 0.10066485*w + 0.061553731*w^2 - 0.041923238*w^3 + O(w^4)
     1668            sage: L1.taylor_series(2, 4, 'w', 30) * L2.taylor_series(2, 4, 'w', 30)
     1669            0.24077647 + 0.10066485*w + 0.061553731*w^2 - 0.041923238*w^3 + O(w^4)
     1670            sage: f = L1 / L2; f
     1671            (Riemann Zeta function viewed as an L-series) * (L-function associated to Ramanujan's Delta (a weight 12 cusp form))^-1
     1672            sage: f.taylor_series(center=2, degree=4, variable='w', prec=30)
     1673            11.237843 - 17.508629*w + 21.688182*w^2 - 24.044641*w^3 + O(w^4)
     1674            sage: L1.taylor_series(2, 4, 'w', 30) / L2.taylor_series(2, 4, 'w', 30)
     1675            11.237843 - 17.508629*w + 21.688182*w^2 - 24.044641*w^3 + O(w^4)
     1676           
     1677        """
     1678        return prod(L.taylor_series(center, degree, variable, prec)**e for L,e in self._factorization)
     1679       
     1680    def analytic_rank(self, prec=53):
     1681        """
     1682        Return sum of the order of vanishing counted with
     1683        multiplicities of each factors at their center point.
     1684
     1685        WARNING: The analytic rank is computed numerically, so is
     1686        definitely not provably correct.
     1687
     1688        EXAMPLES::
     1689
     1690        We compute the analytic rank of the non-simple 32-dimensional modular abelian variety `J_0(389)`.
     1691
     1692            sage: from sage.lfunctions.eulerprod import LSeries
     1693            sage: M = ModularSymbols(389,sign=1).cuspidal_subspace()
     1694            sage: L = LSeries(M); L
     1695            L-series attached to Modular Symbols subspace of dimension 32 of Modular Symbols space of dimension 33 for Gamma_0(389) of weight 2 with sign 1 over Rational Field
     1696
     1697        We first attempt computation of the analytic rank with the default of 53 bits precision::
     1698       
     1699            sage: L.analytic_rank()
     1700            Traceback (most recent call last):
     1701            ...
     1702            RuntimeError: invalid L-series parameters: functional equation not satisfied
     1703
     1704        The above failed because trying to compute one of the degree
     1705        20 newforms resulting in some internal error when double
     1706        checking the functional equation.  So we try with slightly more precision::
     1707
     1708            sage: L.analytic_rank(70)
     1709            13
     1710
     1711        This works, since the factors have dimensions 1,2,3,6,20, and
     1712        the one of degree 1 has rank 2, the ones of degree 2,3,6 have
     1713        rank 2,3,6, respectively, and the one of degree 20 has rank 0::
     1714
     1715            sage: 2*1 + 2 + 3 + 6
     1716            13
     1717        """
     1718        return sum(e*L.analytic_rank(prec=prec) for L,e in self._factorization)
     1719
     1720    def weight(self):
     1721        """
     1722        Return the weight of this L-series, which is the sum of the weights
     1723        of the factors counted with multiplicity.
     1724        """
     1725        return sum(e*L.weight() for L,e in self._factorization)
     1726
     1727    @cached_method
     1728    def _function(self, prec=53):
     1729        """
     1730        Return Dokchitser object that allows for computation of this
     1731        L-series.  This is used behind the scenes for evaluation and
     1732        computation of Taylor series.
     1733        """
     1734        return LSeriesProductEvaluator(self._factorization, prec)
     1735       
     1736       
     1737class LSeriesZeta(LSeriesAbstract):
     1738    """
     1739    EXAMPLES::
     1740
     1741        sage: from sage.lfunctions.eulerprod import LSeries
     1742        sage: L = LSeries('zeta'); L
     1743        Riemann Zeta function viewed as an L-series
     1744        sage: L(2)
     1745        1.64493406684823
     1746        sage: L(2.000000000000000000000000000000000000000000000000000)
     1747        1.64493406684822643647241516664602518921894990120680
     1748        sage: zeta(2.000000000000000000000000000000000000000000000000000)
     1749        1.64493406684822643647241516664602518921894990120680
     1750        sage: L.local_factor(3)
     1751        -T + 1
     1752        sage: L.local_factor(5)
     1753        -T + 1
     1754        sage: L.anlist(30)
     1755        [0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1]       
     1756    """
     1757    def __init__(self):
     1758        LSeriesAbstract.__init__(self, conductor=1, hodge_numbers=[0], weight=1, epsilon=1,
     1759                                 poles=[0,1], residues=[9,8], base_field=QQ)
     1760        T = ZZ['T'].gen()
     1761        self._lf = 1 - T
     1762
     1763    def _cmp(self, right):
     1764        return 0
     1765
     1766    def _local_factor(self, P, prec):
     1767        return self._lf
     1768
     1769    def __repr__(self):
     1770        return "Riemann Zeta function viewed as an L-series"
     1771   
     1772class LSeriesDelta(LSeriesAbstract):
     1773    """
     1774    EXAMPLES::
     1775
     1776        sage: from sage.lfunctions.eulerprod import LSeriesDelta; L = LSeriesDelta()
     1777        sage: L.anlist(10)
     1778        [0, 1, -24, 252, -1472, 4830, -6048, -16744, 84480, -113643, -115920]
     1779        sage: list(delta_qexp(11))
     1780        [0, 1, -24, 252, -1472, 4830, -6048, -16744, 84480, -113643, -115920]
     1781        sage: L.anlist(10^4) == list(delta_qexp(10^4+1))
     1782        True
     1783    """
     1784    def __init__(self):
     1785        LSeriesAbstract.__init__(self, conductor=1, hodge_numbers=[0,1], weight=12, epsilon=1,
     1786                                 poles=[], residues=[], base_field=QQ)
     1787        self._T = ZZ['T'].gen()
     1788        self._lf = {}
     1789
     1790    def _local_factor(self, P, prec):
     1791        """
     1792        EXAMPLES::
     1793
     1794            sage: from sage.lfunctions.eulerprod import LSeriesDelta; L = LSeriesDelta()
     1795            sage: L.local_factor(2)
     1796            2048*T^2 + 24*T + 1
     1797            sage: L._local_factor(11, None)   # really this is called
     1798            285311670611*T^2 - 534612*T + 1
     1799
     1800        The representation is reducible modulo 691::
     1801       
     1802            sage: L.local_factor(2).factor_mod(691)
     1803            (666) * (T + 387) * (T + 690)
     1804            sage: L.local_factor(3).factor_mod(691)
     1805            (251) * (T + 234) * (T + 690)
     1806            sage: L.local_factor(11).factor_mod(691)
     1807            (468) * (T + 471) * (T + 690)
     1808
     1809        ... because of the 691 here::
     1810
     1811            sage: bernoulli(12)
     1812            -691/2730
     1813        """
     1814        try:
     1815            return self._lf[P]
     1816        except KeyError:
     1817            pass
     1818        self._precompute_local_factors(P+1)
     1819        return self._lf[P]
     1820
     1821    def _precompute_local_factors(self, bound, prec=None):
     1822        """
     1823        Precompute local factors up to the given bound.
     1824       
     1825        EXAMPLES::
     1826
     1827            sage: from sage.lfunctions.eulerprod import LSeriesDelta; L = LSeriesDelta()
     1828            sage: L._lf
     1829            {}
     1830            sage: L._precompute_local_factors(10)
     1831            sage: L._lf
     1832            {2: 2048*T^2 + 24*T + 1, 3: 177147*T^2 - 252*T + 1, 5: 48828125*T^2 - 4830*T + 1, 7: 1977326743*T^2 + 16744*T + 1}
     1833        """
     1834        from sage.modular.all import delta_qexp
     1835        T = self._T
     1836        T2 = T**2
     1837        f = delta_qexp(bound)
     1838        for p in prime_range(bound):
     1839            if not self._lf.has_key(p):
     1840                self._lf[p] = 1 - f[p]*T + (p**11)*T2
     1841
     1842    def __repr__(self):
     1843        return "L-function associated to Ramanujan's Delta (a weight 12 cusp form)"
     1844
     1845    def _cmp(self, right):
     1846        return 0
     1847
     1848                 
     1849class LSeriesEllipticCurve(LSeriesAbstract):
     1850    def __init__(self, E, prec=53):
     1851        """
     1852        EXAMPLES::
     1853            sage: from sage.lfunctions.eulerprod import LSeriesEllipticCurve
     1854            sage: L = LSeriesEllipticCurve(EllipticCurve('389a'))
     1855            sage: L(2)
     1856            0.360092863578881
     1857        """
     1858        E = E.global_minimal_model()
     1859        self._E = E
     1860        K = E.base_field()
     1861        d = K.degree()
     1862        self._N = E.conductor()
     1863        LSeriesAbstract.__init__(self, conductor = norm(self._N) * K.discriminant()**2,
     1864                                 hodge_numbers = [0]*d+[1]*d, weight = 2, epsilon = [1,-1],
     1865                                 poles = [], residues=[], base_field = K, prec=prec)
     1866
     1867    def elliptic_curve(self):
     1868        return self._E
     1869
     1870    def _cmp(self, right):
     1871        return cmp(self.elliptic_curve(), right.elliptic_curve())
     1872
     1873    def __repr__(self):
     1874        return "L-series of %s"%self.elliptic_curve()
     1875
     1876    def _local_factor(self, P, prec):
     1877        R = ZZ['T']
     1878        T = R.gen()
     1879        q = norm(P)
     1880        p = prime_below(P)
     1881        f = ZZ(q).ord(p)
     1882        if P.divides(self._N):
     1883            a = self._E.local_data(P).bad_reduction_type()
     1884            return 1 - a*(T**f)
     1885        else:
     1886            a = q + 1 - self._E.reduction(P).count_points()
     1887            return 1 - a*(T**f) + q*(T**(2*f))
     1888
     1889       
     1890class LSeriesEllipticCurveQQ(LSeriesEllipticCurve):
     1891    def __init__(self, E):
     1892        E = E.global_minimal_model()
     1893        self._E = E
     1894        K = E.base_field()
     1895        self._N = E.conductor()
     1896        self._lf = {}
     1897        self._T = ZZ['T'].gen()
     1898        LSeriesAbstract.__init__(self, conductor = self._N,
     1899                                 hodge_numbers = [0,1], weight = 2, epsilon = E.root_number(),
     1900                                 poles = [], residues=[], base_field = QQ)
     1901
     1902    def _lf0(self, p):
     1903        a = self._E.ap(p)
     1904        T = self._T
     1905        if self._N%p == 0:
     1906            if self._N%(p*p) == 0:
     1907                return T.parent()(1)
     1908            else:
     1909                return 1 - a*T
     1910        else:
     1911            return 1 - a*T + p*T*T
     1912
     1913    def _precompute_local_factors(self, bound, prec=None):
     1914        for p in primes(bound):
     1915            if not self._lf.has_key(p):
     1916                self._lf[p] = self._lf0(p)
     1917
     1918    def _local_factor(self, P, prec):
     1919        if self._lf.has_key(P):
     1920            return self._lf[P]
     1921        else:
     1922            return self._lf0(P)
     1923
     1924    def _primes_above(self, p):
     1925        return [p]
     1926
     1927class LSeriesDedekindZeta(LSeriesAbstract):
     1928    """
     1929    EXAMPLES::
     1930
     1931        sage: from sage.lfunctions.eulerprod import LSeries
     1932        sage: K.<a> = NumberField(x^3 - 2)
     1933        sage: L = LSeries(K); L
     1934        Dedekind Zeta function of Number Field in a with defining polynomial x^3 - 2
     1935        sage: L(2)
     1936        1.60266326190044
     1937        sage: L.residues()
     1938        'automatic'
     1939        sage: L.residues(prec=53)
     1940        [-4.77632833933856]
     1941        sage: L.residues(prec=100)
     1942        [-4.7763283393385594030639875094]
     1943    """
     1944    def __init__(self, K):
     1945        if not K.is_absolute():
     1946            K = K.absolute_field(names='a')
     1947        self._K = K
     1948        d = K.degree()
     1949        sigma = K.signature()[1]
     1950        LSeriesAbstract.__init__(self,
     1951                                 conductor = abs(K.discriminant()),
     1952                                 hodge_numbers = [0]*(d-sigma) + [1]*sigma,
     1953                                 weight = 1,
     1954                                 epsilon = 1,
     1955                                 poles = [1],
     1956                                 residues = 'automatic',
     1957                                 base_field = K,
     1958                                 is_selfdual = True)
     1959        self._T = ZZ['T'].gen()
     1960
     1961    def _cmp(self, right):
     1962        return cmp(self.number_field(), right.number_field())
     1963
     1964    def number_field(self):
     1965        return self._K
     1966   
     1967    def __repr__(self):
     1968        return "Dedekind Zeta function of %s"%self._K
     1969
     1970    def _local_factor(self, P, prec):
     1971        T = self._T
     1972        return 1 - T**P.residue_class_degree()
     1973   
     1974
     1975class LSeriesDirichletCharacter(LSeriesAbstract):
     1976    """
     1977    EXAMPLES::
     1978   
     1979        sage: from sage.lfunctions.eulerprod import LSeries;  L = LSeries(DirichletGroup(5).0)
     1980        sage: L(3)
     1981        0.988191681624057 + 0.0891051883457395*I
     1982    """
     1983    def __init__(self, chi):
     1984        if not chi.is_primitive():
     1985            raise NotImplementedError, "chi must be primitive"
     1986        if chi.is_trivial():
     1987            raise NotImplementedError, "chi must be nontrivial"
     1988        if chi.base_ring().characteristic() != 0:
     1989            raise ValueError, "base ring must have characteristic 0"
     1990        self._chi = chi
     1991        LSeriesAbstract.__init__(self, conductor = chi.conductor(),
     1992                                 hodge_numbers = [1] if chi.is_odd() else [0],
     1993                                 weight = 1,
     1994                                 epsilon = None,
     1995                                 poles = [],  # since primitive
     1996                                 residues = [],  # since primitive
     1997                                 base_field = QQ,
     1998                                 is_selfdual = chi.order() <= 2)
     1999        self._T = ZZ['T'].gen()
     2000
     2001    def _cmp(self, right):
     2002        return cmp(self.character(), right.character())
     2003
     2004    def __repr__(self):
     2005        """
     2006        EXAMPLES::
     2007
     2008            sage: from sage.lfunctions.eulerprod import LSeries;  L = LSeries(DirichletGroup(3).0)
     2009            sage: L.__repr__()
     2010            'L-series attached to Dirichlet character modulo 3 of conductor 3 mapping 2 |--> -1'
     2011        """
     2012        return "L-series attached to %s"%self._chi
     2013
     2014    def character(self):
     2015        return self._chi
     2016
     2017    def epsilon(self, prec=None):
     2018        chi = self._chi
     2019        if prec is None:
     2020            # answer in symbolic ring
     2021            return (sqrt(-1) * SR(chi.gauss_sum())) / chi.modulus().sqrt()
     2022        else:
     2023            C = ComplexField(prec)
     2024            x = C(chi.modulus()).sqrt() / chi.gauss_sum_numerical(prec=prec)
     2025            if chi.is_odd():
     2026                x *= C.gen()
     2027            return 1/x
     2028       
     2029    def _local_factor(self, P, prec):
     2030        a = self._chi(P)
     2031        if prec is not None:
     2032            a = ComplexField(prec)(a)
     2033        return 1 - a*self._T
     2034
     2035class LSeriesModularSymbolsAbstract(LSeriesAbstract):
     2036    def _cmp(self, right):
     2037        return cmp(self.modular_symbols(), right.modular_symbols())
     2038       
     2039    def modular_symbols(self):
     2040        return self._M
     2041   
     2042    def __repr__(self):
     2043        """
     2044        EXAMPLES::
     2045
     2046            sage: from sage.lfunctions.eulerprod import LSeries
     2047            sage: f = Newforms(43,2,names='a')[1]; f
     2048            q + a1*q^2 - a1*q^3 + (-a1 + 2)*q^5 + O(q^6)
     2049            sage: LSeries(f).__repr__()
     2050            'L-series of a degree 2 newform of level 43 and weight 2'
     2051        """
     2052        return "L-series of a degree %s newform of level %s and weight %s"%(self._M.dimension(), self._M.level(), self._M.weight())
     2053   
     2054    def _precompute_local_factors(self, bound, prec):
     2055        primes = [p for p in prime_range(bound) if not self._lf.has_key(p) or self._lf[p][0] < prec]
     2056        self._do_precompute(primes, prec)
     2057
     2058    def _do_precompute(self, primes, prec):
     2059        E, v = self._M.compact_system_of_eigenvalues(primes)
     2060        if prec == 53:
     2061            C = CDF
     2062        elif prec is None or prec==oo:
     2063            if v.base_ring() == QQ:
     2064                C = QQ
     2065            else:
     2066                C = CDF
     2067        else:
     2068            C = ComplexField(prec)
     2069           
     2070        phi = v.base_ring().embeddings(C)[self._conjugate]
     2071        v = vector(C, [phi(a) for a in v])
     2072        aplist = E.change_ring(C) * v
     2073        T = C['T'].gen(); T2 = T**2
     2074        chi = self._M.character()
     2075        k = self.weight()
     2076        for i in range(len(primes)):
     2077            p = primes[i]
     2078            s = chi(p)
     2079            if s != 0: s *= p**(k-1)
     2080            F = 1 - aplist[i]*T + s*T2
     2081            self._lf[p] = (prec, F)
     2082       
     2083    def _local_factor(self, P, prec):
     2084        # TODO: ugly -- get rid of all "prec=None" in whole program -- always use oo.
     2085        if prec is None: prec = oo
     2086        if self._lf.has_key(P) and self._lf[P][0] >= prec:
     2087            return self._lf[P][1]
     2088        else:
     2089            self._do_precompute([P],prec)
     2090            return self._lf[P][1]           
     2091
     2092class LSeriesModularSymbolsNewformGamma0(LSeriesModularSymbolsAbstract):
     2093    def _cmp(self, right):
     2094        return cmp((self._M, self._conjugate), (right._M, right._conjugate))
     2095
     2096    def __init__(self, M, conjugate=0, check=True, epsilon=None):
     2097        """
     2098        INPUT:
     2099            - M -- a simple, new, cuspidal modular symbols space with
     2100              sign 1
     2101            - conjugate -- (default: 0), integer between 0 and dim(M)-1
     2102            - check -- (default: True), if True, checks that M is
     2103              simple, new, cuspidal, which can take a very long time,
     2104              depending on how M was defined
     2105            - epsilon -- (default: None), if not None, should be the sign
     2106              in the functional equation, which is -1 or 1.  If this is
     2107              None, then epsilon is computed by computing the sign of
     2108              the main Atkin-Lehner operator on M.  If you have a faster
     2109              way to determine epsilon, use it.
     2110
     2111        EXAMPLES::
     2112
     2113            sage: from sage.lfunctions.eulerprod import LSeriesModularSymbolsNewformGamma0
     2114            sage: M = ModularSymbols(43,sign=1).cuspidal_subspace()[1]; M
     2115            Modular Symbols subspace of dimension 2 of Modular Symbols space of dimension 4 for Gamma_0(43) of weight 2 with sign 1 over Rational Field
     2116            sage: L0 = LSeriesModularSymbolsNewformGamma0(M,0,check=False,epsilon=1); L0
     2117            L-series of a degree 2 newform of level 43 and weight 2
     2118            sage: L0.taylor_series()
     2119            0.620539857407845 + 0.331674007376949*z - 0.226392184536589*z^2 + 0.0960519649929789*z^3 - 0.00451826124421802*z^4 - 0.0203363026933833*z^5 + O(z^6)
     2120            sage: L1 = LSeriesModularSymbolsNewformGamma0(M,1,check=False,epsilon=1); L1
     2121            L-series of a degree 2 newform of level 43 and weight 2
     2122            sage: L1(1)
     2123            0.921328017272472
     2124            sage: L1.taylor_series()
     2125            0.921328017272472 + 0.492443075089339*z - 0.391019352704047*z^2 + 0.113271812405127*z^3 + 0.0213067052584679*z^4 - 0.0344198080536274*z^5 + O(z^6)
     2126        """
     2127        if M.dimension() == 0:
     2128            raise ValueError, "modular symbols space must positive dimension"
     2129        chi = M.character()
     2130        if chi is None or not chi.is_trivial():
     2131            raise ValueError, "modular symbols space must have trivial character"
     2132        self._M = M
     2133        N = M.level()
     2134
     2135        if check:
     2136            if not M.is_simple():
     2137                raise ValueError, "modular symbols space must be simple"
     2138            if not M.is_new():
     2139                raise ValueError, "modular symbols space must be new"
     2140            if not M.is_cuspidal():
     2141                raise ValueError, "modular symbols space must be cuspidal"
     2142
     2143        k = M.weight()
     2144        if epsilon is None:
     2145            w = M.atkin_lehner_operator(N).matrix()
     2146            if w not in [-1, 1]:
     2147                raise ValueError, "modular symbols space must have constant Atkin-Lehner operator"
     2148            epsilon = (-1)**(k/2) * w[0,0]
     2149
     2150        conjugate = ZZ(conjugate)
     2151        if conjugate < 0 or conjugate >= M.dimension():
     2152            raise ValueError, "conjugate must a nonnegative integer less than the dimension"
     2153        self._conjugate = conjugate
     2154        self._lf = {}
     2155       
     2156        LSeriesAbstract.__init__(self,
     2157                                 conductor = N,
     2158                                 hodge_numbers = [0,1],
     2159                                 weight = k,
     2160                                 epsilon = epsilon,
     2161                                 poles = [],  # since primitive
     2162                                 residues = [],  # since primitive
     2163                                 base_field = QQ)
     2164
     2165def _is_valid_modsym_space(M):
     2166    if not is_ModularSymbolsSpace(M):
     2167        raise TypeError, "must be a modular symbols space"
     2168    if M.dimension() == 0:
     2169        raise ValueError, "modular symbols space must positive dimension"
     2170    if M.character() is None:
     2171        raise ValueError, "modular symbols space must have associated character"
     2172    if not M.is_simple():
     2173        raise ValueError, "modular symbols space must be simple"
     2174    if not M.is_new():
     2175        raise ValueError, "modular symbols space must be new"
     2176    if not M.is_cuspidal():
     2177        raise ValueError, "modular symbols space must be cuspidal"
     2178   
     2179
     2180class LSeriesModularSymbolsNewformCharacter(LSeriesModularSymbolsAbstract):
     2181    def _cmp(self, right):
     2182        return cmp((self._M, self._conjugate), (right._M, right._conjugate))
     2183   
     2184    def __init__(self, M, conjugate=0):
     2185        _is_valid_modsym_space(M)
     2186        chi = M.character()
     2187        self._M = M
     2188        N = M.level()
     2189
     2190        # See Remark 5.0.2 in [Diamond-Im] which says: "Let f be a newform of level N
     2191        # which is a common eigenform under all the Hecke operators T_p.  Then
     2192        # w_N(f) = c*fbar, where fbar = sum bar(a_n) q^n and c is a scalar.  The functional
     2193        # equation may be rewritten as   Lambda(s, f) = c * i^k * Lambda(k-s, fbar).
     2194        # That said, this seems hard to compute, so we just solve using the
     2195        # functional equation.
     2196        epsilon = 'solve'
     2197       
     2198        k = M.weight()
     2199        conjugate = ZZ(conjugate)
     2200        if conjugate < 0 or conjugate >= M.dimension():
     2201            raise ValueError, "conjugate must a nonnegative integer less than the dimension"
     2202        self._conjugate = conjugate
     2203
     2204
     2205        LSeriesAbstract.__init__(self,
     2206                                 conductor = N,
     2207                                 hodge_numbers = [0,1],
     2208                                 weight = k,
     2209                                 epsilon = epsilon,
     2210                                 poles = [],  # since primitive
     2211                                 residues = [],  # since primitive
     2212                                 base_field = QQ,
     2213                                 is_selfdual = chi.order() <= 2)
     2214        self._lf = {}
     2215
     2216def _new_modsym_space_with_multiplicity(M):
     2217    """
     2218    Returns a simple new modular symbols space N and an integer d such
     2219    that M is isomorphic to `N^d` as a module over the anemic Hecke
     2220    algebra.
     2221   
     2222    INPUT:
     2223        - M -- a sign=1 modular simple space for the full Hecke
     2224          algebra (including primes dividing the level) that can't be
     2225          decomposed further by the Hecke operators.  None of the
     2226          conditions on M are explicitly checked.
     2227       
     2228    OUTPUT:
     2229        - N -- a simple new modular symbols space
     2230        - d -- a positive integer
     2231    """
     2232    if M.is_new():
     2233        return [(M,1)]
     2234    raise NotImplementedError
     2235
     2236def LSeriesModularSymbolsNewform(M, i=0):
     2237    chi = M.character()
     2238    if chi is None:
     2239        raise NotImplementedError
     2240    elif chi.is_trivial():
     2241        return LSeriesModularSymbolsNewformGamma0(M, i)
     2242    else:
     2243        return LSeriesModularSymbolsNewformCharacter(M, i)
     2244
     2245class LSeriesModularSymbolsMotive(LSeriesProduct):
     2246    """
     2247    The product of L-series attached to the modular symbols space M.
     2248    """
     2249    def __init__(self, M):
     2250        self._M = M
     2251        if not is_ModularSymbolsSpace(M):
     2252            raise TypeError, "X must be a modular symbols space or have a modular symbols method"
     2253        self._M = M
     2254        D = M.decomposition()
     2255        for A in D:
     2256            _is_valid_modsym_space(A)
     2257        F = []
     2258        for A in D:
     2259            for X in _new_modsym_space_with_multiplicity(A):
     2260                N, d = X
     2261                chi = N.character()
     2262                for i in range(N.dimension()):
     2263                    F.append( (LSeriesModularSymbolsNewform(N,i), d))
     2264        LSeriesProduct.__init__(self, F)
     2265
     2266    def modular_symbols(self):
     2267        return self._M
     2268
     2269    def __repr__(self):
     2270        return "L-series attached to %s"%self._M
     2271   
     2272class LSeriesModularAbelianVariety(LSeriesProduct):
     2273    """
     2274    The product of L-series attached to the modular abelian variety A.
     2275
     2276    EXAMPLES::
     2277
     2278        sage: from sage.lfunctions.eulerprod import LSeries
     2279        sage: L = LSeries(J0(54)); L
     2280        L-series attached to Abelian variety J0(54) of dimension 4
     2281        sage: L.factor()
     2282        (L-series of a degree 1 newform of level 27 and weight 2)^2 * (L-series of a degree 1 newform of level 54 and weight 2) * (L-series of a degree 1 newform of level 54 and weight 2)
     2283        sage: L(1)
     2284        0.250717238804658
     2285        sage: L.taylor_series(prec=20)
     2286        0.25072 + 0.59559*z + 0.15099*z^2 - 0.35984*z^3 + 0.056934*z^4 + 0.17184*z^5 + O(z^6)
     2287
     2288    Independent check of L(1)::
     2289   
     2290        sage: prod(EllipticCurve(lbl).lseries()(1) for lbl in ['54a', '54b', '27a', '27a'])
     2291        0.250717238804658
     2292
     2293    Different check that totally avoids using Dokchitser::
     2294   
     2295        sage: prod(EllipticCurve(lbl).lseries().at1()[0] for lbl in ['54a', '54b', '27a', '27a'])
     2296        0.250848605530185
     2297    """
     2298    def __init__(self, A):
     2299        self._A = A
     2300        D = A.decomposition()
     2301        F = None
     2302        for A in D:
     2303            # TODO: This is ugly, but I don't know a cleaner way to do it yet.
     2304            # Could be painfully inefficient in general.
     2305            f = Newform(A.newform_label(), names='a')
     2306            M = f.modular_symbols(sign=1)
     2307            d = ZZ(A.dimension() / M.dimension())
     2308            L = LSeriesModularSymbolsMotive(M)**d
     2309            if F is None:
     2310                F = L
     2311            else:
     2312                F *= L
     2313        if F is None:
     2314            raise ValueError, "abelian variety must have positive dimension"
     2315        LSeriesProduct.__init__(self, F.factor())
     2316
     2317    def abelian_variety(self):
     2318        return self._A
     2319
     2320    def __repr__(self):
     2321        return "L-series attached to %s"%self._A
     2322
     2323
     2324class LSeriesTwist(LSeriesAbstract):
     2325    """
     2326    Twist of an L-series by a character.
     2327    """
     2328    def __init__(self, L, chi, conductor=None, epsilon=None, prec=53):
     2329        """
     2330        INPUT:
     2331            - `L` -- an L-series
     2332            - ``chi`` -- a character of the base field of L
     2333            - ``conductor`` -- None, or a list of conductors to try
     2334            - ``prec`` -- precision to use when trying conductors, if
     2335              conductor is a list
     2336        """
     2337        self._L = L
     2338        self._chi = chi
     2339
     2340        if not chi.is_primitive():
     2341            raise ValueError, "character must be primitive"
     2342       
     2343        A = ZZ(L.conductor())
     2344        B = chi.conductor()
     2345        if conductor is None:
     2346            if A.gcd(B) != 1:
     2347                # Make a list of all possible conductors, and let the
     2348                # functional equation figure it out.
     2349                smallest = ZZ(A)
     2350                while smallest.gcd(B) != 1:
     2351                    smallest = smallest // smallest.gcd(B)
     2352                biggest = A * (B**L.degree())
     2353                assert biggest % smallest == 0
     2354                #
     2355                # TODO: improve this using the theorem stated
     2356                # on page 1 of http://wstein.org/papers/padictwist/
     2357                #
     2358                conductor = [smallest*d for d in divisors(biggest//smallest)]
     2359            else:
     2360                conductor = A * (B**L.degree())
     2361        hodge_numbers = L.hodge_numbers()
     2362        weight = L.weight()
     2363        if epsilon is None:
     2364            if L.epsilon() != 'solve':
     2365                if chi.order() <= 2:
     2366                    if A.gcd(B) == 1:
     2367                        epsilon = L.epsilon() * chi(-A)
     2368                    else:
     2369                        epsilon = [L.epsilon(), -L.epsilon()]
     2370                else:
     2371                    epsilon = 'solve'
     2372            else:
     2373                epsilon = 'solve'
     2374        is_selfdual = L.is_selfdual()
     2375        poles = [] # ???  TODO -- no clue here.
     2376        residues = 'automatic'
     2377        base_field = L.base_field()
     2378        LSeriesAbstract.__init__(self, conductor, hodge_numbers, weight, epsilon,
     2379                                 poles, residues, base_field, is_selfdual, prec)
     2380
     2381    def _local_factor(self, P, prec):
     2382        L0 = self._L.local_factor(P, prec)
     2383        chi = self._chi
     2384        T = L0.parent().gen()
     2385        c = chi(P)
     2386        if prec is not None:
     2387            c = ComplexField(prec)(c)
     2388        return L0(c*T)
     2389
     2390    def __repr__(self):
     2391        return "Twist of %s by %s"%(self._L, self._chi)
     2392
     2393    def _cmp(self, right):
     2394        return cmp((self._L, self._chi), (right._L, right._chi))
     2395
     2396    def untwisted_lseries(self):
     2397        return self._L
     2398
     2399    def twist_character(self):
     2400        return self._chi
     2401
     2402def LSeries(X, *args, **kwds):
     2403    """
     2404    Return the L-series of X, where X can be any of the following:
     2405
     2406        - elliptic curve over a number field (including QQ)
     2407        - Dirichlet character
     2408        - cuspidal newform
     2409        - new cuspidal modular symbols space -- need not be simple
     2410        - string: 'zeta' (Riemann Zeta function), 'delta'
     2411        - modular elliptic curve attached to Hilbert modular forms space
     2412
     2413    For convenience, if L is returned, then L._X is set to X.
     2414
     2415    EXAMPLES::
     2416
     2417    The Dedekind Zeta function of a number field::
     2418
     2419        sage: from sage.lfunctions.eulerprod import LSeries
     2420        sage: K.<a> = NumberField(x^2 + 1)
     2421        sage: L = LSeries(K); L
     2422        Dedekind Zeta function of Number Field in a with defining polynomial x^2 + 1
     2423        sage: L(2)
     2424        1.50670300992299       
     2425
     2426        sage: K.zeta_coefficients(100) == L.anlist(100)[1:]
     2427        True
     2428
     2429        sage: L = LSeries(ModularSymbols(43, weight=2,sign=1).cuspidal_subspace().decomposition()[1])
     2430        sage: L(1)
     2431        0.571720756464112
     2432        sage: L.factor()[0][0](1)
     2433        0.620539857407845
     2434        sage: L.factor()[1][0](1)
     2435        0.921328017272472
     2436        sage: L._X
     2437        Modular Symbols subspace of dimension 2 of Modular Symbols space of dimension 4 for Gamma_0(43) of weight 2 with sign 1 over Rational Field       
     2438        sage: L = LSeries(ModularSymbols(DirichletGroup(13).0^2, weight=2,sign=1).cuspidal_subspace())
     2439        sage: L(1)
     2440        0.298115272465799 - 0.0402203326076733*I
     2441
     2442    The L-series of a modular abelian variety with both new and old parts::
     2443
     2444        sage: L = LSeries(J0(33)); L
     2445        L-series attached to Abelian variety J0(33) of dimension 3
     2446        sage: L.factor()
     2447        (L-series of a degree 1 newform of level 11 and weight 2)^2 * (L-series of a degree 1 newform of level 33 and weight 2)       
     2448        sage: L.local_factor(2, prec=oo)
     2449        8*T^6 + 12*T^5 + 12*T^4 + 8*T^3 + 6*T^2 + 3*T + 1
     2450        sage: L(1)
     2451        0.0481553138900504
     2452
     2453    We check the above computation of L(1) via independent methods (and implementations)::
     2454   
     2455        sage: prod(EllipticCurve(lbl).lseries().at1()[0] for lbl in ['11a', '11a', '33a'])
     2456        0.0481135342926321
     2457        sage: prod(EllipticCurve(lbl).lseries()(1) for lbl in ['11a', '11a', '33a'])
     2458        0.0481553138900504
     2459
     2460    A nonsimple new modular symbols space of level 43::
     2461
     2462        sage: L = LSeries(ModularSymbols(43,sign=1).cuspidal_subspace())
     2463        sage: L
     2464        L-series attached to Modular Symbols subspace of dimension 3 of Modular Symbols space of dimension 4 for Gamma_0(43) of weight 2 with sign 1 over Rational Field
     2465        sage: L(1)
     2466        0
     2467        sage: L.taylor_series()
     2468        0.196399786632435*z + 0.314922741074845*z^2 - 0.0797083673829092*z^3 - 0.161630566287135*z^4 + 0.123939472976207*z^5 + O(z^6)
     2469        sage: L.factor()
     2470        (L-series of a degree 1 newform of level 43 and weight 2) * (L-series of a degree 2 newform of level 43 and weight 2) * (L-series of a degree 2 newform of level 43 and weight 2)       
     2471        sage: L.analytic_rank()
     2472        1
     2473        sage: D = ModularSymbols(43,sign=1).cuspidal_subspace().decomposition()
     2474        sage: L0 = LSeries(D[0]); L1 = LSeries(D[1])
     2475        sage: L0.taylor_series() * L1.taylor_series()
     2476        0.196399786632435*z + 0.314922741074845*z^2 - 0.0797083673829091*z^3 - 0.161630566287135*z^4 + 0.123939472976207*z^5 + O(z^6)
     2477        sage: L0.factor()
     2478        L-series of a degree 1 newform of level 43 and weight 2
     2479        sage: L1.factor()
     2480        (L-series of a degree 2 newform of level 43 and weight 2) * (L-series of a degree 2 newform of level 43 and weight 2)
     2481
     2482    """
     2483    L = _lseries(X, *args, **kwds)
     2484    L._X = X
     2485    return L
     2486
     2487def _lseries(X, *args, **kwds):
     2488    """
     2489    Helper function used by LSeries function.
     2490    """
     2491    if is_EllipticCurve(X):
     2492        K = X.base_ring()
     2493        if is_RationalField(K):
     2494            return LSeriesEllipticCurveQQ(X, *args, **kwds)
     2495        else:
     2496            return LSeriesEllipticCurve(X)
     2497       
     2498    if is_DirichletCharacter(X):
     2499        if X.is_trivial() and X.is_primitive():
     2500            return LSeriesZeta(*args, **kwds)
     2501        else:
     2502            return LSeriesDirichletCharacter(X, *args, **kwds)
     2503
     2504    if is_NumberField(X):
     2505        return LSeriesDedekindZeta(X, *args, **kwds)
     2506
     2507    if isinstance(X, sage.modular.modform.element.Newform):
     2508        return LSeriesModularSymbolsNewform(X.modular_symbols(sign=1), *args, **kwds)
     2509
     2510    if is_ModularSymbolsSpace(X):
     2511        if X.sign() != 1:
     2512            raise NotImplementedError
     2513        return LSeriesModularSymbolsMotive(X, *args, **kwds)
     2514
     2515    if isinstance(X, str):
     2516        y = X.lower()
     2517        if y == 'zeta':
     2518            return LSeriesZeta(*args, **kwds)
     2519        elif y == 'delta':
     2520            return LSeriesDelta(*args, **kwds)
     2521        else:
     2522            raise ValueError, 'unknown L-series "%s"'%y
     2523
     2524    if is_ModularAbelianVariety(X):
     2525        return LSeriesModularAbelianVariety(X, *args, **kwds)
     2526       
     2527    raise NotImplementedError
  • new file sage/lfunctions/eulerprod_fast.pyx

    diff --git a/sage/lfunctions/eulerprod_fast.pyx b/sage/lfunctions/eulerprod_fast.pyx
    new file mode 100644
    - +  
     1
     2from sage.stats.intlist cimport IntList
     3
     4from sage.all import prime_range
     5
     6def prime_powers_intlist(Py_ssize_t B):
     7    """
     8    Return IntList of the prime powers and corresponding primes, up to
     9    B.  The list is *not* in the usual order; instead the order is like
     10    this: 2,4,8,...,3,9,27,..., 5,25,..., etc.
     11
     12    INPUT:
     13
     14        - B -- positive integer
     15
     16    EXAMPLES::
     17   
     18        sage: from sage.lfunctions.eulerprod_fast import prime_powers_intlist
     19        sage: prime_powers_intlist(10)
     20        ([1, 2, 4, 8, 3, 9, 5, 7], [1, 2, 2, 2, 3, 3, 5, 7])
     21        sage: prime_powers_intlist(10)
     22        ([1, 2, 4, 8, 3, 9, 5, 7], [1, 2, 2, 2, 3, 3, 5, 7])
     23        sage: prime_powers_intlist(20)
     24        ([1, 2, 4, 8, 16 ... 7, 11, 13, 17, 19], [1, 2, 2, 2, 2 ... 7, 11, 13, 17, 19])
     25        sage: list(sorted(prime_powers_intlist(10^6)[0])) == prime_powers(10^6)
     26        True
     27        sage: set(prime_powers_intlist(10^6)[1][1:]) == set(primes(10^6))
     28        True
     29    """
     30    v = prime_range(B, py_ints=True)
     31    cdef IntList w = IntList(len(v)*2), w0 = IntList(len(v)*2)
     32    w[0] = 1
     33    w0[0] = 1
     34    # Now fill in prime powers
     35    cdef Py_ssize_t i=1
     36    cdef int p
     37    cdef long long q   # so that the "q < B" test doesn't fail due to overflow.
     38    for p in v:
     39        q = p
     40        while q < B:
     41            w._values[i] = q
     42            w0._values[i] = p
     43            q *= p
     44            i += 1
     45    return w[:i], w0[:i]
     46
     47def extend_multiplicatively(IntList a):
     48    """
     49    Given an IntList a such that the a[p^r] is filled in, for all
     50    prime powers p^r, fill in all the other a[n] multiplicatively.
     51
     52    INPUT:
     53        - a -- IntList with prime-power-index entries set; all other
     54          entries are ignored
     55
     56    OUTPUT:
     57        - the input object a is modified to have all entries set
     58          via multiplicativity.
     59
     60    EXAMPLES::
     61
     62        sage: from sage.lfunctions.eulerprod_fast import extend_multiplicatively
     63        sage: B = 10^5; E = EllipticCurve('389a'); an = stats.IntList(B)
     64        sage: for pp in prime_powers(B):
     65        ...     an[pp] = E.an(pp)
     66        ...
     67        sage: extend_multiplicatively(an)
     68        sage: list(an) == E.anlist(len(an))[:len(an)]
     69        True
     70    """
     71    cdef Py_ssize_t i, B = len(a)
     72    cdef IntList P, P0
     73    P, P0 = prime_powers_intlist(B)
     74
     75    # Known[n] = 1 if a[n] is set.  We use this separate array
     76    # to avoid using a sentinel value in a.
     77    cdef IntList known = IntList(B)
     78    for i in range(len(P)):
     79        known._values[P[i]] = 1
     80       
     81    cdef int k, pp, p, n
     82    # fill in the multiples of pp = prime power
     83    for i in range(len(P)):
     84        pp = P._values[i]; p = P0._values[i]
     85        k = 2
     86        n = k*pp
     87        while n < B:
     88            # only consider n exactly divisible by pp
     89            if k % p and known._values[k]:
     90                a._values[n] = a._values[k] * a._values[pp]
     91                known._values[n] = 1
     92            n += pp
     93            k += 1
     94       
     95       
     96def extend_multiplicatively_generic(list a):
     97    """
     98    Given a list a of numbers such that the a[p^r] is filled in, for
     99    all prime powers p^r, fill in all the other a[n] multiplicatively.
     100
     101    INPUT:
     102        - a -- list with prime-power-index entries set; all other
     103          entries are ignored
     104
     105    OUTPUT:
     106        - the input object a is modified to have all entries set
     107          via multiplicativity.
     108
     109    EXAMPLES::
     110
     111        sage: from sage.lfunctions.eulerprod_fast import extend_multiplicatively_generic
     112        sage: B = 10^5; E = EllipticCurve('389a'); an = [0]*B
     113        sage: for pp in prime_powers(B):
     114        ...     an[pp] = E.an(pp)
     115        ...
     116        sage: extend_multiplicatively_generic(an)
     117        sage: list(an) == E.anlist(len(an))[:len(an)]
     118        True
     119
     120
     121    A test using large integers::
     122   
     123        sage: v = [0, 1, 2**100, 3**100, 4, 5, 0]
     124        sage: extend_multiplicatively_generic(v)
     125        sage: v
     126        [0, 1, 1267650600228229401496703205376, 515377520732011331036461129765621272702107522001, 4, 5, 653318623500070906096690267158057820537143710472954871543071966369497141477376]
     127        sage: v[-1] == v[2]*v[3]
     128        True
     129
     130    A test using variables::
     131
     132        sage: v = list(var(' '.join('x%s'%i for i in [0..30]))); v
     133        [x0, x1, x2, x3, x4, x5, x6, x7, x8, x9, x10, x11, x12, x13, x14, x15, x16, x17, x18, x19, x20, x21, x22, x23, x24, x25, x26, x27, x28, x29, x30]
     134        sage: extend_multiplicatively_generic(v)
     135        sage: v
     136        [x0, x1, x2, x3, x4, x5, x2*x3, x7, x8, x9, x2*x5, x11, x3*x4, x13, x2*x7, x3*x5, x16, x17, x2*x9, x19, x4*x5, x3*x7, x11*x2, x23, x3*x8, x25, x13*x2, x27, x4*x7, x29, x2*x3*x5]
     137    """
     138    cdef Py_ssize_t i, B = len(a)
     139    cdef IntList P, P0
     140    P, P0 = prime_powers_intlist(B)
     141
     142    # Known[n] = 1 if a[n] is set.  We use this separate array
     143    # to avoid using a sentinel value in a.
     144    cdef IntList known = IntList(B)
     145    for i in range(len(P)):
     146        known._values[P[i]] = 1
     147       
     148    cdef int k, pp, p, n
     149    # fill in the multiples of pp = prime power
     150    for i in range(len(P)):
     151        pp = P._values[i]; p = P0._values[i]
     152        k = 2
     153        n = k*pp
     154        while n < B:
     155            # only consider n exactly divisible by pp
     156            if k % p and known._values[k]:
     157                a[n] = a[k] * a[pp]
     158                known._values[n] = 1
     159            n += pp
     160            k += 1
     161   
     162   
     163