Ticket #11024: trac-11024-dokchitser.patch

File trac-11024-dokchitser.patch, 43.3 KB (added by Martin Raum, 12 years ago)
  • sage/lfunctions/dokchitser.py

    diff -r fbb0ca1d026e sage/lfunctions/dokchitser.py
    a b  
    88
    99- William Stein (2006-03-08): Sage interface
    1010
     11- Henri Cohen (2011-03): Convert GP script to C file
     12
     13- Martin Raum (2011-03): Refactor the wrapper and adapt to C file
     14
    1115TODO:
    1216
    1317- add more examples from data/extcode/pari/dokchitser that illustrate
     
    2832import copy
    2933
    3034from sage.structure.sage_object import SageObject
     35from sage.rings.all import ComplexField, RealField, Integer, PowerSeriesRing
     36from sage.misc.all import verbose, sage_eval
     37from sage.schemes.all import is_EllipticCurve
     38import sage.interfaces.gp
     39import os
    3140
    32 from sage.rings.all import ComplexField, RealField, Integer
    33 
    34 from sage.misc.all import verbose, sage_eval
    35 
    36 from sage.schemes.all import is_EllipticCurve
    37 
    38 import sage.interfaces.gp
    39 
    40 class Dokchitser(SageObject):
     41class Dokchitser ( SageObject ) :
    4142    r"""
    4243    Dokchitser's `L`-functions Calculator
    4344   
     
    4647    Dokchitser(conductor, gammaV, weight, eps, poles, residues, init,
    4748    prec)
    4849   
    49     where
     50    where:
    5051   
    51     - ``conductor`` - integer, the conductor
     52        - ``conductor`` - integer, the conductor
    5253
    53     - ``gammaV`` - list of Gamma-factor parameters, e.g. [0] for
    54       Riemann zeta, [0,1] for ell.curves, (see examples).
     54        - ``gammaV`` - list of Gamma-factor parameters, e.g. [0] for
     55          Riemann zeta, [0,1] for ell.curves, (see examples).
    5556
    56     - ``weight`` - positive real number, usually an integer e.g. 1 for
    57       Riemann zeta, 2 for `H^1` of curves/`\QQ`
     57        - ``weight`` - positive real number, usually an integer e.g. 1 for
     58          Riemann zeta, 2 for `H^1` of curves/`\QQ`
    5859
    59     - ``eps`` - complex number; sign in functional equation
     60        - ``eps`` - complex number; sign in functional equation
    6061
    61     - ``poles`` - (default: []) list of points where `L^*(s)` has
    62       (simple) poles; only poles with `Re(s)>weight/2` should be
    63       included
     62        - ``poles`` - (default: []) list of points where `L^*(s)` has
     63          (simple) poles; only poles with `Re(s)>weight/2` should be
     64          included
    6465
    65     - ``residues`` - vector of residues of `L^*(s)` in those poles or
    66       set residues='automatic' (default value)
     66        - ``residues`` - vector of residues of `L^*(s)` in those poles or
     67          set residues='automatic' (default value)
    6768
    68     - ``prec`` - integer (default: 53) number of *bits* of precision
    69    
    70     RIEMANN ZETA FUNCTION:
     69        - ``prec`` - integer (default: 53) number of *bits* of precision
     70
     71    EXAMPLES:
    7172   
    7273    We compute with the Riemann Zeta function.
    7374   
    7475    ::
    7576   
    7677        sage: L = Dokchitser(conductor=1, gammaV=[0], weight=1, eps=1, poles=[1], residues=[-1], init='1')
    77         sage: L
    78         Dokchitser L-series of conductor 1 and weight 1
    7978        sage: L(1)
     79        ...
    8080        Traceback (most recent call last):
    8181        ...
    82         ArithmeticError
     82        ArithmeticError: pole occured in calculations
    8383        sage: L(2)
    8484        1.64493406684823
    8585        sage: L(2, 1.1)
     
    9292        sage: L.taylor_series(2, k=5)
    9393        1.64493406684823 - 0.937548254315844*z + 0.994640117149451*z^2 - 1.00002430047384*z^3 + 1.00006193307...*z^4 + O(z^5)
    9494   
    95     RANK 1 ELLIPTIC CURVE:
    96    
    9795    We compute with the `L`-series of a rank `1`
    9896    curve.
    99    
     97
    10098    ::
    10199   
    102100        sage: E = EllipticCurve('37a')
     
    109107        sage: L.derivative(1,2)
    110108        0.373095594536324
    111109        sage: L.num_coeffs()
    112         48
     110        57
    113111        sage: L.taylor_series(1,4)
    114112        0.305999773834052*z + 0.186547797268162*z^2 - 0.136791463097188*z^3 + O(z^4)
    115113        sage: L.check_functional_equation()
    116114        6.11218974800000e-18                            # 32-bit
    117115        6.04442711160669e-18                            # 64-bit
    118116   
    119     RANK 2 ELLIPTIC CURVE:
    120    
    121117    We compute the leading coefficient and Taylor expansion of the
    122118    `L`-series of a rank `2` curve.
    123119   
     
    126122        sage: E = EllipticCurve('389a')
    127123        sage: L = E.lseries().dokchitser()
    128124        sage: L.num_coeffs ()
    129         156
     125        187
    130126        sage: L.derivative(1,E.rank())
    131127        1.51863300057685
    132128        sage: L.taylor_series(1,4)
    133129        2.90759778535572e-20 + (-1.64772676916085e-20)*z + 0.759316500288427*z^2 - 0.430302337583362*z^3 + O(z^4)      # 32-bit
    134130        -3.11623283109075e-21 + (1.76595961125962e-21)*z + 0.759316500288427*z^2 - 0.430302337583362*z^3 + O(z^4)      # 64-bit
    135131   
    136     RAMANUJAN DELTA L-FUNCTION:
    137    
    138132    The coefficients are given by Ramanujan's tau function::
    139133   
    140134        sage: L = Dokchitser(conductor=1, gammaV=[0,1], weight=12, eps=1)
     
    149143    ::
    150144   
    151145        sage: L.num_coeffs()
    152         12
     146        15
    153147        sage: L.set_coeff_growth('2*n^(11/2)')
    154148        sage: L.num_coeffs()
    155         11
     149        13
    156150   
    157151    Now we're ready to evaluate, etc.
    158152   
     
    165159        sage: L.taylor_series(1,3)
    166160        0.0374412812685155 + 0.0709221123619322*z + 0.0380744761270520*z^2 + O(z^3)
    167161    """
    168     def __init__(self, conductor, gammaV, weight, eps, \
    169                        poles=[], residues='automatic', prec=53,
    170                        init=None):
     162    def __init__(self, conductor, gammaV, weight, eps,
     163                       poles = [], residues = 'automatic', prec = 53,
     164                       init = None) :
    171165        """
    172         Initialization of Dokchitser calculator EXAMPLES::
     166        Initialization of Dokchitser calculator
     167
     168        EXAMPLES::
    173169       
    174170            sage: L = Dokchitser(conductor=1, gammaV=[0], weight=1, eps=1, poles=[1], residues=[-1], init='1')
    175171            sage: L.num_coeffs()
    176             4
     172            5
    177173        """
     174        self.__CC   = ComplexField(prec)
     175        self.__RR   = self.__CC._real_field()
     176
    178177        self.conductor = conductor
    179178        self.gammaV = gammaV
    180179        self.weight = weight
     
    182181        self.poles  = poles
    183182        self.residues = residues
    184183        self.prec   = prec
    185         self.__CC   = ComplexField(self.prec)
    186         self.__RR   = self.__CC._real_field()
     184
     185        self.__max_imaginary_part = 0
     186        self.__max_asymp_coeffs = 40
     187
     188        self._init_parameters()
     189
     190        self.__initialized = False
     191        self.__init_coeff_data = None
    187192        if not init is None:
    188193            self.init_coeffs(init)
    189             self.__init = init
    190         else:
    191             self.__init = False
    192194
    193     def __reduce__(self):
    194         D = copy.copy(self.__dict__)
    195         if D.has_key('_Dokchitser__gp'):
    196             del D['_Dokchitser__gp']
    197         return reduce_load_dokchitser, (D, )
    198        
    199     def _repr_(self):
    200         z = "Dokchitser L-series of conductor %s and weight %s"%(
    201                    self.conductor, self.weight)
    202         return z
     195        self.__coeff_growth = None
    203196
    204     def __del__(self):
     197    def __del__(self) :
     198        """
     199        TESTS::
     200
     201            sage: L = Dokchitser(conductor=1, gammaV=[0], weight=1, eps=1, poles=[1], residues=[-1], init='1')
     202            sage: del L
     203        """
    205204        self.gp().quit()
    206205
    207     def gp(self):
     206    def __reduce__(self) :
     207        """
     208        TESTS::
     209
     210            sage: L = Dokchitser(conductor=1, gammaV=[0], weight=1, eps=1, poles=[1], residues=[-1], init='1')
     211            sage: L = loads(dumps(L)) ## indirect doctest
     212            sage: L(2)
     213            1.64493406684823
     214        """
     215        D = copy.copy(self.__dict__)
     216        if '_Dokchitser__gp' in D :
     217            del D['_Dokchitser__gp']
     218
     219        return (reduce_load_dokchitser, (D, ))
     220       
     221    def _repr_(self) :
     222        """
     223        TESTS::
     224           
     225            sage: Dokchitser(conductor=1, gammaV=[0], weight=1, eps=1, poles=[1], residues=[-1], init='1')
     226            Dokchitser L-series of conductor 1 and weight 1
     227        """
     228        return "Dokchitser L-series of conductor %s and weight %s" \
     229                   % ( self.conductor, self.weight )
     230
     231    ###########################################################################
     232    ### Initialization of variables for the Dokchitser calculator
     233    ###########################################################################
     234
     235    def _init_parameters(self, max_imaginary_part = 0, max_asymp_coeffs = 40) :
     236        """
     237        INPUT:
     238
     239            -  ``max_imaginary_part`` - (default: 0): redefine if
     240               you want to compute L(s) for s having large imaginary part,
     241           
     242            -  ``max_asymp_coeffs`` - (default: 40): at most this
     243               many terms are generated in asymptotic series for phi(t) and
     244               G(s,t).
     245
     246        TESTS::
     247       
     248            sage: L = Dokchitser(conductor=1, gammaV=[0], weight=1, eps=1, poles=[1], residues=[-1], init='1') ## indirect doctest
     249            sage: L._init_parameters( max_imaginary_part = 2, max_asymp_coeffs = 100)
     250        """
     251        if self.__max_imaginary_part != max_imaginary_part or \
     252           self.__max_asymp_coeffs != max_asymp_coeffs :
     253            self.__max_imaginary_part = max_imaginary_part
     254            self.__max_asymp_coeffs = max_asymp_coeffs
     255           
     256        self._gp_eval('default(realprecision, %s)'%(self.prec//3 + 2))
     257        self._gp_eval('conductor = %s'%self.conductor)
     258        self._gp_eval('gammaV = %s'%self.gammaV)
     259        self._gp_eval('weight = %s'%self.weight)
     260        self._gp_eval('sgn = %s'%self.eps)
     261        self._gp_eval('Lpoles = %s'%self.poles)
     262        self._gp_eval('Lresidues = %s'%self.residues)
     263
     264        self._gp_eval('MaxImaginaryPart = %s' % self.__max_imaginary_part)
     265        self._gp_eval('MaxAsymptoticCoeffs = %s' % self.__max_asymp_coeffs)
     266       
     267        self._gp_eval('initall([conductor,gammaV,weight,sgn, Lpoles, Lresidues, MaxImaginaryPart, MaxAsymptoticCoeffs])')
     268
     269        self.__initialized = False
     270
     271    def set_coeff_growth(self, coefgrow) :
     272        r"""
     273        You might have to redefine the coefficient growth function if the
     274        `a_n` of the `L`-series are not given by the
     275        following PARI function::
     276       
     277            coefgrow(n) = if(length(Lpoles),   
     278                              1.5*n^(vecmax(real(Lpoles))-1), 
     279                              sqrt(4*n)^(weight-1));
     280                   
     281       
     282        INPUT:
     283       
     284            -  ``coefgrow`` - string that evaluates to a PARI
     285               function of n that defines a coefgrow function.
     286       
     287        EXAMPLE::
     288       
     289            sage: L = Dokchitser(conductor=1, gammaV=[0,1], weight=12, eps=1)
     290            sage: pari_precode = 'tau(n)=(5*sigma(n,3)+7*sigma(n,5))*n/12 - 35*sum(k=1,n-1,(6*k-4*(n-k))*sigma(k,3)*sigma(n-k,5))'
     291            sage: L.init_coeffs('tau(k)', pari_precode=pari_precode)
     292            sage: L.set_coeff_growth('2*n^(11/2)')
     293            sage: L(1)
     294            0.0374412812685155
     295        """
     296        if coefgrow is None :
     297            return
     298        if not isinstance(coefgrow, str) :
     299            raise TypeError( "coefgrow must be a string" )
     300
     301        self.__coeff_growth = coefgrow
     302
     303        self._gp_eval('grow_closure = (n -> (%s))' % coefgrow.replace('\n',' '))
     304        self._gp_eval('setmycoefgrow(grow_closure)')
     305
     306    def init_coeffs(self, v, cutoff = 1,
     307                             w = None,
     308                             pari_precode = '',
     309                             max_imaginary_part = None,
     310                             max_asymp_coeffs = None) :
     311        """
     312        Set the coefficients `a_n` of the `L`-series. If
     313        `L(s)` is not equal to its dual, pass the coefficients of
     314        the dual as the second optional argument.
     315        If `v` is ``None`` and the `L`-series has already been
     316        initialized, you may reset the maximal imaginary part or
     317        the maximal number of asymptotic coefficients.
     318       
     319        INPUT:
     320       
     321            -  `v` - list of complex numbers or string (pari
     322               function of k) or None
     323           
     324            -  ``cutoff`` - real number = 1 (default: 1)
     325           
     326            -  `w` - list of complex numbers or string (pari
     327               function of k)
     328           
     329            -  ``pari_precode`` - some code to execute in pari
     330               before calling initLdata
     331
     332            -  ``max_imaginary_part`` - (default: ``None``): redefine if
     333               you want to compute L(s) for s having large imaginary part,
     334           
     335            -  ``max_asymp_coeffs`` - (default: ``None``): at most this
     336               many terms are generated in asymptotic series for phi(t) and
     337               G(s,t).
     338       
     339        EXAMPLES::
     340       
     341            sage: L = Dokchitser(conductor=1, gammaV=[0,1], weight=12, eps=1)
     342            sage: pari_precode = 'tau(n)=(5*sigma(n,3)+7*sigma(n,5))*n/12 - 35*sum(k=1,n-1,(6*k-4*(n-k))*sigma(k,3)*sigma(n-k,5))'
     343            sage: L.init_coeffs('tau(k)', pari_precode=pari_precode)
     344
     345        Evaluate the resulting L-function at a point, and compare with the answer that
     346        one gets "by definition" (of L-function attached to a modular form)::
     347
     348        ::
     349
     350            sage: L(14)
     351            0.998583063162746
     352            sage: a = delta_qexp(1000)
     353            sage: sum(a[n]/float(n)^14 for n in range(1,1000))
     354            0.99858306316274592
     355
     356        Illustrate that one can give a list of complex numbers for v (see trac 10937)::
     357
     358            sage: L2 = Dokchitser(conductor=1, gammaV=[0,1], weight=12, eps=1)
     359            sage: L2.init_coeffs(list(delta_qexp(1000))[1:])
     360            sage: L2(14)
     361            0.998583063162746
     362
     363        TESTS:
     364
     365            Verify that setting the w parameter doesn't raise an error
     366            (see trac 10937).  Note that the meaning of w does not seem to
     367            be documented anywhere in Dokchitser's package yet, so there is
     368            no claim that the example below is meaningful!
     369
     370        ::
     371
     372            sage: L2 = Dokchitser(conductor=1, gammaV=[0,1], weight=12, eps=1)
     373            sage: L2.init_coeffs(list(delta_qexp(1000))[1:], w=[1..1000])
     374
     375        ::
     376           
     377            sage: L2 = Dokchitser(conductor=1, gammaV=[0,1], weight=12, eps=1)
     378            sage: L2.init_coeffs(list(delta_qexp(1000))[1:])
     379            sage: L2.init_coeffs(None, max_imaginary_part = 2, max_asymp_coeffs = 50)
     380            sage: L2(14)
     381            0.998583063162746
     382        """
     383        if v is None and self.__initialized :
     384            (v, cutoff, w, pari_precode) = self.__init_coeff_data
     385        else :
     386            self.__init_coeff_data = (v, cutoff, w, pari_precode)
     387            self.__initialized = True
     388
     389        params = dict()
     390        if not max_imaginary_part is None :
     391            params['max_imaginary_part'] = max_imaginary_part
     392        if not max_asymp_coeffs is None :
     393            params['max_asymp_coeffs'] = max_asymp_coeffs
     394        if len(params) != 0 :
     395            self._init_parameters(**params)
     396            self.__initialized = True
     397       
     398
     399        CC = self.__CC
     400        cutoff = self.__RR(cutoff)
     401
     402        # Execute the precode
     403        if pari_precode != '':
     404            self._gp_eval(pari_precode)
     405
     406        # Pari code is passed as the coefficient function
     407        if isinstance(v, str):
     408            # to maintain compatibility with the Dokchitser scripts we convert
     409            # the expression into a closure
     410            v = "(k -> (%s))" % v
     411            self._gp_eval('vcoeff_closure = %s;' % v)
     412
     413            # if the dual coefficients are not passed
     414            if w is None :
     415                self._gp_eval('initLdata(vcoeff_closure, %s);' % cutoff)
     416                return
     417
     418            w = "(k -> (%s))" % w
     419            self._gp_eval('wcoeff_closure = %s;' % w)
     420
     421            self._gp_eval('initLdata(vcoeff_closure, %s, wcoeff_closure);' % cutoff)
     422            return
     423
     424        # check that v is a list or a tuple of coefficients
     425        if not isinstance(v, (list, tuple)):
     426            raise TypeError( "v (=%s) must be a list, tuple, or string" % v )
     427
     428
     429        self._gp_eval('Avec = [%s];' % ','.join([CC(a)._pari_init_() for a in v]))
     430        self._gp_eval('Avec_closure = (k -> Avec[k]);')
     431
     432        # if the dual coefficients are not passed
     433        if w is None:
     434            self._gp_eval('initLdata(Avec_closure, %s);' % cutoff)
     435            return
     436
     437        self._gp_eval('Bvec = [%s];' % ','.join([CC(a)._pari_init_() for a in w]))
     438        self._gp_eval('Bvec_closure = (k -> Bvec[k]);')
     439
     440        self._gp_eval('initLdata(Avec_closure, %s, Bvec_closure);' % cutoff)
     441
     442    def __check_init(self, s = None) :
     443        """
     444        Check wheter the class has been initialized to compute L-values with
     445        imaginary parts like `s` has.
     446
     447        TESTS::
     448           
     449            sage: L = Dokchitser(conductor=1, gammaV=[0], weight=1, eps=1, poles=[1], residues=[-1])
     450            sage: L(1) ## indirect doctest
     451            Traceback (most recent call last):
     452            ...
     453            ValueError: you must call init_coeffs on the L-function first
     454            sage: L.init_coeffs('1')
     455            sage: L(1 + I) ## indirect doctest
     456            Traceback (most recent call last):
     457            ...
     458            ValueError: maximal imaginary part was set to 0 exceeding abs(Im(s)) = 1.00000000000000; use init_coeffs(None, max_imaginary_part = 1.00000000000000)
     459           """
     460        if not self.__initialized :
     461            raise ValueError( "you must call init_coeffs on the L-function first" )
     462
     463        if not s is None and self.__max_imaginary_part < abs(s.imag()) :
     464            raise ValueError( (   "maximal imaginary part was set to %s exceeding abs(Im(s)) = %s; "
     465                                + "use init_coeffs(None, max_imaginary_part = %s)" ) \
     466                              % (self.__max_imaginary_part, s.imag(), abs(s.imag())) )
     467
     468    ###########################################################################
     469    ### PARI/GP service methods
     470    ###########################################################################
     471
     472    def gp(self) :
    208473        """
    209474        Return the gp interpreter that is used to implement this Dokchitser
    210475        L-function.
     
    222487            return self.__gp
    223488        except AttributeError:
    224489            g = sage.interfaces.gp.Gp(script_subdirectory='dokchitser',
    225                                       logfile=None) 
    226             g.read('computel.gp')
     490                                      logfile=None)
    227491            self.__gp = g
    228             self._gp_eval('default(realprecision, %s)'%(self.prec//3 + 2))
    229             self._gp_eval('conductor = %s'%self.conductor)
    230             self._gp_eval('gammaV = %s'%self.gammaV)
    231             self._gp_eval('weight = %s'%self.weight)
    232             self._gp_eval('sgn = %s'%self.eps)
    233             self._gp_eval('Lpoles = %s'%self.poles)
    234             self._gp_eval('Lresidues = %s'%self.residues)
    235             g._dokchitser = True
     492
     493            libraries = ["computel5.gp." + suffix for suffix in ["so", "dll", "dylib", "dylib64"]]
     494            while True :
     495                library = libraries.pop()
     496                try :
     497                    g('GP;install("init_computel5","v","init_computel5","./%s");' % library)
     498                except :
     499                    continue
     500                break
     501           
     502            map(g, ['GP;install("initglobs","","initglobs","./%s");' % library,
     503                    'GP;install("getconductor","","getconductor","./%s");' % library,
     504                    'GP;install("getgammaV","","getgammaV","./%s");' % library,
     505                    'GP;install("getweight","","getweight","./%s");' % library,
     506                    'GP;install("getsgn","","getsgn","./%s");' % library,
     507                    'GP;install("getLpoles","","getLpoles","./%s");' % library,
     508                    'GP;install("getLresidues","","getLresidues","./%s");' % library,
     509                    'GP;install("getMaxImaginaryPart","","getMaxImaginaryPart","./%s");' % library,
     510                    'GP;install("getMaxAsympCoeffs","","getMaxAsympCoeffs","./%s");' % library,
     511                    'GP;install("setglobs","D0,G,","setglobs","./%s");' % library,
     512                    'GP;install("initall","D0,G,","initall","./%s");' % library,
     513                    'GP;install("setconductor","D0,G,","setconductor","./%s");' % library,
     514                    'GP;install("setgammaV","D0,G,","setgammaV","./%s");' % library,
     515                    'GP;install("setweight","D0,G,","setweight","./%s");' % library,
     516                    'GP;install("setsgn","D0,G,","setsgn","./%s");' % library,
     517                    'GP;install("setLpoles","D0,G,","setLpoles","./%s");' % library,
     518                    'GP;install("setLresidues","D0,G,","setLresidues","./%s");' % library,
     519                    'GP;install("setMaxImaginaryPart","D0,G,","setMaxImaginaryPart","./%s");' % library,
     520                    'GP;install("setMaxAsympCoeffs","D0,G,","setMaxAsympCoeffs","./%s");' % library,
     521                    'GP;install("setmycoefgrow","D0,G,","setmycoefgrow","./%s");' % library,
     522                    'GP;install("coefgrow","D0,G,p","coefgrow","./%s");' % library,
     523                    'GP;install("vecleft","D0,G,D0,G,","vecleft","./%s");' % library,
     524                    'GP;install("vecright","D0,G,D0,G,","vecright","./%s");' % library,
     525                    'GP;install("StrTab","D0,G,D0,G,","StrTab","./%s");' % library,
     526                    'GP;install("concatstr","DGDGDGDG","concatstr","./%s");' % library,
     527                    'GP;install("errprint","D0,G,p","errprint","./%s");' % library,
     528                    'GP;install("gammaseries","D0,G,D0,G,p","gammaseries","./%s");' % library,
     529                    'GP;install("fullgamma","D0,G,p","fullgamma","./%s");' % library,
     530                    'GP;install("fullgammaseries","D0,G,D0,G,p","fullgammaseries","./%s");' % library,
     531                    'GP;install("RecursionsAtInfinity","p","RecursionsAtInfinity","./%s");' % library,
     532                    'GP;install("SeriesToContFrac","D0,G,p","SeriesToContFrac","./%s");' % library,
     533                    'GP;install("EvaluateContFrac","D0,G,D0,G,D0,G,p","EvaluateContFrac","./%s");' % library,
     534                    'GP;install("cflength","DGp","cflength","./%s");' % library,
     535                    'GP;install("initLdata","lD0,G,DGDGp","initLdata","./%s");' % library,
     536                    'GP;install("phi","D0,G,p","phi","./%s");' % library,
     537                    'GP;install("RecursephiV","v","RecursephiV","./%s");' % library,
     538                    'GP;install("phi0","D0,G,p","phi0","./%s");' % library,
     539                    'GP;install("phiinf","D0,G,DGp","phiinf","./%s");' % library,
     540                    'GP;install("G","D0,G,D0,G,DGp","G","./%s");' % library,
     541                    'GP;install("LogInt","D0,G,D0,G,D0,G,D0,G,p","LogInt","./%s");' % library,
     542                    'GP;install("MakeLogSum","D0,G,D0,G,p","MakeLogSum","./%s");' % library,
     543                    'GP;install("G0","D0,G,D0,G,D0,G,p","G0","./%s");' % library,
     544                    'GP;install("Ginf","D0,G,D0,G,D0,G,DGp","Ginf","./%s");' % library,
     545                    'GP;install("initGinf","lD0,G,D0,G,p","initGinf","./%s");' % library,
     546                    'GP;install("ltheta","D0,G,p","ltheta","./%s");' % library,
     547                    'GP;install("ldualtheta","D0,G,p","ldualtheta","./%s");' % library,
     548                    'GP;install("checkfeq","DGp","checkfeq","./%s");' % library,
     549                    'GP;install("L","D0,G,DGDGp","L","./%s");' % library,
     550                    'GP;install("Lseries","D0,G,DGDGp","Lseries","./%s");' % library,
     551                    'GP;install("Lstar","D0,G,DGDGp","Lstar","./%s");' % library] )
     552
    236553            return g
    237554
    238     def _gp_eval(self, s):
     555    def _gp_eval(self, s) :
     556        """
     557        TESTS::
     558           
     559            sage: L = Dokchitser(conductor=1, gammaV=[0], weight=1, eps=1, poles=[1], residues=[-1])
     560            sage: L._gp_eval('a := 2;')
     561            Traceback (most recent call last):
     562            ...
     563            RuntimeError: PARI Error.
     564              ***   syntax error, unexpected '=', expecting variable name: a:=2;
     565              ***                                                            ^---
     566        """
    239567        try:
    240568            t = self.gp().eval(s)
    241         except (RuntimeError, TypeError):
    242             raise RuntimeError, "Unable to create L-series, due to precision or other limits in PARI."
    243         if '***' in t:
    244             raise RuntimeError, "Unable to create L-series, due to precision or other limits in PARI."
     569        except (RuntimeError, TypeError) as e:
     570            raise RuntimeError( "PARI Error.\n" + str(e))
     571        if '***' in t and \
     572           not 'division by zero' in t :
     573            raise RuntimeError( "PARI Error.\n" + str(t))
     574
    245575        return t
    246                
    247     def __check_init(self):
    248         if not self.__init:
    249             raise ValueError, "you must call init_coeffs on the L-function first"
    250576
    251     def num_coeffs(self, T=1):
     577    ###########################################################################
     578    ### information about the L-series
     579    ###########################################################################
     580
     581    def number_of_needed_coeffs(self, T = 1.2) :
    252582        """
    253583        Return number of coefficients `a_n` that are needed in
    254584        order to perform most relevant `L`-function computations to
     
    259589            sage: E = EllipticCurve('11a')
    260590            sage: L = E.lseries().dokchitser()
    261591            sage: L.num_coeffs()
    262             26
     592            31
    263593            sage: E = EllipticCurve('5077a')
    264594            sage: L = E.lseries().dokchitser()
    265595            sage: L.num_coeffs()
    266             568
     596            683
    267597            sage: L = Dokchitser(conductor=1, gammaV=[0], weight=1, eps=1, poles=[1], residues=[-1], init='1')
    268598            sage: L.num_coeffs()
    269             4
     599            5
    270600        """
    271         return Integer(self.gp().eval('cflength(%s)'%T))
     601        return Integer( self._gp_eval('cflength(%s)' % self.__RR(T)) )
    272602
    273     def init_coeffs(self, v, cutoff=1,
    274                              w=None,
    275                              pari_precode='',
    276                              max_imaginary_part=0,
    277                              max_asymp_coeffs=40):
    278         """
    279         Set the coefficients `a_n` of the `L`-series. If
    280         `L(s)` is not equal to its dual, pass the coefficients of
    281         the dual as the second optional argument.
    282        
    283         INPUT:
    284        
    285        
    286         -  ``v`` - list of complex numbers or string (pari
    287            function of k)
    288        
    289         -  ``cutoff`` - real number = 1 (default: 1)
    290        
    291         -  ``w`` - list of complex numbers or string (pari
    292            function of k)
    293        
    294         -  ``pari_precode`` - some code to execute in pari
    295            before calling initLdata
    296        
    297         -  ``max_imaginary_part`` - (default: 0): redefine if
    298            you want to compute L(s) for s having large imaginary part,
    299        
    300         -  ``max_asymp_coeffs`` - (default: 40): at most this
    301            many terms are generated in asymptotic series for phi(t) and
    302            G(s,t).
    303        
    304        
    305         EXAMPLES::
    306        
    307             sage: L = Dokchitser(conductor=1, gammaV=[0,1], weight=12, eps=1)
    308             sage: pari_precode = 'tau(n)=(5*sigma(n,3)+7*sigma(n,5))*n/12 - 35*sum(k=1,n-1,(6*k-4*(n-k))*sigma(k,3)*sigma(n-k,5))'
    309             sage: L.init_coeffs('tau(k)', pari_precode=pari_precode)
    310         """
    311         if isinstance(v, tuple) and w is None:
    312             v, cutoff, w, pari_precode, max_imaginary_part, max_asymp_coeffs = v
     603    ## We maintain the old name for compatibility, but unify the name
     604    ## style throughout the class
     605    num_coeffs = number_of_needed_coeffs
    313606
    314         self.__init = (v, cutoff, w, pari_precode, max_imaginary_part, max_asymp_coeffs)
    315         gp = self.gp()
    316         if pari_precode != '':
    317             self._gp_eval(pari_precode)
    318         RR = self.__CC._real_field()
    319         cutoff = RR(cutoff)
    320         if isinstance(v, str):
    321             if w is None:
    322                 self._gp_eval('initLdata("%s", %s)'%(v, cutoff))
    323                 return
    324             self._gp_eval('initLdata("%s",%s,"%s")'%(v,cutoff,w))
    325             return
    326         if not isinstance(v, (list, tuple)):
    327             raise TypeError, "v (=%s) must be a list, tuple, or string"%v
    328         CC = self.__CC
    329         v = [CC(a)._pari_init_() for a in v]
    330         self._gp_eval('Avec = %s'%v)
    331         if w is None:
    332             self._gp_eval('initLdata("Avec[k]", %s)'%cutoff)
    333             return
    334         w = [CC(a)._pari_init_() for a in w]
    335         self._gp_eval('Bvec = %s'%w)
    336         self._gp_eval('initLdata("Avec[k]"),%s,"Bvec[k]"'%cutoff)
    337 
    338     def __to_CC(self, s):
    339         s = s.replace('.E','.0E').replace(' ','')
    340         return self.__CC(sage_eval(s, locals={'I':self.__CC.gen(0)}))
    341 
    342     def _clear_value_cache(self):
    343         del self.__values
    344 
    345     def __call__(self, s, c=None):
    346         r"""
    347         INPUT:
    348        
    349         -  ``s`` - complex number
    350        
    351         .. note::
    352 
    353            Evaluation of the function takes a long time, so each
    354            evaluation is cached. Call ``self._clear_value_cache()`` to
    355            clear the evaluation cache.
    356        
    357         EXAMPLES::
    358        
    359             sage: E = EllipticCurve('5077a')
    360             sage: L = E.lseries().dokchitser(100)
    361             sage: L(1)
    362             0
    363             sage: L(1+I)
    364             -1.3085436607849493358323930438 + 0.81298000036784359634835412129*I
    365         """
    366         self.__check_init()
    367         s = self.__CC(s)
    368         try:
    369             return self.__values[s]
    370         except AttributeError:
    371             self.__values = {}
    372         except KeyError:
    373             pass
    374         z = self.gp().eval('L(%s)'%s)
    375         if 'pole' in z:
    376             print z
    377             raise ArithmeticError
    378         elif '***' in z:
    379             print z
    380             raise RuntimeError
    381         elif 'Warning' in z:
    382             i = z.rfind('\n')
    383             msg = z[:i].replace('digits','decimal digits')
    384             verbose(msg, level=-1)
    385             ans = self.__to_CC(z[i+1:])
    386             self.__values[s] = ans
    387             return ans
    388         ans = self.__to_CC(z)
    389         self.__values[s] = ans
    390         return ans
    391        
    392     def derivative(self, s, k=1):
    393         r"""
    394         Return the `k`-th derivative of the `L`-series at
    395         `s`.
    396        
    397         .. warning::
    398 
    399            If `k` is greater than the order of vanishing of
    400            `L` at `s` you may get nonsense.
    401        
    402         EXAMPLES::
    403        
    404             sage: E = EllipticCurve('389a')
    405             sage: L = E.lseries().dokchitser()
    406             sage: L.derivative(1,E.rank())
    407             1.51863300057685
    408         """
    409         self.__check_init()       
    410         s = self.__CC(s)
    411         k = Integer(k)
    412         z = self.gp().eval('L(%s,,%s)'%(s,k))
    413         if 'pole' in z:
    414             raise ArithmeticError, z
    415         elif 'Warning' in z:
    416             i = z.rfind('\n')
    417             msg = z[:i].replace('digits','decimal digits')
    418             verbose(msg, level=-1)
    419             return self.__CC(z[i:])
    420         return self.__CC(z)
    421 
    422 
    423     def taylor_series(self, a=0, k=6, var='z'):
    424         r"""
    425         Return the first `k` terms of the Taylor series expansion
    426         of the `L`-series about `a`.
    427        
    428         This is returned as a series in ``var``, where you
    429         should view ``var`` as equal to `s-a`. Thus
    430         this function returns the formal power series whose coefficients
    431         are `L^{(n)}(a)/n!`.
    432        
    433         INPUT:
    434        
    435        
    436         -  ``a`` - complex number (default: 0); point about
    437            which to expand
    438        
    439         -  ``k`` - integer (default: 6), series is
    440            `O(``var``^k)`
    441        
    442         -  ``var`` - string (default: 'z'), variable of power
    443            series
    444        
    445        
    446         EXAMPLES::
    447        
    448             sage: L = Dokchitser(conductor=1, gammaV=[0], weight=1, eps=1, poles=[1], residues=[-1], init='1')
    449             sage: L.taylor_series(2, 3)
    450             1.64493406684823 - 0.937548254315844*z + 0.994640117149451*z^2 + O(z^3)
    451             sage: E = EllipticCurve('37a')
    452             sage: L = E.lseries().dokchitser()
    453             sage: L.taylor_series(1)
    454             0.305999773834052*z + 0.186547797268162*z^2 - 0.136791463097188*z^3 + 0.0161066468496401*z^4 + 0.0185955175398802*z^5 + O(z^6)
    455        
    456         We compute a Taylor series where each coefficient is to high
    457         precision.
    458        
    459         ::
    460        
    461             sage: E = EllipticCurve('389a')
    462             sage: L = E.lseries().dokchitser(200)
    463             sage: L.taylor_series(1,3)
    464             6.2240188634103774348273446965620801288836328651973234573133e-73 + (-3.527132447498646306292650465494647003849868770...e-73)*z + 0.75931650028842677023019260789472201907809751649492435158581*z^2 + O(z^3)
    465         """
    466         self.__check_init()       
    467         a = self.__CC(a)
    468         k = Integer(k)
    469         try:
    470             z = self.gp()('Vec(Lseries(%s,,%s))'%(a,k-1))
    471         except TypeError, msg:
    472             raise RuntimeError, "%s\nUnable to compute Taylor expansion (try lowering the number of terms)"%msg
    473         r = repr(z)
    474         if 'pole' in r:
    475             raise ArithmeticError, r
    476         elif 'Warning' in r:
    477             i = r.rfind('\n')
    478             msg = r[:i].replace('digits','decimal digits')
    479             verbose(msg, level=-1)
    480         v = list(z)
    481         K = self.__CC
    482         v = [K(repr(x)) for x in v]
    483         R = self.__CC[[var]]
    484         return R(v,len(v))
    485        
    486     def check_functional_equation(self, T=1.2):
     607    def check_functional_equation(self, T = 1.2) :
    487608        r"""
    488609        Verifies how well numerically the functional equation is satisfied,
    489         and also determines the residues if ``self.poles !=
    490         []`` and residues='automatic'.
     610        and also determines the residues if ``self.poles != []``
     611        and residues='automatic'.
    491612       
    492613        More specifically: for `T>1` (default 1.2),
    493614        ``self.check_functional_equation(T)`` should ideally
     
    530651            -9.73967861488124
    531652        """
    532653        self.__check_init()
    533         z = self.gp().eval('checkfeq(%s)'%T).replace(' ','')
     654
     655        z = self._gp_eval('checkfeq(%s)' % self.__RR(T)).replace(' ','')
     656
    534657        return self.__CC(z)
     658
     659
     660    ###########################################################################
     661    ### Evaluation
     662    ###########################################################################
     663
     664    def __call__(self, s, c = None) :
     665        r"""
     666        INPUT:
    535667       
    536     def set_coeff_growth(self, coefgrow):
     668            -  ``s`` - complex number
     669       
     670        .. NOTE::
     671
     672           Evaluation of the function takes a long time, so each
     673           evaluation is cached. Call ``self._clear_value_cache()`` to
     674           clear the evaluation cache.
     675       
     676        EXAMPLES::
     677       
     678            sage: E = EllipticCurve('5077a')
     679            sage: L = E.lseries().dokchitser(100)
     680            sage: L(1)
     681            0
     682            sage: L.init_coeffs(None, max_imaginary_part = 1)
     683            sage: L(1+I)
     684            -1.3085436607849493358323930438 + 0.81298000036784359634835412129*I
     685        """
     686        s = self.__CC(s)
     687        self.__check_init(s)
     688
     689        # first check whether we have cached the result
     690        try:
     691            return self.__values[s]
     692        except AttributeError:
     693            self.__values = {}
     694        except KeyError:
     695            pass
     696   
     697       
     698        z = self._gp_eval('L(%s)' % s)
     699        ans = self.__to_CC( self._parse_Lresult(z) )
     700        self.__values[s] = ans
     701 
     702        return ans
     703       
     704    def derivative(self, s, k = 1) :
    537705        r"""
    538         You might have to redefine the coefficient growth function if the
    539         `a_n` of the `L`-series are not given by the
    540         following PARI function::
     706        Return the `k`-th derivative of the `L`-series at
     707        `s`.
    541708       
    542                         coefgrow(n) = if(length(Lpoles),   
    543                                           1.5*n^(vecmax(real(Lpoles))-1), 
    544                                           sqrt(4*n)^(weight-1));
    545                    
     709        WARNING:
     710
     711           If `k` is greater than the order of vanishing of
     712           `L` at `s` you may get nonsense.
     713       
     714        EXAMPLES::
     715       
     716            sage: E = EllipticCurve('389a')
     717            sage: L = E.lseries().dokchitser()
     718            sage: L.derivative(1,E.rank())
     719            1.51863300057685
     720        """
     721        s = self.__CC(s)
     722        self.__check_init(s)
     723
     724        k = Integer(k)
     725
     726        z = self._gp_eval('L(%s,,%s)' % (s, k))
     727
     728        return self.__to_CC( self._parse_Lresult(z) )
     729
     730    def taylor_series(self, a = 0, k = 6, var = 'z') :
     731        r"""
     732        Return the first `k` terms of the Taylor series expansion
     733        of the `L`-series about `a`.
     734       
     735        This is returned as a series in ``var``, where you
     736        should view ``var`` as equal to `s-a`. Thus
     737        this function returns the formal power series whose coefficients
     738        are `L^{(n)}(a)/n!`.
    546739       
    547740        INPUT:
    548741       
    549742       
    550         -  ``coefgrow`` - string that evaluates to a PARI
    551            function of n that defines a coefgrow function.
     743        -  ``a`` - complex number (default: 0); point about
     744           which to expand
    552745       
     746        -  ``k`` - integer (default: 6), series is
     747           `O(``var``^k)`
    553748       
    554         EXAMPLE::
     749        -  ``var`` - string (default: 'z'), variable of power
     750           series
    555751       
    556             sage: L = Dokchitser(conductor=1, gammaV=[0,1], weight=12, eps=1)
    557             sage: pari_precode = 'tau(n)=(5*sigma(n,3)+7*sigma(n,5))*n/12 - 35*sum(k=1,n-1,(6*k-4*(n-k))*sigma(k,3)*sigma(n-k,5))'
    558             sage: L.init_coeffs('tau(k)', pari_precode=pari_precode)
    559             sage: L.set_coeff_growth('2*n^(11/2)')
    560             sage: L(1)
    561             0.0374412812685155
     752       
     753        EXAMPLES::
     754       
     755            sage: L = Dokchitser(conductor=1, gammaV=[0], weight=1, eps=1, poles=[1], residues=[-1], init='1')
     756            sage: L.taylor_series(2, 3)
     757            1.64493406684823 - 0.937548254315844*z + 0.994640117149451*z^2 + O(z^3)
     758            sage: E = EllipticCurve('37a')
     759            sage: L = E.lseries().dokchitser(200)
     760            sage: L.taylor_series(1)
     761            0.30599977383405230182048368332167647445263777459077199853454*z + 0.18654779726816196417381736877950759145408244520015683918099*z^2 - 0.13679146309718766630258221642815857076919470783750297065352*z^3 + 0.016106646849640053505097729456527078518195741393406783301132*z^4 + 0.018595517539880219615300779471612660822550180175154932080477*z^5 + O(z^6)
     762       
     763        We compute a Taylor series where each coefficient is to high
     764        precision.
     765       
     766        ::
     767       
     768            sage: E = EllipticCurve('389a')
     769            sage: L = E.lseries().dokchitser(100)
     770            sage: L.taylor_series(1,3)
     771            5.5808206662082541971182689498e-39 + (-3.1626339970101520275831845164e-39)*z + 0.75931650028842677023019260789*z^2 + O(z^3)
    562772        """
    563         if not isinstance(coefgrow, str):
    564             raise TypeError, "coefgrow must be a string"
    565         g = self.gp()
    566         g.eval('coefgrow(n) = %s'%(coefgrow.replace('\n',' ')))
     773        a = self.__CC(a)
     774        self.__check_init(a)
     775
     776        k = Integer(k)
     777
     778        try:
     779            z = self._gp_eval('Vec(Lseries(%s,,%s))' % (a, k-1))
     780        except TypeError as msg:
     781            raise RuntimeError( "%s\nUnable to compute Taylor expansion (try lowering the number of terms)" % msg )
     782
     783        # We need to convert back and forth in case there is a warning, that has
     784        # to be parsed before proceding
     785        r = self._parse_Lresult(repr(z))
     786        v = z[1:-1].split(',')
     787
     788        CC = self.__CC
     789        v = map(lambda x: CC(repr(x)), v)
     790
     791        R = PowerSeriesRing(CC, var)
     792
     793        return R(v, len(v))
     794
     795    @staticmethod
     796    def _parse_Lresult( z ) :
     797        """
     798        TESTS::
     799            sage: L = Dokchitser(conductor=1, gammaV=[0], weight=1, eps=1, poles=[1], residues=[-1], init='1')
     800            sage: L(1) ## indirect doctest
     801            ...
     802            Traceback (most recent call last):
     803            ...
     804            ArithmeticError: pole occured in calculations
     805        """
     806        if 'pole' in z or \
     807           'division by zero' in z :
     808            print z
     809            raise ArithmeticError( "pole occured in calculations" )
     810        elif 'Warning' in z:
     811            i = z.rfind('\n')
     812            msg = z[:i].replace('digits','decimal digits')
     813            verbose(msg, level=-1)
     814
     815            return z[i+1:]
     816
     817        return z
     818
     819    def __to_CC(self, s) :
     820        """
     821            TESTS::
     822                sage: L = Dokchitser(conductor=1, gammaV=[0], weight=1, eps=1, poles=[1], residues=[-1], init='1')
     823                sage: L(1.000000001) ## indirect doctest
     824                1.00000000056147e9
     825        """
     826        s = s.replace('.E','.0E').replace(' ','')
     827
     828        return self.__CC(sage_eval(s, locals={'I':self.__CC.gen(0)}))
     829
     830    def _clear_value_cache(self) :
     831        """
     832        TESTS::
     833                sage: L = Dokchitser(conductor=1, gammaV=[0], weight=1, eps=1, poles=[1], residues=[-1], init='1')
     834                sage: L(2)
     835                1.64493406684823
     836                sage: L._clear_value_cache()
     837                sage: L(2)
     838                1.64493406684823
     839        """
     840        del self.__values       
     841               
     842
     843def reduce_load_dokchitser(D) :
     844    """
     845    TESTS::
    567846       
    568        
    569 def reduce_load_dokchitser(D):
     847        sage: L = Dokchitser(conductor=1, gammaV=[0], weight=1, eps=1, poles=[1], residues=[-1], init='1')
     848        sage: L = loads(dumps(L)) ## indirect doctest
     849        sage: L(2)
     850        1.64493406684823
     851    """
    570852    X = Dokchitser(1,1,1,1)
    571853    X.__dict__ = D
    572     X.init_coeffs(X._Dokchitser__init)
     854    X._init_parameters()
     855
     856    if True :
     857        X.set_coeff_growth(X._Dokchitser__coeff_growth)
     858        X.init_coeffs(*X._Dokchitser__init_coeff_data)
     859     
     860
     861    if False :
     862        # We have an old version of the Dokchitser class, so convert it
     863        (v, cutoff, w, pari_precode, max_imaginary_part, max_asymp_coeffs) = X._Dokchitser__init
     864
     865        X._init_parameters(max_imaginary_part, max_asymp_coeffs)
     866        X.init_coeffs(v, cutoff, w, pari_precode)
     867
    573868    return X
     869
  • sage/schemes/elliptic_curves/lseries_ell.py

    diff -r fbb0ca1d026e sage/schemes/elliptic_curves/lseries_ell.py
    a b  
    108108        If the curve has too large a conductor, it isn't possible to
    109109        compute with the L-series using this command.  Instead a
    110110        RuntimeError is raised:
     111
    111112            sage: e = EllipticCurve([1,1,0,-63900,-1964465932632])
    112113            sage: L = e.lseries().dokchitser(15)
    113114            Traceback (most recent call last):
    114115            ...
    115             RuntimeError: Unable to create L-series, due to precision or other limits in PARI.
     116            RuntimeError: PARI Error.
     117              ***   at top-level: initLdata(vcoeff_clo
     118              ***                 ^--------------------
     119              *** initLdata: the PARI stack overflows !
     120            ...
    116121        """
    117122        if algorithm == 'magma':
    118123            from sage.interfaces.all import magma
     
    493498            0
    494499            sage: L(1.1)
    495500            0.285491007678148
     501            sage: L.dokchitser().init_coeffs(None, max_imaginary_part = 1)
    496502            sage: L(1.1 + I)
    497503            0.174851377216615 + 0.816965038124457*I
    498504        """