"""
General L-series
AUTHOR:
- William Stein
TODO:
- Symmetric powers (and modular degree -- see trac 9758)
- Triple product L-functions: Gross-Kudla, Zhang, etc -- see the code in triple_prod/triple.py
- Support L-calc L-function
- Make it so we use exactly one GP session for *all* of the Dokchitser L-functions
- Tensor products
- Genus 2 curves, via smalljac and genus2reduction
- Fast L-series of elliptic curves over number fields (not just sqrt(5)), via smalljac
- Inverse of number_of_coefficients function.
"""
import copy, math, types
from sage.all import prime_range, cached_method, sqrt, SR, vector, gcd
from sage.rings.all import ZZ, Integer, QQ, O, ComplexField, CDF, infinity as oo
from sage.arith.all import primes
from sage.rings.rational_field import is_RationalField
from sage.schemes.elliptic_curves.ell_generic import is_EllipticCurve
from psage.ellcurve.lseries.helper import extend_multiplicatively_generic
from sage.misc.all import prod
from sage.modular.abvar.abvar import is_ModularAbelianVariety
from sage.modular.dirichlet import is_DirichletCharacter
from sage.lfunctions.dokchitser import Dokchitser
from sage.modular.modsym.space import is_ModularSymbolsSpace
from sage.modular.abvar.abvar import is_ModularAbelianVariety
from sage.rings.number_field.number_field_base import is_NumberField
import sage.modular.modform.element
from sage.modular.all import Newform
from sage.structure.factorization import Factorization
from sage.misc.mrange import cartesian_product_iterator
I = sqrt(-1)
def prec(s):
"""
Return precision of s, if it has a precision attribute. Otherwise
return 53. This is useful internally in this module.
EXAMPLES::
sage: from psage.lseries.eulerprod import prec
sage: prec(ComplexField(100)(1))
100
sage: prec(RealField(125)(1))
125
sage: prec(1/3)
53
"""
if hasattr(s, 'prec'):
return s.prec()
return 53
def norm(a):
"""
Return the norm of a, for a in either a number field or QQ.
This is a function used internally in this module, mainly because
elements of QQ and ZZ have no norm method.
EXAMPLES::
sage: from psage.lseries.eulerprod import norm
sage: K.<a> = NumberField(x^2-x-1)
sage: norm(a+5)
29
sage: (a+5).norm()
29
sage: norm(17)
17
"""
try:
return a.norm()
except AttributeError:
return a
def tiny(prec):
"""
Return a number that we consider tiny to the given precision prec
in bits. This is used in various places as "zero" to the given
precision for various checks, e.g., of correctness of the
functional equation.
"""
return max(1e-8, 1.0/2**(prec-1))
def prime_below(P):
"""
Return the prime in ZZ below the prime P (or element of QQ).
EXAMPLES::
sage: from psage.lseries.eulerprod import prime_below
sage: K.<a> = NumberField(x^2-x-1)
sage: prime_below(K.prime_above(11))
11
sage: prime_below(K.prime_above(5))
5
sage: prime_below(K.prime_above(3))
3
sage: prime_below(7)
7
"""
try:
return P.smallest_integer()
except AttributeError:
return ZZ(P)
class LSeriesDerivative(object):
"""
The formal derivative of an L-series.
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeries
sage: L = LSeries('delta')
sage: L.derivative()
First derivative of L-function associated to Ramanujan's Delta (a weight 12 cusp form)
sage: L.derivative()(11/2)
0.125386233743526
We directly create an instance of the class (users shouldn't need to do this)::
sage: from psage.lseries.eulerprod import LSeriesDerivative
sage: Ld = LSeriesDerivative(L, 2); Ld
Second derivative of L-function associated to Ramanujan's Delta (a weight 12 cusp form)
sage: type(Ld)
<class 'psage.lseries.eulerprod.LSeriesDerivative'>
"""
def __init__(self, lseries, k):
"""
INPUT:
- lseries -- any LSeries object (derives from LseriesAbstract)
- k -- positive integer
"""
k = ZZ(k)
if k <= 0:
raise ValueError, "k must be a positive integer"
self._lseries = lseries
self._k = k
def __cmp__(self, right):
return cmp((self._lseries,self._k), (right._lseries, right._k))
def __call__(self, s):
"""
Return the value of this derivative at s, which must coerce to a
complex number. The computation is to the same precision as s,
or to 53 bits of precision if s is exact.
As usual, if s has large imaginary part, then there could be
substantial precision loss (a warning is printed in that case).
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeries
sage: L = LSeries(EllipticCurve('389a'))
sage: L1 = L.derivative(1); L1(1)
-1.94715429754927e-20
sage: L1(2)
0.436337613850735
sage: L1(I)
-19.8890471908356 + 31.2633280771869*I
sage: L2 = L.derivative(2); L2(1)
1.51863300057685
sage: L2(I)
134.536162459604 - 62.6542402272310*I
"""
return self._lseries._function(prec(s)).derivative(s, self._k)
def __repr__(self):
"""
String representation of this derivative of an L-series.
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeries
sage: L = LSeries('delta')
sage: L.derivative(1).__repr__()
"First derivative of L-function associated to Ramanujan's Delta (a weight 12 cusp form)"
sage: L.derivative(2).__repr__()
"Second derivative of L-function associated to Ramanujan's Delta (a weight 12 cusp form)"
sage: L.derivative(3).__repr__()
"Third derivative of L-function associated to Ramanujan's Delta (a weight 12 cusp form)"
sage: L.derivative(4).__repr__()
"4-th derivative of L-function associated to Ramanujan's Delta (a weight 12 cusp form)"
sage: L.derivative(2011).__repr__()
"2011-th derivative of L-function associated to Ramanujan's Delta (a weight 12 cusp form)"
"""
k = self._k
if k == 1:
kth = 'First'
elif k == 2:
kth = 'Second'
elif k == 3:
kth = 'Third'
else:
kth = '%s-th'%k
return "%s derivative of %s"%(kth, self._lseries)
def derivative(self, k=1):
"""
Return the k-th derivative of this derivative object.
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeries
sage: f = Newforms(43,2,names='a')[1]; f
q + a1*q^2 - a1*q^3 + (-a1 + 2)*q^5 + O(q^6)
sage: L = LSeries(f); L1 = L.derivative()
sage: L1(1)
0.331674007376949
sage: L(1)
0.620539857407845
sage: L = LSeries(f); L1 = L.derivative(); L1
First derivative of L-series of a degree 2 newform of level 43 and weight 2
sage: L1.derivative()
Second derivative of L-series of a degree 2 newform of level 43 and weight 2
sage: L1.derivative(3)
4-th derivative of L-series of a degree 2 newform of level 43 and weight 2
"""
if k == 0:
return self
return LSeriesDerivative(self._lseries, self._k + k)
class LSeriesParentClass(object):
def __contains__(self, x):
return isinstance(x, (LSeriesAbstract, LSeriesProduct))
def __call__(self, x):
"""
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeries
sage: L = LSeries('zeta')
sage: P = L.parent(); P
All L-series objects and their products
sage: P(L) is L
True
sage: P(L^3)
(Riemann Zeta function viewed as an L-series)^3
sage: (L^3).parent()
All L-series objects and their products
You can make the L-series attached to an object by coercing
that object into the parent::
sage: P(EllipticCurve('11a'))
L-series of Elliptic Curve defined by y^2 + y = x^3 - x^2 - 10*x - 20 over Rational Field
We also allow coercing in 1, because this is very useful for
formal factorizations and products::
sage: P(1)
1
Other numbers do not coerce in, of course::
sage: P(2)
Traceback (most recent call last):
...
TypeError
"""
if isinstance(x, LSeriesAbstract):
return x
elif isinstance(x, LSeriesProduct):
return x
elif x == 1:
return x
else:
try:
return LSeries(x)
except NotImplementedError:
raise TypeError
def __repr__(self):
"""
Return string representation of this parent object.
sage: from psage.lseries.eulerprod import LSeriesParent
sage: LSeriesParent.__repr__()
'All L-series objects and their products'
"""
return "All L-series objects and their products"
def __cmp__(self, right):
"""
Returns equality if right is an instance of
LSeriesParentClass; otherwise compare memory locations.
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeriesParent, LSeriesParentClass
sage: cmp(LSeriesParent, LSeriesParentClass())
0
sage: cmp(LSeriesParent, 7) != 0
True
sage: cmp(7, LSeriesParent) == -cmp(LSeriesParent, 7)
True
"""
if isinstance(right, LSeriesParentClass):
return 0
return cmp(id(self), id(right))
LSeriesParent = LSeriesParentClass()
class LSeriesAbstract(object):
r"""
L-series defined by an Euler product.
The parameters that define the 'shape' of the L-series are:
conductor, hodge_numbers, weight, epsilon, poles, base_field
Let gamma(s) = prod( Gamma((s+h)/2) for h in hodge_numbers ).
Denote this L-series by L(s), and let `L^*(s) = A^s \gamma(s) L(s)`, where
`A = \sqrt(N)/\pi^{d/2}`, where d = len(hodge_numbers) and N = conductor.
Then the functional equation is
Lstar(s) = epsilon * Lstar(weight - s).
To actually use this class we create a derived class that in
addition implements a method _local_factor(P), that takes as input
a prime ideal P of K=base_field, and returns a polynomial, which
is typically the reversed characteristic polynomial of Frobenius
at P of Gal(Kbar/K) acting on the maximal unramified quotient of
some Galois representation. This class automatically computes the
Dirichlet series coefficients `a_n` from the local factors of the
`L`-function.
The derived class may optionally -- and purely as an optimization
-- define a method self._precompute_local_factors(bound,
prec=None), which is typically called before
[_local_factor(P) for P with norm(P) < bound]
is called in the course of various computations.
"""
def __init__(self,
conductor,
hodge_numbers,
weight,
epsilon,
poles,
residues,
base_field,
is_selfdual=True,
prec=53):
"""
INPUT:
- ``conductor`` -- list or number (in a subset of the
positive real numbers); if the conductor is a list, then
each conductor is tried in order (to the precision prec
below) until we find one that works.
- ``hodge_numbers`` -- list of numbers (in a subring of the complex numbers)
- ``weight`` -- number (in a subset of the positive real numbers)
- ``epsilon`` -- number (in a subring of the complex numbers)
- ``poles`` -- list of numbers (in subring of complex numbers); poles of the *completed* L-function
- ``residues`` -- list of residues at each pole given in poles or string "automatic"
- ``base_field`` -- QQ or a number field; local L-factors
correspond to nonzero prime ideals of this field.
- ``is_selfdual`` -- bool (default: True)
- ``prec`` -- integer (default: 53); precision to use when trying to figure
out parameters using the functional equation
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeriesAbstract
sage: L = LSeriesAbstract(conductor=1, hodge_numbers=[0], weight=1, epsilon=1, poles=[1], residues=[-1], base_field=QQ)
sage: type(L)
<class 'psage.lseries.eulerprod.LSeriesAbstract'>
sage: L._conductor
1
sage: L._hodge_numbers
[0]
sage: L._weight
1
sage: L._epsilon
1
sage: L._poles
[1]
sage: L._residues
[-1]
sage: L._base_field
Rational Field
sage: L
Euler Product L-series with conductor 1, Hodge numbers [0], weight 1, epsilon 1, poles [1], residues [-1] over Rational Field
"""
self._anlist = {None:[]}
(self._conductor, self._hodge_numbers, self._weight, self._epsilon,
self._poles, self._residues, self._base_field, self._is_selfdual) = (
conductor, hodge_numbers, weight, epsilon, poles, residues,
base_field, is_selfdual)
v = []
if isinstance(conductor, list):
v.append('_conductor')
if len(hodge_numbers)>0 and isinstance(hodge_numbers[0], list):
v.append('_hodge_numbers')
if isinstance(weight, list):
v.append('_weight')
if len(poles) > 0 and isinstance(poles[0], list):
v.append('_poles')
if isinstance(epsilon, list):
v.append('_epsilon')
if len(v) > 0:
found_params = False
for X in cartesian_product_iterator([getattr(self, attr) for attr in v]):
kwds = dict((v[i],X[i]) for i in range(len(v)))
if self._is_valid_parameters(prec=prec, save=True, **kwds):
found_params = True
break
if not found_params:
raise RuntimeError, "no choice of values for %s works"%(', '.join(v))
def _is_valid_parameters(self, prec=53, save=True, **kwds):
valid = False
try:
old = [(k, getattr(self, k)) for k in kwds.keys()]
for k,v in kwds.iteritems():
setattr(self, k, v)
self._function.clear_cache()
self._function(prec=prec)
try:
self._function(prec=prec)
valid = True
except RuntimeError:
pass
finally:
if not save:
for k, v in old:
setattr(self, k, v)
return valid
def __cmp__(self, right):
if self is right:
return 0
for a in ['degree', 'weight', 'conductor', 'epsilon', 'base_field']:
c = cmp(getattr(self, a)(), getattr(right, a)())
if c: return c
c = cmp(type(self), type(right))
return self._cmp(right)
def _cmp(self, right):
raise TypeError
def parent(self):
"""
Return parent of this L-series, which is the collection of all
L-series and their products.
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeries
sage: L = LSeries('delta'); P = L.parent(); P
All L-series objects and their products
sage: L in P
True
"""
return LSeriesParent
def __pow__(self, n):
"""
Return the n-th power of this L-series, where n can be any integer.
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeries
sage: L = LSeries('delta');
sage: L3 = L^3; L3
(L-function associated to Ramanujan's Delta (a weight 12 cusp form))^3
sage: L3(1)
0.0000524870430366548
sage: L(1)^3
0.0000524870430366548
sage: M = L^(-3); M(1)
19052.3211471761
sage: L(1)^(-3)
19052.3211471761
Higher precision::
sage: M = L^(-3); M(RealField(100)(1))
19052.321147176093380952680193
sage: L(RealField(100)(1))^(-3)
19052.321147176093380952680193
Special case -- 0th power -- is not allowed::
sage: L^0
Traceback (most recent call last):
...
ValueError: product must be nonempty
"""
return LSeriesProduct([(self, ZZ(n))])
def __mul__(self, right):
"""
Multiply two L-series, or an L-series times a formal product of L-series.
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeries
sage: d = LSeries('delta'); z = LSeries('zeta')
sage: d * z
(Riemann Zeta function viewed as an L-series) * (L-function associated to Ramanujan's Delta (a weight 12 cusp form))
sage: d * (d * z)
(Riemann Zeta function viewed as an L-series) * (L-function associated to Ramanujan's Delta (a weight 12 cusp form))^2
"""
if isinstance(right, LSeriesAbstract):
return LSeriesProduct([(self, 1), (right, 1)])
elif isinstance(right, LSeriesProduct):
return right * self
raise TypeError
def __div__(self, right):
"""
Divide two L-series or formal L-series products.
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeries
sage: d = LSeries('delta'); z = LSeries('zeta')
sage: d / z
(Riemann Zeta function viewed as an L-series)^-1 * (L-function associated to Ramanujan's Delta (a weight 12 cusp form))
sage: d / (z^3)
(Riemann Zeta function viewed as an L-series)^-3 * (L-function associated to Ramanujan's Delta (a weight 12 cusp form))
"""
if isinstance(right, LSeriesAbstract):
return LSeriesProduct([(self, 1), (right, -1)])
elif isinstance(right, LSeriesProduct):
return LSeriesProduct(Factorization([(self, 1)]) / right._factorization)
raise TypeError
def conductor(self):
"""
Return the conductor of this L-series.
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeries
sage: LSeries('zeta').conductor()
1
sage: LSeries('delta').conductor()
1
sage: LSeries(EllipticCurve('11a')).conductor()
11
sage: LSeries(Newforms(33)[0]).conductor()
33
sage: LSeries(DirichletGroup(37).0).conductor()
37
sage: LSeries(kronecker_character(7)).conductor()
28
sage: kronecker_character(7).conductor()
28
sage: L = LSeries(EllipticCurve('11a3').base_extend(QQ[sqrt(2)]), prec=5)
sage: L.conductor().factor()
2^6 * 11^2
"""
return self._conductor
def hodge_numbers(self):
"""
Return the Hodge numbers of this L-series. These define the local Gamma factors.
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeries
sage: LSeries('zeta').hodge_numbers()
[0]
sage: LSeries('delta').hodge_numbers()
[0, 1]
sage: LSeries(EllipticCurve('11a')).hodge_numbers()
[0, 1]
sage: LSeries(Newforms(43,names='a')[1]).hodge_numbers()
[0, 1]
sage: LSeries(DirichletGroup(37).0).hodge_numbers()
[1]
sage: LSeries(EllipticCurve(QQ[sqrt(-1)],[1,2]), prec=5).hodge_numbers() # long time
[0, 0, 1, 1]
"""
return list(self._hodge_numbers)
def weight(self):
"""
Return the weight of this L-series.
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeries
sage: LSeries('zeta').weight()
1
sage: LSeries('delta').weight()
12
sage: LSeries(EllipticCurve('389a')).weight()
2
sage: LSeries(Newforms(43,names='a')[1]).weight()
2
sage: LSeries(Newforms(6,4)[0]).weight()
4
sage: LSeries(DirichletGroup(37).0).weight()
1
sage: L = LSeries(EllipticCurve('11a3').base_extend(QQ[sqrt(2)]),prec=5); L.weight()
2
"""
return self._weight
def poles(self):
"""
Poles of the *completed* L-function with the extra Gamma
factors included.
WARNING: These are not just the poles of self.
OUTPUT:
- list of numbers
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeries
sage: LSeries('zeta').poles()
[0, 1]
sage: LSeries('delta').poles()
[]
"""
return list(self._poles)
def residues(self, prec=None):
"""
Residues of the *completed* L-function at each pole.
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeries
sage: LSeries('zeta').residues()
[9, 8]
sage: LSeries('delta').residues()
[]
sage: v = LSeries('zeta').residues()
sage: v.append(10)
sage: LSeries('zeta').residues()
[9, 8]
The residues of the Dedekind Zeta function of a field are dynamically computed::
sage: K.<a> = NumberField(x^2 + 1)
sage: L = LSeries(K); L
Dedekind Zeta function of Number Field in a with defining polynomial x^2 + 1
If you just call residues you get back that they are automatically computed::
sage: L.residues()
'automatic'
But if you call with a specific precision, they are computed using that precision::
sage: L.residues(prec=53)
[-0.886226925452758]
sage: L.residues(prec=200)
[-0.88622692545275801364908374167057259139877472806119356410690]
"""
if self._residues == 'automatic':
if prec is None:
return self._residues
else:
C = ComplexField(prec)
return [C(a) for a in self._function(prec=prec).gp()('Lresidues')]
else:
return list(self._residues)
def base_field(self):
"""
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeries
sage: LSeries('zeta').base_field()
Rational Field
sage: L = LSeries(EllipticCurve('11a3').base_extend(QQ[sqrt(2)]), prec=5); L.base_field()
Number Field in sqrt2 with defining polynomial x^2 - 2
"""
return self._base_field
def epsilon(self, prec=None):
"""
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeries
sage: LSeries('zeta').epsilon()
1
sage: LSeries('delta').epsilon()
1
sage: LSeries(EllipticCurve('389a')).epsilon()
1
sage: LSeries(EllipticCurve('37a')).epsilon()
-1
sage: LSeries(Newforms(389,names='a')[1]).epsilon()
-1
sage: LSeries(Newforms(6,4)[0]).epsilon()
1
sage: LSeries(DirichletGroup(7).0).epsilon()
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)
sage: LSeries(DirichletGroup(7).0).epsilon(prec=53)
0.386513572759156 + 0.922283718859307*I
In this example, the epsilon factor is computed when the curve
is created. The prec parameter determines the floating point precision
used in computing the epsilon factor::
sage: L = LSeries(EllipticCurve('11a3').base_extend(QQ[sqrt(2)]), prec=5); L.epsilon()
-1
Here is extra confirmation that the rank is really odd over the quadratic field::
sage: EllipticCurve('11a').quadratic_twist(2).rank()
1
We can compute with the L-series too::
sage: L(RealField(5)(2))
0.53
For L-functions of newforms with nontrivial character, the
epsilon factor is harder to find (we don't have a good
algorithm implemented to find it) and might not even be 1 or
-1, so it is set to 'solve'. In this case, the functional
equation is used to determine the solution.::
sage: f = Newforms(kronecker_character_upside_down(7),3)[0]; f
q - 3*q^2 + 5*q^4 + O(q^6)
sage: L = LSeries(f)
sage: L.epsilon()
'solve'
sage: L(3/2)
0.332981771482934
sage: L.epsilon()
1
Here is an example with nontrivial character::
sage: f = Newforms(DirichletGroup(7).0, 5, names='a')[0]; f
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)
sage: L = LSeries(f)
First trying to evaluate with the default (53 bits) of
precision fails, since for some reason (that I do not
understand) the program is unable to find a valid epsilon
factor::
sage: L(0)
Traceback (most recent call last):
...
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
However, when we evaluate to 100 bits of precision it works::
sage: L(RealField(100)(0))
0
The epsilon factor is *not* known to infinite precision::
sage: L.epsilon()
'solve'
But it is now known to 100 bits of precision, and here it is::
sage: L.epsilon(100)
0.42563106101692403875896879406 - 0.90489678963824790765479396740*I
When we try to compute to higher precision, again Sage solves for the epsilon factor
numerically::
sage: L(RealField(150)(1))
0.26128389551787271923496480408992971337929665 - 0.29870133769674001421149135036267324347896657*I
And now it is known to 150 bits of precision. Notice that
this is consistent with the value found above, and has
absolute value (very close to) 1.
sage: L.epsilon(150)
0.42563106101692403875896879406038776338921622 - 0.90489678963824790765479396740501409301704122*I
sage: abs(L.epsilon(150))
1.0000000000000000000000000000000000000000000
"""
if self._epsilon == 'solve' or (
hasattr(self._epsilon, 'prec') and (prec is None or self._epsilon.prec() < prec)):
return 'solve'
if prec is not None:
C = ComplexField(prec)
if isinstance(self._epsilon, list):
return [C(x) for x in self._epsilon]
return C(self._epsilon)
return self._epsilon
def is_selfdual(self):
"""
Return True if this L-series is self dual; otherwise, return False.
EXAMPLES::
Many L-series are self dual::
sage: from psage.lseries.eulerprod import LSeries
sage: LSeries('zeta').is_selfdual()
True
sage: LSeries('delta').is_selfdual()
True
sage: LSeries(Newforms(6,4)[0]).is_selfdual()
True
Nonquadratic characters have non-self dual L-series::
sage: LSeries(DirichletGroup(7).0).is_selfdual()
False
sage: LSeries(kronecker_character(7)).is_selfdual()
True
Newforms with non-quadratic characters also have non-self dual L-seris::
sage: L = LSeries(Newforms(DirichletGroup(7).0, 5, names='a')[0]); L.is_selfdual()
False
"""
return self._is_selfdual
@cached_method
def degree(self):
"""
Return the degree of this L-function, which is by definition the number of Gamma
factors (e.g., the number of Hodge numbers) divided by the degree of the base field.
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeries
sage: LSeries('zeta').degree()
1
sage: LSeries(DirichletGroup(5).0).degree()
1
sage: LSeries(EllipticCurve('11a')).degree()
2
The L-series attached to this modular symbols space of dimension 2 is a product
of 2 degree 2 L-series, hence has degree 4::
sage: M = ModularSymbols(43,2,sign=1).cuspidal_subspace()[1]; M.dimension()
2
sage: L = LSeries(M); L
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
sage: L.factor()
(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)
sage: L.degree()
4
sage: x = var('x'); K.<a> = NumberField(x^2-x-1); LSeries(EllipticCurve([0,-a,a,0,0])).degree()
2
"""
n = len(self.hodge_numbers())
d = self.base_field().degree()
assert n % d == 0, "degree of base field must divide the number of Hodge numbers"
return n//d
def twist(self, chi, conductor=None, epsilon=None, prec=53):
r"""
Return the quadratic twist of this L-series by the character chi, which
must be a character of self.base_field(). Thus chi should take as input
prime ideals (or primes) of the ring of integers of the base field, and
output something that can be coerced to the complex numbers.
INPUT:
- `\chi` -- 1-dimensional character
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeries
sage: E = EllipticCurve('11a')
sage: L = LSeries(E); L
L-series of Elliptic Curve defined by y^2 + y = x^3 - x^2 - 10*x - 20 over Rational Field
sage: L3 = L.twist(DirichletGroup(3).0); L3
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
sage: L3._chi
Dirichlet character modulo 3 of conductor 3 mapping 2 |--> -1
sage: L3._L
L-series of Elliptic Curve defined by y^2 + y = x^3 - x^2 - 10*x - 20 over Rational Field
sage: L3(1)
1.68449633297548
sage: F = E.quadratic_twist(-3)
sage: L3.conductor()
99
sage: F.conductor()
99
sage: F.lseries()(1)
1.68449633297548
sage: L3.anlist(20)
[0, 1, 2, 0, 2, -1, 0, -2, 0, 0, -2, -1, 0, 4, -4, 0, -4, 2, 0, 0, -2]
sage: F.anlist(20)
[0, 1, 2, 0, 2, -1, 0, -2, 0, 0, -2, -1, 0, 4, -4, 0, -4, 2, 0, 0, -2]
sage: L3.anlist(1000) == F.anlist(1000)
True
sage: L3.local_factor(11)
T + 1
A higher degree twist::
sage: L = LSeries(EllipticCurve('11a'))
sage: L5 = L.twist(DirichletGroup(5).0); L5
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
sage: L5(1)
1.28009593569230 - 0.681843202124309*I
sage: L5.epsilon()
'solve'
sage: L5.epsilon(53)
0.989989082587826 + 0.141143956147310*I
sage: L5.conductor()
275
sage: L5.taylor_series(center=1, degree=3)
1.28009593569230 - 0.681843202124309*I + (-0.536450338806282 + 0.166075270978779*I)*z + (0.123743053129226 + 0.320802890011298*I)*z^2 + O(z^3)
WARNING!! Twisting is not implemented in full generality when
the conductors are not coprime. One case where we run into
trouble is when twisting lowers the level of a newform. Below
we take the form of level 11 and weight 2, twist it by the
character chi of conductor 3 to get a form of level 99. Then
we take the L-series of the level 99 form, and twist that by
chi, which should be the L-series attached to the form of
level 11. Unfortunately, our code for working out the local
L-factors doesn't succeed in this case, hence the local factor
is wrong, so the functional equation is not satisfied.
sage: f = Newform('11a'); f
q - 2*q^2 - q^3 + 2*q^4 + q^5 + O(q^6)
sage: L = LSeries(f)
sage: chi = DirichletGroup(3).0
sage: Lc = L.twist(chi); Lc
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
sage: Lc.anlist(20)
[0, 1, 2, 0, 2, -1, 0, -2, 0, 0, -2, -1, 0, 4, -4, 0, -4, 2, 0, 0, -2]
sage: g = Newform('99d'); g
q + 2*q^2 + 2*q^4 - q^5 + O(q^6)
sage: list(g.qexp(20))
[0, 1, 2, 0, 2, -1, 0, -2, 0, 0, -2, -1, 0, 4, -4, 0, -4, 2]
sage: Lt = Lc.twist(chi, conductor=11)
Traceback (most recent call last):
...
RuntimeError: no choice of values for _epsilon works
sage: Lt = Lc.twist(chi,conductor=11, epsilon=1)
sage: Lt(1)
Traceback (most recent call last):
...
RuntimeError: invalid L-series parameters: functional equation not satisfied
This is because the local factor is wrong::
sage: Lt.local_factor(3)
1
sage: L.local_factor(3)
3*T^2 + T + 1
"""
return LSeriesTwist(self, chi=chi, conductor=conductor, epsilon=epsilon, prec=prec)
@cached_method
def local_factor(self, P, prec=None):
"""
Return the local factor of the L-function at the prime P of
self._base_field. The result is cached.
INPUT:
- a prime P of the ring of integers of the base_field
- prec -- None or positive integer (bits of precision)
OUTPUT:
- a polynomial, e.g., something like "1-a*T+p*T^2".
EXAMPLES:
You must overload this in the derived class::
sage: from psage.lseries.eulerprod import LSeriesAbstract
sage: L = LSeriesAbstract(conductor=1, hodge_numbers=[0], weight=1, epsilon=1, poles=[1], residues=[-1], base_field=QQ)
sage: L.local_factor(2)
Traceback (most recent call last):
...
NotImplementedError: must be implemented in the derived class
"""
return self._local_factor(P, prec=prec)
def _local_factor(self, P, prec=None):
"""
Compute local factor at prime P. This must be overwritten in the derived class.
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeriesAbstract
sage: L = LSeriesAbstract(conductor=1, hodge_numbers=[0], weight=1, epsilon=1, poles=[1], residues=[-1], base_field=QQ)
sage: L.local_factor(2)
Traceback (most recent call last):
...
NotImplementedError: must be implemented in the derived class
"""
raise NotImplementedError, "must be implemented in the derived class"
def __repr__(self):
"""
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeriesAbstract
sage: L = LSeriesAbstract(conductor=1, hodge_numbers=[0], weight=1, epsilon=1, poles=[1], residues=[-1], base_field=QQ)
sage: L.__repr__()
'Euler Product L-series with conductor 1, Hodge numbers [0], weight 1, epsilon 1, poles [1], residues [-1] over Rational Field'
"""
return "Euler Product L-series with conductor %s, Hodge numbers %s, weight %s, epsilon %s, poles %s, residues %s over %s"%(
self._conductor, self._hodge_numbers, self._weight, self.epsilon(), self._poles, self._residues, self._base_field)
def _precompute_local_factors(self, bound, prec):
"""
Derived classes may use this as a 'hint' that _local_factors
will soon get called for primes of norm less than the bound.
In the base class, this is a no-op, and it is not necessary to
overload this class.
INPUT:
- ``bound`` -- integer
- ``prec`` -- integer
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeriesAbstract
sage: L = LSeriesAbstract(conductor=1, hodge_numbers=[0], weight=1, epsilon=1, poles=[1], residues=[-1], base_field=QQ)
sage: L._precompute_local_factors(100, 53)
"""
pass
def _primes_above(self, p):
"""
Return the primes of the ring of integers of the base field above the integer p.
INPUT:
- p -- prime integer (no type checking necessarily done)
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeriesAbstract
sage: L = LSeriesAbstract(conductor=1, hodge_numbers=[0], weight=1, epsilon=1, poles=[1], residues=[-1], base_field=QQ)
sage: L._primes_above(3)
[3]
sage: L = LSeriesAbstract(conductor=1, hodge_numbers=[0], weight=1, epsilon=1, poles=[1], residues=[-1], base_field=QQ[sqrt(-1)])
sage: L._primes_above(5)
[Fractional ideal (I + 2), Fractional ideal (-I + 2)]
sage: L._primes_above(3)
[Fractional ideal (3)]
"""
K = self._base_field
if is_RationalField(K):
return [p]
else:
return K.primes_above(p)
def zeros(self, n):
"""
Return the imaginary parts of the first n nontrivial zeros of
this L-function on the critical line in the upper half plane,
as 32-bit real numbers.
INPUT:
- n -- nonnegative integer
EXAMPLES::
"""
return self._lcalc().zeros(n)
def _lcalc(self):
"""
Return Rubinstein Lcalc object attached to this L-series. This
is useful both for evaluating the L-series, especially when
the imaginary part is large, and for computing the zeros in
the critical strip.
EXAMPLES::
"""
raise NotImplementedError
def anlist(self, bound, prec=None):
"""
Return list `v` of Dirichlet series coefficients `a_n`for `n` up to and including bound,
where `v[n] = a_n`. If at least one of the coefficients is ambiguous, e.g., the
local_factor method returns a list of possibilities, then this function instead
returns a generator over all possible `a_n` lists.
In particular, we include v[0]=0 as a convenient place holder
to avoid having `v[n-1] = a_n`.
INPUT:
- ``bound`` -- nonnegative integer
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeries; L = LSeries('zeta')
sage: L.anlist(30)
[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]
sage: from psage.lseries.eulerprod import LSeries; L = LSeries(EllipticCurve('11a'))
sage: L.anlist(30)
[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]
sage: K.<a> = NumberField(x^2-x-1); L = LSeries(EllipticCurve([0,-a,a,0,0]))
sage: L.anlist(30)
[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]
"""
if len(self._anlist[None]) > bound:
if prec is None:
return self._anlist[None][:bound+1]
else:
C = ComplexField(prec)
return [C(a) for a in self._anlist[None]]
if prec is not None:
t = [z for z in self._anlist.iteritems() if z[0] >= prec and len(z[1]) > bound]
if len(t) > 0:
C = ComplexField(prec)
return [C(a) for a in t[0][1][:bound+1]]
self._precompute_local_factors(bound+1, prec=prec)
compute_anlist_multiple = False
LF = []
for p in prime_range(bound+1):
lf = []
for P in self._primes_above(p):
if norm(P) <= bound:
F = self._local_factor(P, prec)
if isinstance(F, list):
compute_anlist_multiple = True
lf.append(F)
LF.append((p, lf))
coefficients = self._compute_anlist(LF, bound, prec)
if not compute_anlist_multiple:
coefficients = list(coefficients)[0]
self._anlist[prec] = coefficients
return coefficients
def _compute_anlist(self, LF, bound, prec):
"""
Iterator over possible anlists, given LF, bound, and prec.
INPUT:
- ``LF`` -- list of pairs (p, [local factors (or lists of them) at primes over p])
- ``bound`` -- positive integer
- ``prec`` -- positive integer (bits of precision)
"""
K = self._base_field
coefficients = [0,1] + [0]*(bound-1)
for i, (p, v) in enumerate(LF):
if len(v) > 0:
some_list = False
for j, lf in enumerate(v):
if isinstance(lf, list):
some_list = True
for f in list(lf):
LF0 = copy.deepcopy(LF)
LF0[i][1][j] = f
for z in self._compute_anlist(LF0, bound, prec):
yield z
if some_list:
return
f = prod(v)
accuracy_p = int(math.floor(math.log(bound)/math.log(p))) + 1
T = f.parent().gen()
series_p = (f + O(T**accuracy_p))**(-1)
for j in range(1, accuracy_p):
coefficients[p**j] = series_p[j]
extend_multiplicatively_generic(coefficients)
yield list(coefficients)
def _symbolic_(self, R, bound=10, prec=None):
"""
Convert self into the symbolic ring as a truncated Dirichleter series, including
terms up to `n^s` where n=bound.
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeries
sage: L = LSeries(EllipticCurve('37a'))
sage: SR(L)
-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
sage: L._symbolic_(SR, 20)
-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
"""
s = R.var('s')
a = self.anlist(bound, prec)
return sum(a[n]/n**s for n in range(1,bound+1))
def __call__(self, s):
"""
Return value of this L-function at s. If s is a real or
complex number to prec bits of precision, then the result is
also computed to prec bits of precision. If s has infinite or
unknown precision, then the L-value is computed to 53 bits of
precision.
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeries
sage: L = LSeries(EllipticCurve('37a'))
sage: L(1)
0
sage: L(2)
0.381575408260711
sage: z = L(RealField(100)(2)); z
0.38157540826071121129371040958
sage: z.prec()
100
sage: L(RealField(150)(2))
0.38157540826071121129371040958008663667709753
WARNING: There will be precision loss (with a warning) if the imaginary
part of the input is large::
sage: L = LSeries('zeta')
sage: L(1/2 + 40*I)
verbose -1 (...: dokchitser.py, __call__) Warning: Loss of 14 decimal digits due to cancellation
0.793046013671137 - 1.04127377821427*I
sage: L(ComplexField(200)(1/2 + 40*I))
verbose -1 (...: dokchitser.py, __call__) Warning: Loss of 14 decimal digits due to cancellation
0.79304495256192867196489258889793696080572220439302833315881 - 1.0412746146510650200518905953910554313275550685861559488384*I
An example with a pole::
sage: L = LSeries('zeta')
sage: L(2) # good
1.64493406684823
sage: L(1) # a pole!
Traceback (most recent call last):
...
ZeroDivisionError: pole at 1
"""
if s in self._poles:
raise ZeroDivisionError, "pole at %s"%s
return self._function(prec(s))(s)
def derivative(self, k=1):
"""
Return the k-th derivative of self.
INPUT:
- k -- (default: 1) nonnegative integer
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeries
sage: L = LSeries('zeta')
sage: L.derivative()
First derivative of Riemann Zeta function viewed as an L-series
We numerically approximate the derivative at two points and compare
with evaluating the derivative object::
sage: eps=1e-10; (L(2+eps) - L(2))/eps
-0.937547817159157
sage: Lp = L.derivative(); Lp(2)
-0.937548254315844
sage: eps=1e-10; (L(2+I+eps) - L(2+I))/eps
0.0624900131640516 + 0.489033813444451*I
sage: Lp(2+I)
0.0624900021906470 + 0.489033591679899*I
Higher derivatives::
sage: L.derivative(2)
Second derivative of Riemann Zeta function viewed as an L-series
sage: L.derivative(2)(2)
1.98928023429890
sage: L.derivative(3)
Third derivative of Riemann Zeta function viewed as an L-series
sage: L.derivative(4)
4-th derivative of Riemann Zeta function viewed as an L-series
sage: L.derivative(5)
5-th derivative of Riemann Zeta function viewed as an L-series
Derivative of derivative::
sage: L.derivative().derivative()
Second derivative of Riemann Zeta function viewed as an L-series
Using the derivative function in Sage works::
sage: derivative(L)
First derivative of Riemann Zeta function viewed as an L-series
sage: derivative(L,2)
Second derivative of Riemann Zeta function viewed as an L-series
"""
if k == 0:
return self
return LSeriesDerivative(self, k)
def taylor_series(self, center=None, degree=6, variable='z', prec=53):
"""
Return the Taylor series expansion of self about the given
center point to the given degree in the specified variable
numerically computed to the precision prec in bits. If the
center is not specified it defaults to weight / 2.
INPUT:
- ``center`` -- None or number that coerces to the complex numbers
- ``degree`` -- integer
- ``variable`` -- string or symbolic variable
- ``prec`` -- positive integer (floating point bits of precision)
EXAMPLES:::
sage: from psage.lseries.eulerprod import LSeries; L = LSeries('zeta')
sage: L.taylor_series()
-1.46035450880959 - 3.92264613920915*z - 8.00417850696433*z^2 - 16.0005515408865*z^3 - 31.9998883216853*z^4 - 64.0000050055172*z^5 + O(z^6)
sage: RealField(53)(zeta(1/2))
-1.46035450880959
sage: L.taylor_series(center=2, degree=4, variable='t', prec=30)
1.6449341 - 0.93754825*t + 0.99464012*t^2 - 1.0000243*t^3 + O(t^4)
sage: RealField(30)(zeta(2))
1.6449341
"""
if center is None:
center = ComplexField(prec)(self._weight)/2
return self._function(prec).taylor_series(center, degree, variable)
def analytic_rank(self, tiny=1e-8, prec=53):
center = ComplexField(prec)(self._weight) / 2
degree = 4
while True:
f = self.taylor_series(center, degree, prec=prec)
i = 0
while i < degree:
if abs(f[i]) > tiny:
return i
i += 1
degree += 2
@cached_method
def _function(self, prec=53, T=1.2):
"""
Return Dokchitser object that allows for computation of this
L-series computed to enough terms so that the functional
equation checks out with the given value of T and precision.
This is used behind the scenes for evaluation and computation
of Taylor series.
"""
eps = self.epsilon(prec)
return self._dokchitser(prec, eps, T=T)
def _dokchitser_unitialized(self, prec, epsilon):
if epsilon == 'solve':
eps = 'X'
else:
eps = epsilon
return Dokchitser(conductor = self.conductor(), gammaV = self.hodge_numbers(), weight = self.weight(),
eps = eps, poles = self.poles(), residues = self.residues(),
prec = prec)
def number_of_coefficients(self, prec=53, T=1.2):
"""
Return the number of Dirichlet series coefficients that will
be needed in order to evaluate this L-series (near the real
line) to prec bits of precision and have the functional
equation test pass with the given value of T.
INPUT:
- prec -- integer
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeries
sage: L = LSeries(DirichletGroup(5).0)
sage: L.number_of_coefficients(20)
8
sage: L.number_of_coefficients()
11
sage: L.number_of_coefficients(1000)
43
"""
return self._dokchitser_unitialized(prec, self.epsilon(prec)).num_coeffs(T)
def _dokchitser(self, prec, epsilon, T=1.2):
L = self._dokchitser_unitialized(prec, epsilon)
n = L.num_coeffs(T=T)
X = self.anlist(n, prec)
if isinstance(X, types.GeneratorType):
coeff_lists = X
else:
coeff_lists = [X]
tiny0 = tiny(prec)
for coeffs in coeff_lists:
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:]])
if self.is_selfdual():
L.init_coeffs('a(k)', pari_precode = s)
else:
L._Dokchitser__init = True
L._gp_eval(s)
L._gp_eval('initLdata("a(k)",1,"conj(a(k))")')
if epsilon == 'solve':
cmd = "sgneq = Vec(checkfeq()); sgn = -sgneq[2]/sgneq[1]; sgn"
epsilon = ComplexField(prec)(L._gp_eval(cmd))
if abs(abs(epsilon)-1) > tiny0:
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)
if epsilon == 1:
self._epsilon = 1
elif epsilon == -1:
self._epsilon = -1
else:
self._epsilon = epsilon
fe = L.check_functional_equation()
if abs(fe) <= tiny0:
self._anlist[prec] = coeffs
return L
else:
pass
raise RuntimeError, "invalid L-series parameters: functional equation not satisfied"
def check_functional_equation(self, T, prec=53):
return self._function(prec=prec).check_functional_equation(T)
class LSeriesProductEvaluator(object):
def __init__(self, factorization, prec):
self._factorization = factorization
self._prec = prec
def __call__(self, s):
try:
v = self._functions
except AttributeError:
self._functions = [(L._function(self._prec),e) for L,e in self._factorization]
v = self._functions
return prod(f(s)**e for f,e in v)
class LSeriesProduct(object):
"""
A formal product of L-series.
"""
def __init__(self, F):
"""
INPUT:
- `F` -- list of pairs (L,e) where L is an L-function and e is a nonzero integer.
"""
if not isinstance(F, Factorization):
F = Factorization(F)
F.sort(cmp)
if len(F) == 0:
raise ValueError, "product must be nonempty"
self._factorization = F
def __cmp__(self, right):
return cmp(self._factorization, right._factorization)
def is_selfdual(self):
"""
Return True if every factor of self is self dual; otherwise, return False.
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeries
sage: L1 = LSeries('zeta'); L1.is_selfdual()
True
sage: L2 = LSeries(DirichletGroup(9).0); L2.is_selfdual()
False
sage: (L1*L1).is_selfdual()
True
sage: (L1*L2).is_selfdual()
False
sage: (L2*L2).is_selfdual()
False
"""
is_selfdual = set(L.is_selfdual() for L,e in self._factorization)
if len(is_selfdual) > 1:
return False
else:
return list(is_selfdual)[0]
def conductor(self):
"""
Return the conductor of this product, which we define to be
the product with multiplicities of the conductors of the
factors of self.
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeries
sage: J = J0(33); L = LSeries(J)
sage: L.conductor()
3993
sage: L.conductor().factor()
3 * 11^3
sage: J.decomposition()
[
Simple abelian subvariety 11a(1,33) of dimension 1 of J0(33),
Simple abelian subvariety 11a(3,33) of dimension 1 of J0(33),
Simple abelian subvariety 33a(1,33) of dimension 1 of J0(33)
]
Of course, the conductor as we have defined it need not be an integer::
sage: L = LSeries(J[0])/LSeries(J[2])^2; L
(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
sage: L.conductor()
1/99
"""
return prod(L.conductor()**e for L, e in self._factorization)
def __pow__(self, n):
"""
Return the n-th power of this formal L-series product, where n can be any integer.
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeries
sage: L = LSeries('zeta') * LSeries(DirichletGroup(9).0)
sage: L(2)
1.80817853715812 + 0.369298119218816*I
sage: L = LSeries('zeta') * LSeries(DirichletGroup(9).0)
sage: L3 = L^3; L3
(Riemann Zeta function viewed as an L-series)^3 * (L-series attached to Dirichlet character modulo 9 of conductor 9 mapping 2 |--> zeta6)^3
sage: L3(2)
5.17205298762567 + 3.57190597873829*I
sage: CC((zeta(2)*LSeries(DirichletGroup(9).0)(2))^3)
5.17205298762567 + 3.57190597873829*I
sage: L^(-1999)
(Riemann Zeta function viewed as an L-series)^-1999 * (L-series attached to Dirichlet character modulo 9 of conductor 9 mapping 2 |--> zeta6)^-1999
sage: (L^(-1999)) (2)
8.90248311986228e-533 - 6.20123089437732e-533*I
"""
return LSeriesProduct(self._factorization**ZZ(n))
def __mul__(self, right):
if isinstance(right, LSeriesAbstract):
return LSeriesProduct(self._factorization * Factorization([(right,1)]))
elif isinstance(right, LSeriesProduct):
return LSeriesProduct(self._factorization * right._factorization)
else:
raise TypeError
def __div__(self, right):
if isinstance(right, LSeriesAbstract):
return LSeriesProduct(self._factorization * Factorization([(right,-1)]))
elif isinstance(right, LSeriesProduct):
return LSeriesProduct(self._factorization / right._factorization)
raise TypeError
def parent(self):
return LSeriesParent
def factor(self):
return self._factorization
def hodge_numbers(self):
"""
Return the Hodge numbers of this product of L-series.
"""
def degree(self):
return sum(e*L.degree() for L,e in self._factorization)
def local_factor(self, P, prec=None):
return prod(L.local_factor(P,prec)**e for L,e in self._factorization)
def __getitem__(self, *args, **kwds):
return self._factorization.__getitem__(*args, **kwds)
def __repr__(self):
return self.factor().__repr__()
def __call__(self, s):
return self._function(prec(s))(s)
def derivative(self, k=1):
raise NotImplementedError
def taylor_series(self, center=None, degree=6, variable='z', prec=53):
"""
EXAMPLE::
sage: from psage.lseries.eulerprod import LSeries
sage: L1 = LSeries('zeta'); L2 = LSeries('delta')
sage: f = L1 * L2; f
(Riemann Zeta function viewed as an L-series) * (L-function associated to Ramanujan's Delta (a weight 12 cusp form))
sage: f.taylor_series(center=2, degree=4, variable='w', prec=30)
0.24077647 + 0.10066485*w + 0.061553731*w^2 - 0.041923238*w^3 + O(w^4)
sage: L1.taylor_series(2, 4, 'w', 30) * L2.taylor_series(2, 4, 'w', 30)
0.24077647 + 0.10066485*w + 0.061553731*w^2 - 0.041923238*w^3 + O(w^4)
sage: f = L1 / L2; f
(Riemann Zeta function viewed as an L-series) * (L-function associated to Ramanujan's Delta (a weight 12 cusp form))^-1
sage: f.taylor_series(center=2, degree=4, variable='w', prec=30)
11.237843 - 17.508629*w + 21.688182*w^2 - 24.044641*w^3 + O(w^4)
sage: L1.taylor_series(2, 4, 'w', 30) / L2.taylor_series(2, 4, 'w', 30)
11.237843 - 17.508629*w + 21.688182*w^2 - 24.044641*w^3 + O(w^4)
"""
return prod(L.taylor_series(center, degree, variable, prec)**e for L,e in self._factorization)
def analytic_rank(self, prec=53):
"""
Return sum of the order of vanishing counted with
multiplicities of each factors at their center point.
WARNING: The analytic rank is computed numerically, so is
definitely not provably correct.
EXAMPLES::
We compute the analytic rank of the non-simple 32-dimensional modular abelian variety `J_0(389)`.
sage: from psage.lseries.eulerprod import LSeries
sage: M = ModularSymbols(389,sign=1).cuspidal_subspace()
sage: L = LSeries(M); L
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
We first attempt computation of the analytic rank with the default of 53 bits precision::
sage: L.analytic_rank()
Traceback (most recent call last):
...
RuntimeError: invalid L-series parameters: functional equation not satisfied
The above failed because trying to compute one of the degree
20 newforms resulting in some internal error when double
checking the functional equation. So we try with slightly more precision::
sage: L.analytic_rank(70)
13
This works, since the factors have dimensions 1,2,3,6,20, and
the one of degree 1 has rank 2, the ones of degree 2,3,6 have
rank 2,3,6, respectively, and the one of degree 20 has rank 0::
sage: 2*1 + 2 + 3 + 6
13
"""
return sum(e*L.analytic_rank(prec=prec) for L,e in self._factorization)
def weight(self):
"""
Return the weight of this L-series, which is the sum of the weights
of the factors counted with multiplicity.
"""
return sum(e*L.weight() for L,e in self._factorization)
@cached_method
def _function(self, prec=53):
"""
Return Dokchitser object that allows for computation of this
L-series. This is used behind the scenes for evaluation and
computation of Taylor series.
"""
return LSeriesProductEvaluator(self._factorization, prec)
class LSeriesZeta(LSeriesAbstract):
"""
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeries
sage: L = LSeries('zeta'); L
Riemann Zeta function viewed as an L-series
sage: L(2)
1.64493406684823
sage: L(2.000000000000000000000000000000000000000000000000000)
1.64493406684822643647241516664602518921894990120680
sage: zeta(2.000000000000000000000000000000000000000000000000000)
1.64493406684822643647241516664602518921894990120680
sage: L.local_factor(3)
-T + 1
sage: L.local_factor(5)
-T + 1
sage: L.anlist(30)
[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]
"""
def __init__(self):
LSeriesAbstract.__init__(self, conductor=1, hodge_numbers=[0], weight=1, epsilon=1,
poles=[0,1], residues=[9,8], base_field=QQ)
T = ZZ['T'].gen()
self._lf = 1 - T
def _cmp(self, right):
return 0
def _local_factor(self, P, prec):
return self._lf
def __repr__(self):
return "Riemann Zeta function viewed as an L-series"
class LSeriesDelta(LSeriesAbstract):
"""
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeriesDelta; L = LSeriesDelta()
sage: L.anlist(10)
[0, 1, -24, 252, -1472, 4830, -6048, -16744, 84480, -113643, -115920]
sage: list(delta_qexp(11))
[0, 1, -24, 252, -1472, 4830, -6048, -16744, 84480, -113643, -115920]
sage: L.anlist(10^4) == list(delta_qexp(10^4+1))
True
"""
def __init__(self):
LSeriesAbstract.__init__(self, conductor=1, hodge_numbers=[0,1], weight=12, epsilon=1,
poles=[], residues=[], base_field=QQ)
self._T = ZZ['T'].gen()
self._lf = {}
def _local_factor(self, P, prec):
"""
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeriesDelta; L = LSeriesDelta()
sage: L.local_factor(2)
2048*T^2 + 24*T + 1
sage: L._local_factor(11, None) # really this is called
285311670611*T^2 - 534612*T + 1
The representation is reducible modulo 691::
sage: L.local_factor(2).factor_mod(691)
(666) * (T + 387) * (T + 690)
sage: L.local_factor(3).factor_mod(691)
(251) * (T + 234) * (T + 690)
sage: L.local_factor(11).factor_mod(691)
(468) * (T + 471) * (T + 690)
... because of the 691 here::
sage: bernoulli(12)
-691/2730
"""
try:
return self._lf[P]
except KeyError:
pass
self._precompute_local_factors(P+1)
return self._lf[P]
def _precompute_local_factors(self, bound, prec=None):
"""
Precompute local factors up to the given bound.
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeriesDelta; L = LSeriesDelta()
sage: L._lf
{}
sage: L._precompute_local_factors(10)
sage: L._lf
{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}
"""
from sage.modular.all import delta_qexp
T = self._T
T2 = T**2
f = delta_qexp(bound)
for p in prime_range(bound):
if not self._lf.has_key(p):
self._lf[p] = 1 - f[p]*T + (p**11)*T2
def __repr__(self):
return "L-function associated to Ramanujan's Delta (a weight 12 cusp form)"
def _cmp(self, right):
return 0
class LSeriesEllipticCurve(LSeriesAbstract):
def __init__(self, E, prec=53):
"""
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeriesEllipticCurve
sage: L = LSeriesEllipticCurve(EllipticCurve('389a'))
sage: L(2)
0.360092863578881
"""
E = E.global_minimal_model()
self._E = E
K = E.base_field()
d = K.degree()
self._lf = {}
self._T = ZZ['T'].gen()
self._N = E.conductor()
if gcd(norm(self._N), K.discriminant())!=1:
raise NotImplementedError, "Computation of conductor only implemented when the discrimenant of the base field and the norm of the conductor are coprime"
LSeriesAbstract.__init__(self, conductor = norm(self._N) * K.discriminant()**2,
hodge_numbers = [0]*d+[1]*d, weight = 2, epsilon = 1,
poles = [], residues=[], base_field = K, prec=prec)
def elliptic_curve(self):
return self._E
def _cmp(self, right):
return cmp(self.elliptic_curve(), right.elliptic_curve())
def __repr__(self):
return "L-series of %s"%self.elliptic_curve()
def _lf0(self, P):
R = ZZ['T']
T = R.gen()
q = norm(P)
p = prime_below(P)
f = ZZ(q).ord(p)
if P.divides(self._N):
a = self._E.local_data(P).bad_reduction_type()
return 1 - a*(T**f)
else:
a = q + 1 - self._E.reduction(P).count_points()
return 1 - a*(T**f) + q*(T**(2*f))
def _local_factor(self, P, prec):
if not self._lf.has_key(P):
self._lf[P]=self._lf0(P)
return self._lf[P]
def _precompute_local_factors(self, bound, prec=None):
for p in primes(bound):
for q in self._primes_above(p):
if q.norm()>=bound:
continue
if not self._lf.has_key(q):
self._lf[q] = self._lf0(q)
class LSeriesEllipticCurveQQ(LSeriesEllipticCurve):
def __init__(self, E):
E = E.global_minimal_model()
self._E = E
K = E.base_field()
self._N = E.conductor()
self._lf = {}
self._T = ZZ['T'].gen()
LSeriesAbstract.__init__(self, conductor = self._N,
hodge_numbers = [0,1], weight = 2, epsilon = E.root_number(),
poles = [], residues=[], base_field = QQ)
def _lf0(self, p):
a = self._E.ap(p)
T = self._T
if self._N%p == 0:
if self._N%(p*p) == 0:
return T.parent()(1)
else:
return 1 - a*T
else:
return 1 - a*T + p*T*T
def _precompute_local_factors(self, bound, prec=None):
for p in primes(bound):
if not self._lf.has_key(p):
self._lf[p] = self._lf0(p)
class LSeriesEllipticCurveSqrt5(LSeriesEllipticCurve):
def _precompute_local_factors(self, bound, prec=None):
E = self._E
from psage.number_fields.sqrt5.prime import primes_of_bounded_norm, Prime
primes = primes_of_bounded_norm(bound)
from psage.ellcurve.lseries.aplist_sqrt5 import aplist
v = aplist(E, bound)
bad_primes = set([Prime(a.prime()) for a in E.local_data()])
P = ZZ['T']
T = P.gen()
Tp = [T**i for i in range(5)]
if not hasattr(self, '_lf'):
self._lf = {}
for i, P in enumerate(primes):
inertial_deg = 2 if P.is_inert() else 1
a_p = v[i]
if P in bad_primes:
f = 1 - a_p*Tp[inertial_deg]
else:
q = P.norm()
f = 1 - a_p*Tp[inertial_deg] + q*Tp[2*inertial_deg]
self._lf[P] = f
def _local_factor(self, P, prec):
from psage.number_fields.sqrt5.prime import primes_of_bounded_norm, Prime
if not isinstance(P, Prime):
P = Prime(P)
if self._lf.has_key(P):
return self._lf[P]
else:
return LSeriesEllipticCurve._local_factor(self, P.sage_ideal())
def _primes_above(self, p):
"""
Return the primes above p. This function returns a special
optimized prime of the ring of integers of Q(sqrt(5)).
"""
from prime import primes_above
return primes_above(p)
class LSeriesDedekindZeta(LSeriesAbstract):
"""
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeries
sage: K.<a> = NumberField(x^3 - 2)
sage: L = LSeries(K); L
Dedekind Zeta function of Number Field in a with defining polynomial x^3 - 2
sage: L(2)
1.60266326190044
sage: L.residues()
'automatic'
sage: L.residues(prec=53)
[-4.77632833933856]
sage: L.residues(prec=100)
[-4.7763283393385594030639875094]
"""
def __init__(self, K):
if not K.is_absolute():
K = K.absolute_field(names='a')
self._K = K
d = K.degree()
sigma = K.signature()[1]
LSeriesAbstract.__init__(self,
conductor = abs(K.discriminant()),
hodge_numbers = [0]*(d-sigma) + [1]*sigma,
weight = 1,
epsilon = 1,
poles = [1],
residues = 'automatic',
base_field = K,
is_selfdual = True)
self._T = ZZ['T'].gen()
def _cmp(self, right):
return cmp(self.number_field(), right.number_field())
def number_field(self):
return self._K
def __repr__(self):
return "Dedekind Zeta function of %s"%self._K
def _local_factor(self, P, prec):
T = self._T
return 1 - T**P.residue_class_degree()
class LSeriesDirichletCharacter(LSeriesAbstract):
"""
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeries; L = LSeries(DirichletGroup(5).0)
sage: L(3)
0.988191681624057 + 0.0891051883457395*I
"""
def __init__(self, chi):
if not chi.is_primitive():
raise NotImplementedError, "chi must be primitive"
if chi.is_trivial():
raise NotImplementedError, "chi must be nontrivial"
if chi.base_ring().characteristic() != 0:
raise ValueError, "base ring must have characteristic 0"
self._chi = chi
LSeriesAbstract.__init__(self, conductor = chi.conductor(),
hodge_numbers = [1] if chi.is_odd() else [0],
weight = 1,
epsilon = None,
poles = [],
residues = [],
base_field = QQ,
is_selfdual = chi.order() <= 2)
self._T = ZZ['T'].gen()
def _cmp(self, right):
return cmp(self.character(), right.character())
def __repr__(self):
"""
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeries; L = LSeries(DirichletGroup(3).0)
sage: L.__repr__()
'L-series attached to Dirichlet character modulo 3 of conductor 3 mapping 2 |--> -1'
"""
return "L-series attached to %s"%self._chi
def character(self):
return self._chi
def epsilon(self, prec=None):
chi = self._chi
if prec is None:
return (sqrt(-1) * SR(chi.gauss_sum())) / chi.modulus().sqrt()
else:
C = ComplexField(prec)
x = C(chi.modulus()).sqrt() / chi.gauss_sum_numerical(prec=prec)
if chi.is_odd():
x *= C.gen()
return 1/x
def _local_factor(self, P, prec):
a = self._chi(P)
if prec is not None:
a = ComplexField(prec)(a)
return 1 - a*self._T
class LSeriesModularSymbolsAbstract(LSeriesAbstract):
def _cmp(self, right):
return cmp(self.modular_symbols(), right.modular_symbols())
def modular_symbols(self):
return self._M
def __repr__(self):
"""
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeries
sage: f = Newforms(43,2,names='a')[1]; f
q + a1*q^2 - a1*q^3 + (-a1 + 2)*q^5 + O(q^6)
sage: LSeries(f).__repr__()
'L-series of a degree 2 newform of level 43 and weight 2'
"""
return "L-series of a degree %s newform of level %s and weight %s"%(self._M.dimension(), self._M.level(), self._M.weight())
def _precompute_local_factors(self, bound, prec):
primes = [p for p in prime_range(bound) if not self._lf.has_key(p) or self._lf[p][0] < prec]
self._do_precompute(primes, prec)
def _do_precompute(self, primes, prec):
E, v = self._M.compact_system_of_eigenvalues(primes)
if prec == 53:
C = CDF
elif prec is None or prec==oo:
if v.base_ring() == QQ:
C = QQ
else:
C = CDF
else:
C = ComplexField(prec)
phi = v.base_ring().embeddings(C)[self._conjugate]
v = vector(C, [phi(a) for a in v])
aplist = E.change_ring(C) * v
T = C['T'].gen(); T2 = T**2
chi = self._M.character()
k = self.weight()
for i in range(len(primes)):
p = primes[i]
s = chi(p)
if s != 0: s *= p**(k-1)
F = 1 - aplist[i]*T + s*T2
self._lf[p] = (prec, F)
def _local_factor(self, P, prec):
if prec is None: prec = oo
if self._lf.has_key(P) and self._lf[P][0] >= prec:
return self._lf[P][1]
else:
self._do_precompute([P],prec)
return self._lf[P][1]
class LSeriesModularSymbolsNewformGamma0(LSeriesModularSymbolsAbstract):
def _cmp(self, right):
return cmp((self._M, self._conjugate), (right._M, right._conjugate))
def __init__(self, M, conjugate=0, check=True, epsilon=None):
"""
INPUT:
- M -- a simple, new, cuspidal modular symbols space with
sign 1
- conjugate -- (default: 0), integer between 0 and dim(M)-1
- check -- (default: True), if True, checks that M is
simple, new, cuspidal, which can take a very long time,
depending on how M was defined
- epsilon -- (default: None), if not None, should be the sign
in the functional equation, which is -1 or 1. If this is
None, then epsilon is computed by computing the sign of
the main Atkin-Lehner operator on M. If you have a faster
way to determine epsilon, use it.
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeriesModularSymbolsNewformGamma0
sage: M = ModularSymbols(43,sign=1).cuspidal_subspace()[1]; M
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
sage: L0 = LSeriesModularSymbolsNewformGamma0(M,0,check=False,epsilon=1); L0
L-series of a degree 2 newform of level 43 and weight 2
sage: L0.taylor_series()
0.620539857407845 + 0.331674007376949*z - 0.226392184536589*z^2 + 0.0960519649929789*z^3 - 0.00451826124421802*z^4 - 0.0203363026933833*z^5 + O(z^6)
sage: L1 = LSeriesModularSymbolsNewformGamma0(M,1,check=False,epsilon=1); L1
L-series of a degree 2 newform of level 43 and weight 2
sage: L1(1)
0.921328017272472
sage: L1.taylor_series()
0.921328017272472 + 0.492443075089339*z - 0.391019352704047*z^2 + 0.113271812405127*z^3 + 0.0213067052584679*z^4 - 0.0344198080536274*z^5 + O(z^6)
"""
if M.dimension() == 0:
raise ValueError, "modular symbols space must positive dimension"
chi = M.character()
if chi is None or not chi.is_trivial():
raise ValueError, "modular symbols space must have trivial character"
self._M = M
N = M.level()
if check:
if not M.is_simple():
raise ValueError, "modular symbols space must be simple"
if not M.is_new():
raise ValueError, "modular symbols space must be new"
if not M.is_cuspidal():
raise ValueError, "modular symbols space must be cuspidal"
k = M.weight()
if epsilon is None:
w = M.atkin_lehner_operator(N).matrix()
if w not in [-1, 1]:
raise ValueError, "modular symbols space must have constant Atkin-Lehner operator"
epsilon = (-1)**(k/2) * w[0,0]
conjugate = ZZ(conjugate)
if conjugate < 0 or conjugate >= M.dimension():
raise ValueError, "conjugate must a nonnegative integer less than the dimension"
self._conjugate = conjugate
self._lf = {}
LSeriesAbstract.__init__(self,
conductor = N,
hodge_numbers = [0,1],
weight = k,
epsilon = epsilon,
poles = [],
residues = [],
base_field = QQ)
def _is_valid_modsym_space(M):
if not is_ModularSymbolsSpace(M):
raise TypeError, "must be a modular symbols space"
if M.dimension() == 0:
raise ValueError, "modular symbols space must positive dimension"
if M.character() is None:
raise ValueError, "modular symbols space must have associated character"
if not M.is_simple():
raise ValueError, "modular symbols space must be simple"
if not M.is_new():
raise ValueError, "modular symbols space must be new"
if not M.is_cuspidal():
raise ValueError, "modular symbols space must be cuspidal"
class LSeriesModularSymbolsNewformCharacter(LSeriesModularSymbolsAbstract):
def _cmp(self, right):
return cmp((self._M, self._conjugate), (right._M, right._conjugate))
def __init__(self, M, conjugate=0):
_is_valid_modsym_space(M)
chi = M.character()
self._M = M
N = M.level()
epsilon = 'solve'
k = M.weight()
conjugate = ZZ(conjugate)
if conjugate < 0 or conjugate >= M.dimension():
raise ValueError, "conjugate must a nonnegative integer less than the dimension"
self._conjugate = conjugate
LSeriesAbstract.__init__(self,
conductor = N,
hodge_numbers = [0,1],
weight = k,
epsilon = epsilon,
poles = [],
residues = [],
base_field = QQ,
is_selfdual = chi.order() <= 2)
self._lf = {}
class LSeriesModularEllipticCurveSqrt5(LSeriesAbstract):
def __init__(self, E, **kwds):
self._E = E
self._N = E.conductor()
self._T = ZZ['T'].gen()
LSeriesAbstract.__init__(self,
conductor = norm(self._N) * 25,
hodge_numbers = [0,0,1,1],
weight = 2,
epsilon = [1,-1],
poles = [], residues = [],
base_field = E.base_field(),
is_selfdual = True,
**kwds)
def elliptic_curve(self):
return self._E
def _cmp(self, right):
return cmp(self.elliptic_curve(), right.elliptic_curve())
def __repr__(self):
return "L-series attached to %s"%self._E
def _local_factor(self, P, prec):
T = self._T
v = self._N.valuation(P)
if v >= 2:
return T.parent()(1)
elif v >= 1:
d = P.residue_class_degree()
return [1 - T**d, 1 + T**d]
else:
a = self._E.ap(P)
d = P.residue_class_degree()
q = P.norm()
return 1 - a*T**d + q*T**(2*d)
def _new_modsym_space_with_multiplicity(M):
"""
Returns a simple new modular symbols space N and an integer d such
that M is isomorphic to `N^d` as a module over the anemic Hecke
algebra.
INPUT:
- M -- a sign=1 modular simple space for the full Hecke
algebra (including primes dividing the level) that can't be
decomposed further by the Hecke operators. None of the
conditions on M are explicitly checked.
OUTPUT:
- N -- a simple new modular symbols space
- d -- a positive integer
"""
if M.is_new():
return [(M,1)]
raise NotImplementedError
def LSeriesModularSymbolsNewform(M, i=0):
chi = M.character()
if chi is None:
raise NotImplementedError
elif chi.is_trivial():
return LSeriesModularSymbolsNewformGamma0(M, i)
else:
return LSeriesModularSymbolsNewformCharacter(M, i)
class LSeriesModularSymbolsMotive(LSeriesProduct):
"""
The product of L-series attached to the modular symbols space M.
"""
def __init__(self, M):
self._M = M
if not is_ModularSymbolsSpace(M):
raise TypeError, "X must be a modular symbols space or have a modular symbols method"
self._M = M
D = M.decomposition()
for A in D:
_is_valid_modsym_space(A)
F = []
for A in D:
for X in _new_modsym_space_with_multiplicity(A):
N, d = X
chi = N.character()
for i in range(N.dimension()):
F.append( (LSeriesModularSymbolsNewform(N,i), d))
LSeriesProduct.__init__(self, F)
def modular_symbols(self):
return self._M
def __repr__(self):
return "L-series attached to %s"%self._M
class LSeriesModularAbelianVariety(LSeriesProduct):
"""
The product of L-series attached to the modular abelian variety A.
EXAMPLES::
sage: from psage.lseries.eulerprod import LSeries
sage: L = LSeries(J0(54)); L
L-series attached to Abelian variety J0(54) of dimension 4
sage: L.factor()
(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)
sage: L(1)
0.250717238804658
sage: L.taylor_series(prec=20)
0.25072 + 0.59559*z + 0.15099*z^2 - 0.35984*z^3 + 0.056934*z^4 + 0.17184*z^5 + O(z^6)
Independent check of L(1)::
sage: prod(EllipticCurve(lbl).lseries()(1) for lbl in ['54a', '54b', '27a', '27a'])
0.250717238804658
Different check that totally avoids using Dokchitser::
sage: prod(EllipticCurve(lbl).lseries().at1()[0] for lbl in ['54a', '54b', '27a', '27a'])
0.250848605530185
"""
def __init__(self, A):
self._A = A
D = A.decomposition()
F = None
for A in D:
f = Newform(A.newform_label(), names='a')
M = f.modular_symbols(sign=1)
d = ZZ(A.dimension() / M.dimension())
L = LSeriesModularSymbolsMotive(M)**d
if F is None:
F = L
else:
F *= L
if F is None:
raise ValueError, "abelian variety must have positive dimension"
LSeriesProduct.__init__(self, F.factor())
def abelian_variety(self):
return self._A
def __repr__(self):
return "L-series attached to %s"%self._A
class LSeriesTwist(LSeriesAbstract):
"""
Twist of an L-series by a character.
"""
def __init__(self, L, chi, conductor=None, epsilon=None, prec=53):
"""
INPUT:
- `L` -- an L-series
- ``chi`` -- a character of the base field of L
- ``conductor`` -- None, or a list of conductors to try
- ``prec`` -- precision to use when trying conductors, if
conductor is a list
"""
self._L = L
self._chi = chi
if not chi.is_primitive():
raise ValueError, "character must be primitive"
A = ZZ(L.conductor())
B = chi.conductor()
if conductor is None:
if A.gcd(B) != 1:
smallest = ZZ(A)
while smallest.gcd(B) != 1:
smallest = smallest // smallest.gcd(B)
biggest = A * (B**L.degree())
assert biggest % smallest == 0
conductor = [smallest*d for d in divisors(biggest//smallest)]
else:
conductor = A * (B**L.degree())
hodge_numbers = L.hodge_numbers()
weight = L.weight()
if epsilon is None:
if L.epsilon() != 'solve':
if chi.order() <= 2:
if A.gcd(B) == 1:
epsilon = L.epsilon() * chi(-A)
else:
epsilon = [L.epsilon(), -L.epsilon()]
else:
epsilon = 'solve'
else:
epsilon = 'solve'
is_selfdual = L.is_selfdual()
poles = []
residues = 'automatic'
base_field = L.base_field()
LSeriesAbstract.__init__(self, conductor, hodge_numbers, weight, epsilon,
poles, residues, base_field, is_selfdual, prec)
def _local_factor(self, P, prec):
L0 = self._L.local_factor(P, prec)
chi = self._chi
T = L0.parent().gen()
c = chi(P)
if prec is not None:
c = ComplexField(prec)(c)
return L0(c*T)
def __repr__(self):
return "Twist of %s by %s"%(self._L, self._chi)
def _cmp(self, right):
return cmp((self._L, self._chi), (right._L, right._chi))
def untwisted_lseries(self):
return self._L
def twist_character(self):
return self._chi
def LSeries(X, *args, **kwds):
"""
Return the L-series of X, where X can be any of the following:
- elliptic curve over a number field (including QQ)
- Dirichlet character
- cuspidal newform
- new cuspidal modular symbols space -- need not be simple
- string: 'zeta' (Riemann Zeta function), 'delta'
- modular elliptic curve attached to Hilbert modular forms space
For convenience, if L is returned, then L._X is set to X.
EXAMPLES::
The Dedekind Zeta function of a number field::
sage: from psage.lseries.eulerprod import LSeries
sage: K.<a> = NumberField(x^2 + 1)
sage: L = LSeries(K); L
Dedekind Zeta function of Number Field in a with defining polynomial x^2 + 1
sage: L(2)
1.50670300992299
sage: K.zeta_coefficients(100) == L.anlist(100)[1:]
True
sage: L = LSeries(ModularSymbols(43, weight=2,sign=1).cuspidal_subspace().decomposition()[1])
sage: L(1)
0.571720756464112
sage: L.factor()[0][0](1)
0.620539857407845
sage: L.factor()[1][0](1)
0.921328017272472
sage: L._X
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
sage: L = LSeries(ModularSymbols(DirichletGroup(13).0^2, weight=2,sign=1).cuspidal_subspace())
sage: L(1)
0.298115272465799 - 0.0402203326076733*I
sage: from psage.modform.hilbert.sqrt5.hmf import F, HilbertModularForms
sage: D = HilbertModularForms(-13*F.0+5).elliptic_curve_factors()
sage: D
[
Isogeny class of elliptic curves over QQ(sqrt(5)) attached to form number 0 in Hilbert modular forms of dimension 4, level -13*a+5 (of norm 209=11*19) over QQ(sqrt(5)),
Isogeny class of elliptic curves over QQ(sqrt(5)) attached to form number 1 in Hilbert modular forms of dimension 4, level -13*a+5 (of norm 209=11*19) over QQ(sqrt(5)),
Isogeny class of elliptic curves over QQ(sqrt(5)) attached to form number 2 in Hilbert modular forms of dimension 4, level -13*a+5 (of norm 209=11*19) over QQ(sqrt(5))
]
sage: L = LSeries(D[0])
sage: L(RealField(10)(1))
0
sage: L.epsilon()
-1
The L-series of a modular abelian variety with both new and old parts::
sage: L = LSeries(J0(33)); L
L-series attached to Abelian variety J0(33) of dimension 3
sage: L.factor()
(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)
sage: L.local_factor(2, prec=oo)
8*T^6 + 12*T^5 + 12*T^4 + 8*T^3 + 6*T^2 + 3*T + 1
sage: L(1)
0.0481553138900504
We check the above computation of L(1) via independent methods (and implementations)::
sage: prod(EllipticCurve(lbl).lseries().at1()[0] for lbl in ['11a', '11a', '33a'])
0.0481135342926321
sage: prod(EllipticCurve(lbl).lseries()(1) for lbl in ['11a', '11a', '33a'])
0.0481553138900504
A nonsimple new modular symbols space of level 43::
sage: L = LSeries(ModularSymbols(43,sign=1).cuspidal_subspace())
sage: L
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
sage: L(1)
0
sage: L.taylor_series()
0.196399786632435*z + 0.314922741074845*z^2 - 0.0797083673829092*z^3 - 0.161630566287135*z^4 + 0.123939472976207*z^5 + O(z^6)
sage: L.factor()
(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)
sage: L.analytic_rank()
1
sage: D = ModularSymbols(43,sign=1).cuspidal_subspace().decomposition()
sage: L0 = LSeries(D[0]); L1 = LSeries(D[1])
sage: L0.taylor_series() * L1.taylor_series()
0.196399786632435*z + 0.314922741074845*z^2 - 0.0797083673829091*z^3 - 0.161630566287135*z^4 + 0.123939472976207*z^5 + O(z^6)
sage: L0.factor()
L-series of a degree 1 newform of level 43 and weight 2
sage: L1.factor()
(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)
"""
L = _lseries(X, *args, **kwds)
L._X = X
return L
def _lseries(X, *args, **kwds):
"""
Helper function used by LSeries function.
"""
if is_EllipticCurve(X):
K = X.base_ring()
if is_RationalField(K):
return LSeriesEllipticCurveQQ(X, *args, **kwds)
elif list(K.defining_polynomial()) == [-1,-1,1]:
return LSeriesEllipticCurveSqrt5(X, *args, **kwds)
else:
return LSeriesEllipticCurve(X)
if is_DirichletCharacter(X):
if X.is_trivial() and X.is_primitive():
return LSeriesZeta(*args, **kwds)
else:
return LSeriesDirichletCharacter(X, *args, **kwds)
if is_NumberField(X):
return LSeriesDedekindZeta(X, *args, **kwds)
if isinstance(X, sage.modular.modform.element.Newform):
return LSeriesModularSymbolsNewform(X.modular_symbols(sign=1), *args, **kwds)
if is_ModularSymbolsSpace(X):
if X.sign() != 1:
raise NotImplementedError
return LSeriesModularSymbolsMotive(X, *args, **kwds)
if isinstance(X, str):
y = X.lower()
if y == 'zeta':
return LSeriesZeta(*args, **kwds)
elif y == 'delta':
return LSeriesDelta(*args, **kwds)
else:
raise ValueError, 'unknown L-series "%s"'%y
if isinstance(X, psage.modform.hilbert.sqrt5.hmf.EllipticCurveFactor):
return LSeriesModularEllipticCurveSqrt5(X, *args, **kwds)
if is_ModularAbelianVariety(X):
return LSeriesModularAbelianVariety(X, *args, **kwds)
raise NotImplementedError
