"""
Numerical computation of newforms
"""
from sage.structure.sage_object import SageObject
from sage.structure.sequence import Sequence
from sage.modular.modsym.all import ModularSymbols
from sage.modular.arithgroup.all import Gamma0
from sage.modules.all import vector
from sage.misc.misc import verbose
from sage.rings.all import CDF, Integer, QQ, next_prime, prime_range
from sage.misc.prandom import randint
from sage.matrix.constructor import matrix
scipy=None
class NumericalEigenforms(SageObject):
"""
numerical_eigenforms(group, weight=2, eps=1e-20, delta=1e-2, tp=[2,3,5])
INPUT:
- ``group`` - a congruence subgroup of a Dirichlet character of
order 1 or 2
- ``weight`` - an integer >= 2
- ``eps`` - a small float; abs( ) < eps is what "equal to zero" is
interpreted as for floating point numbers.
- ``delta`` - a small-ish float; eigenvalues are considered distinct
if their difference has absolute value at least delta
- ``tp`` - use the Hecke operators T_p for p in tp when searching
for a random Hecke operator with distinct Hecke eigenvalues.
OUTPUT:
A numerical eigenforms object, with the following useful methods:
- :meth:`ap` - return all eigenvalues of $T_p$
- :meth:`eigenvalues` - list of eigenvalues corresponding
to the given list of primes, e.g.,::
[[eigenvalues of T_2],
[eigenvalues of T_3],
[eigenvalues of T_5], ...]
- :meth:`systems_of_eigenvalues` - a list of the systems of
eigenvalues of eigenforms such that the chosen random linear
combination of Hecke operators has multiplicity 1 eigenvalues.
EXAMPLES::
sage: n = numerical_eigenforms(23)
sage: n == loads(dumps(n))
True
sage: n.ap(2)
[3.0, 0.61803398875, -1.61803398875]
sage: n.systems_of_eigenvalues(7)
[
[-1.61803398875, 2.2360679775, -3.2360679775],
[0.61803398875, -2.2360679775, 1.2360679775],
[3.0, 4.0, 6.0]
]
sage: n.systems_of_abs(7)
[
[0.6180339887..., 2.236067977..., 1.236067977...],
[1.6180339887..., 2.236067977..., 3.236067977...],
[3.0, 4.0, 6.0]
]
sage: n.eigenvalues([2,3,5])
[[3.0, 0.61803398875, -1.61803398875],
[4.0, -2.2360679775, 2.2360679775],
[6.0, 1.2360679775, -3.2360679775]]
"""
def __init__(self, group, weight=2, eps=1e-20,
delta=1e-2, tp=[2,3,5]):
"""
Create a new space of numerical eigenforms.
EXAMPLES::
sage: numerical_eigenforms(61) # indirect doctest
Numerical Hecke eigenvalues for Congruence Subgroup Gamma0(61) of weight 2
"""
if isinstance(group, (int, long, Integer)):
group = Gamma0(Integer(group))
self._group = group
self._weight = Integer(weight)
self._tp = tp
if self._weight < 2:
raise ValueError, "weight must be at least 2"
self._eps = eps
self._delta = delta
def __cmp__(self, other):
"""
Compare two spaces of numerical eigenforms. Currently
returns 0 if they come from the same space of modular
symbols, and -1 otherwise.
EXAMPLES::
sage: n = numerical_eigenforms(23)
sage: n.__cmp__(loads(dumps(n)))
0
"""
if not isinstance( other, NumericalEigenforms ):
raise ValueError, "%s is not a space of numerical eigenforms"%other
if self.modular_symbols() == other.modular_symbols():
return 0
else:
return -1
def level(self):
"""
Return the level of this set of modular eigenforms.
EXAMPLES::
sage: n = numerical_eigenforms(61) ; n.level()
61
"""
return self._group.level()
def weight(self):
"""
Return the weight of this set of modular eigenforms.
EXAMPLES::
sage: n = numerical_eigenforms(61) ; n.weight()
2
"""
return self._weight
def _repr_(self):
"""
Print string representation of self.
EXAMPLES::
sage: n = numerical_eigenforms(61) ; n
Numerical Hecke eigenvalues for Congruence Subgroup Gamma0(61) of weight 2
sage: n._repr_()
'Numerical Hecke eigenvalues for Congruence Subgroup Gamma0(61) of weight 2'
"""
return "Numerical Hecke eigenvalues for %s of weight %s"%(
self._group, self._weight)
def modular_symbols(self):
"""
Return the space of modular symbols used for computing this
set of modular eigenforms.
EXAMPLES::
sage: n = numerical_eigenforms(61) ; n.modular_symbols()
Modular Symbols space of dimension 5 for Gamma_0(61) of weight 2 with sign 1 over Rational Field
"""
try:
return self.__modular_symbols
except AttributeError:
M = ModularSymbols(self._group,
self._weight, sign=1)
if M.base_ring() != QQ:
raise ValueError, "modular forms space must be defined over QQ"
self.__modular_symbols = M
return M
def _eigenvectors(self):
r"""
Find numerical approximations to simultaneous eigenvectors in
self.modular_symbols() for all T_p in self._tp.
EXAMPLES::
sage: n = numerical_eigenforms(61)
sage: n._eigenvectors() # random order
[ 1.0 0.289473640239 0.176788851952 0.336707726757 2.4182243084e-16]
[ 0 -0.0702748344418 0.491416161212 0.155925712173 0.707106781187]
[ 0 0.413171180356 0.141163094698 0.0923242547901 0.707106781187]
[ 0 0.826342360711 0.282326189397 0.18464850958 6.79812569682e-16]
[ 0 0.2402380858 0.792225196393 0.905370774276 4.70805946682e-16]
TESTS:
This tests if this routine selects only eigenvectors with
multiplicity one. Two of the eigenvalues are
(roughly) -92.21 and -90.30 so if we set ``eps = 2.0``
then they should compare as equal, causing both eigenvectors
to be absent from the matrix returned. The remaining eigenvalues
(ostensibly unique) are visible in the test, which should be
indepedent of which eigenvectors are returned, but it does presume
an ordering of these eigenvectors for the test to succeed.
This exercises a correction in Trac 8018. ::
sage: n = numerical_eigenforms(61, eps=2.0)
sage: evectors = n._eigenvectors()
sage: evalues = diagonal_matrix(CDF, [-283.0, 108.522012456, 142.0])
sage: diff = n._hecke_matrix*evectors - evectors*evalues
sage: sum([abs(diff[i,j]) for i in range(5) for j in range(3)]) < 1.0e-9
True
"""
try:
return self.__eigenvectors
except AttributeError:
pass
verbose('Finding eigenvector basis')
M = self.modular_symbols()
N = self.level()
tp = self._tp
p = tp[0]
t = M.T(p).matrix()
for p in tp[1:]:
t += randint(-50,50)*M.T(p).matrix()
self._hecke_matrix = t
global scipy
if scipy is None:
import scipy
import scipy.linalg
evals,eig = scipy.linalg.eig(self._hecke_matrix.numpy(), right=True, left=False)
B = matrix(eig)
v = [CDF(evals[i]) for i in range(len(evals))]
eps = self._eps
w = []
for i in range(len(v)):
e = v[i]
uniq = True
for j in range(len(v)):
if uniq and i != j and abs(e-v[j]) < eps:
uniq = False
if uniq:
w.append(i)
self.__eigenvectors = B.matrix_from_columns(w)
return self.__eigenvectors
def _easy_vector(self):
"""
Return a very sparse vector v such that v times the eigenvector matrix
has all entries nonzero.
ALGORITHM:
1. Choose row with the most nonzero entries. (put 1 there)
2. Consider submatrix of columns corresponding to zero entries
in row chosen in 1.
3. Find row of submatrix with most nonzero entries, and add
appropriate multiple. Repeat.
EXAMPLES::
sage: n = numerical_eigenforms(37)
sage: n._easy_vector() # slightly random output
(1.0, 1.0, 0)
sage: n = numerical_eigenforms(43)
sage: n._easy_vector() # slightly random output
(1.0, 0, 1.0, 0)
sage: n = numerical_eigenforms(125)
sage: n._easy_vector() # slightly random output
(0, 0, 0, 1.0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
"""
try:
return self.__easy_vector
except AttributeError:
pass
E = self._eigenvectors()
delta = self._delta
x = (CDF**E.nrows()).zero_vector()
if E.nrows() == 0:
return x
def best_row(M):
"""
Find the best row among rows of M, i.e. the row
with the most entries supported outside [-delta, delta].
EXAMPLES::
sage: numerical_eigenforms(61)._easy_vector() # indirect doctest
(1.0, 1.0, 0.0, 0.0, 0.0)
"""
R = M.rows()
v = [len(support(r, delta)) for r in R]
m = max(v)
i = v.index(m)
return i, R[i]
i, e = best_row(E)
x[i] = 1
while True:
s = set(support(e, delta))
zp = [i for i in range(e.degree()) if not i in s]
if len(zp) == 0:
break
C = E.matrix_from_columns(zp)
i, f = best_row(C)
x[i] += 1
e = x*E
self.__easy_vector = x
return x
def _eigendata(self):
"""
Return all eigendata for self._easy_vector().
EXAMPLES::
sage: numerical_eigenforms(61)._eigendata() # random order
((1.0, 0.668205013164, 0.219198805797, 0.49263343893, 0.707106781187), (1.0, 1.49654668896, 4.5620686498, 2.02990686579, 1.41421356237), [0, 1], (1.0, 1.0))
"""
try:
return self.__eigendata
except AttributeError:
pass
x = self._easy_vector()
B = self._eigenvectors()
def phi(y):
"""
Take coefficients and a basis, and return that
linear combination of basis vectors.
EXAMPLES::
sage: n = numerical_eigenforms(61) # indirect doctest
sage: n._eigendata() # random order
((1.0, 0.668205013164, 0.219198805797, 0.49263343893, 0.707106781187), (1.0, 1.49654668896, 4.5620686498, 2.02990686579, 1.41421356237), [0, 1], (1.0, 1.0))
"""
return y.element() * B
phi_x = phi(x)
V = phi_x.parent()
phi_x_inv = V([a**(-1) for a in phi_x])
eps = self._eps
nzp = support(x, eps)
x_nzp = vector(CDF, x.list_from_positions(nzp))
self.__eigendata = (phi_x, phi_x_inv, nzp, x_nzp)
return self.__eigendata
def ap(self, p):
"""
Return a list of the eigenvalues of the Hecke operator `T_p`
on all the computed eigenforms. The eigenvalues match up
between one prime and the next.
INPUT:
- ``p`` - integer, a prime number
OUTPUT:
- ``list`` - a list of double precision complex numbers
EXAMPLES::
sage: n = numerical_eigenforms(11,4)
sage: n.ap(2) # random order
[9.0, 9.0, 2.73205080757, -0.732050807569]
sage: n.ap(3) # random order
[28.0, 28.0, -7.92820323028, 5.92820323028]
sage: m = n.modular_symbols()
sage: x = polygen(QQ, 'x')
sage: m.T(2).charpoly(x).factor()
(x - 9)^2 * (x^2 - 2*x - 2)
sage: m.T(3).charpoly(x).factor()
(x - 28)^2 * (x^2 + 2*x - 47)
"""
p = Integer(p)
if not p.is_prime():
raise ValueError, "p must be a prime"
try:
return self._ap[p]
except AttributeError:
self._ap = {}
except KeyError:
pass
a = Sequence(self.eigenvalues([p])[0], immutable=True)
self._ap[p] = a
return a
def eigenvalues(self, primes):
"""
Return the eigenvalues of the Hecke operators corresponding
to the primes in the input list of primes. The eigenvalues
match up between one prime and the next.
INPUT:
- ``primes`` - a list of primes
OUTPUT:
list of lists of eigenvalues.
EXAMPLES::
sage: n = numerical_eigenforms(1,12)
sage: n.eigenvalues([3,5,13])
[[177148.0, 252.0], [48828126.0, 4830.0], [1.79216039404e+12, -577737.999...]]
"""
primes = [Integer(p) for p in primes]
for p in primes:
if not p.is_prime():
raise ValueError, 'each element of primes must be prime.'
phi_x, phi_x_inv, nzp, x_nzp = self._eigendata()
B = self._eigenvectors()
def phi(y):
"""
Take coefficients and a basis, and return that
linear combination of basis vectors.
EXAMPLES::
sage: n = numerical_eigenforms(1,12) # indirect doctest
sage: n.eigenvalues([3,5,13])
[[177148.0, 252.0], [48828126.0, 4830.0], [1.79216039404e+12, -577737.999...]]
"""
return y.element() * B
ans = []
m = self.modular_symbols().ambient_module()
for p in primes:
t = m._compute_hecke_matrix_prime(p, nzp)
w = phi(x_nzp*t)
ans.append([w[i]*phi_x_inv[i] for i in range(w.degree())])
return ans
def systems_of_eigenvalues(self, bound):
"""
Return all systems of eigenvalues for self for primes
up to bound.
EXAMPLES::
sage: numerical_eigenforms(61).systems_of_eigenvalues(10)
[
[-1.48119430409..., 0.806063433525..., 3.15632517466..., 0.675130870567...],
[-1.0..., -2.0..., -3.0..., 1.0...],
[0.311107817466..., 2.90321192591..., -2.52542756084..., -3.21431974338...],
[2.17008648663..., -1.70927535944..., -1.63089761382..., -0.460811127189...],
[3.0, 4.0, 6.0, 8.0]
]
"""
P = prime_range(bound)
e = self.eigenvalues(P)
v = Sequence([], cr=True)
if len(e) == 0:
return v
for i in range(len(e[0])):
v.append([e[j][i] for j in range(len(e))])
v.sort()
v.set_immutable()
return v
def systems_of_abs(self, bound):
"""
Return the absolute values of all systems of eigenvalues for
self for primes up to bound.
EXAMPLES::
sage: numerical_eigenforms(61).systems_of_abs(10)
[
[0.311107817466, 2.90321192591, 2.52542756084, 3.21431974338],
[1.0, 2.0, 3.0, 1.0],
[1.48119430409, 0.806063433525, 3.15632517466, 0.675130870567],
[2.17008648663, 1.70927535944, 1.63089761382, 0.460811127189],
[3.0, 4.0, 6.0, 8.0]
]
"""
P = prime_range(bound)
e = self.eigenvalues(P)
v = Sequence([], cr=True)
if len(e) == 0:
return v
for i in range(len(e[0])):
v.append([abs(e[j][i]) for j in range(len(e))])
v.sort()
v.set_immutable()
return v
def support(v, eps):
"""
Given a vector `v` and a threshold eps, return all
indices where `|v|` is larger than eps.
EXAMPLES::
sage: sage.modular.modform.numerical.support( numerical_eigenforms(61)._easy_vector(), 1.0 )
[]
sage: sage.modular.modform.numerical.support( numerical_eigenforms(61)._easy_vector(), 0.5 )
[0, 1]
"""
return [i for i in range(v.degree()) if abs(v[i]) > eps]
