Hot-keys on this page
r m x p toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
""" Heilbronn matrix computation """
#***************************************************************************** # Copyright (C) 2004 William Stein <wstein@gmail.com> # # Distributed under the terms of the GNU General Public License (GPL) # # This code is distributed in the hope that it will be useful, # but WITHOUT ANY WARRANTY; without even the implied warranty of # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU # General Public License for more details. # # The full text of the GPL is available at: # # http://www.gnu.org/licenses/ #***************************************************************************** from __future__ import absolute_import
from cysignals.memory cimport check_allocarray, sig_malloc, sig_free from cysignals.signals cimport sig_on, sig_off
import sage.arith.all
import sage.misc.misc
from sage.libs.gmp.mpz cimport * from sage.libs.gmp.mpq cimport * from sage.libs.flint.fmpz cimport fmpz_add from sage.libs.flint.fmpz_poly cimport * from sage.libs.flint.fmpq_mat cimport fmpq_mat_entry_num
cdef extern from "<math.h>": float roundf(float x)
cimport sage.modular.modsym.p1list as p1list from . import p1list cdef p1list.export export export = p1list.export()
from .apply cimport Apply cdef Apply PolyApply= Apply()
from sage.rings.integer cimport Integer from sage.matrix.matrix_rational_dense cimport Matrix_rational_dense from sage.matrix.matrix_cyclo_dense cimport Matrix_cyclo_dense
ctypedef long long llong
cdef int llong_prod_mod(int a, int b, int N): cdef int c c = <int> ( ((<llong> a) * (<llong> b)) % (<llong> N) ) if c < 0: c = c + N return c
cdef struct list: int *v int i # how many positions of list are filled int n # how much memory has been allocated
cdef int* expand(int *v, int n, int new_length) except NULL: cdef int *w cdef int i
cdef int list_append(list* L, int a) except -1: cdef int j
cdef int list_append4(list* L, int a, int b, int c, int d) except -1:
cdef void list_clear(list L):
cdef void list_init(list* L):
cdef class Heilbronn: cdef int length cdef list list
def __dealloc__(self):
def _initialize_list(self): """ Initialize the list of matrices corresponding to self. (This function is automatically called during initialization.)
.. note::
This function must be overridden by all derived classes!
EXAMPLES::
sage: H = sage.modular.modsym.heilbronn.Heilbronn() sage: H._initialize_list() Traceback (most recent call last): ... NotImplementedError """
def __getitem__(self, int n): """ Return the nth matrix in self.
EXAMPLES::
sage: H = HeilbronnCremona(11) sage: H[17] [-1, 0, 1, -11] sage: H[98234] Traceback (most recent call last): ... IndexError """
def __len__(self): """ Return the number of matrices in self.
EXAMPLES::
sage: HeilbronnCremona(2).__len__() 4 """
def to_list(self): """ Return the list of Heilbronn matrices corresponding to self. Each matrix is given as a list of four ints.
EXAMPLES::
sage: H = HeilbronnCremona(2); H The Cremona-Heilbronn matrices of determinant 2 sage: H.to_list() [[1, 0, 0, 2], [2, 0, 0, 1], [2, 1, 0, 1], [1, 0, 1, 2]] """ cdef int i
cdef apply_only(self, int u, int v, int N, int* a, int* b): """ INPUT:
- ``u, v, N`` - integers
- ``a, b`` - preallocated int arrays of the length self.
OUTPUT: sets the entries of a,b
EXAMPLES::
sage: M = ModularSymbols(33,2,1) # indirect test sage: sage.modular.modsym.heilbronn.hecke_images_gamma0_weight2(1,0,33,[2,3],M.manin_gens_to_basis()) [ 3 0 1 0 -1 1] [ 3 2 2 0 -2 2] sage: z = M((1,0)) sage: [M.T(n)(z).element() for n in [2,2]] [(3, 0, 1, 0, -1, 1), (3, 0, 1, 0, -1, 1)] """ cdef Py_ssize_t i elif N < 46340: # use ints but reduce mod N so can add two for i in range(self.length): a[i] = (u * self.list.v[4*i])%N + (v * self.list.v[4*i+2])%N b[i] = (u * self.list.v[4*i+1])%N + (v * self.list.v[4*i+3])%N else: for i in range(self.length): a[i] = llong_prod_mod(u,self.list.v[4*i],N) + llong_prod_mod(v,self.list.v[4*i+2], N) b[i] = llong_prod_mod(u,self.list.v[4*i+1],N) + llong_prod_mod(v,self.list.v[4*i+3], N)
cdef apply_to_polypart(self, fmpz_poly_t* ans, int i, int k): """ INPUT:
- ``ans`` - fmpz_poly_t\*; pre-allocated an initialized array of self.length fmpz_poly_t's - ``i`` - integer - ``k`` - integer
OUTPUT: sets entries of ans """ self.list.v[4*j], self.list.v[4*j+1], self.list.v[4*j+2], self.list.v[4*j+3])
def apply(self, int u, int v, int N): """ Return a list of pairs ((c,d),m), which is obtained as follows: 1) Compute the images (a,b) of the vector (u,v) (mod N) acted on by each of the HeilbronnCremona matrices in self. 2) Reduce each (a,b) to canonical form (c,d) using p1normalize 3) Sort. 4) Create the list ((c,d),m), where m is the number of times that (c,d) appears in the list created in steps 1-3 above. Note that the pairs ((c,d),m) are sorted lexicographically by (c,d).
INPUT:
- ``u, v, N`` - integers
OUTPUT: list
EXAMPLES::
sage: H = sage.modular.modsym.heilbronn.HeilbronnCremona(2); H The Cremona-Heilbronn matrices of determinant 2 sage: H.apply(1,2,7) [((1, 5), 1), ((1, 6), 1), ((1, 1), 1), ((1, 4), 1)] """ cdef int i, a, b, c, d, s cdef object X else: elif N < 46340: # use ints but reduce mod N so can add two for i in range(self.length): a = (u * self.list.v[4*i])%N + (v * self.list.v[4*i+2])%N b = (u * self.list.v[4*i+1])%N + (v * self.list.v[4*i+3])%N export.c_p1_normalize_int(N, a, b, &c, &d, &s, 0) X = (c,d) if X in M: M[X] = M[X] + 1 else: M[X] = 1 else: for i in range(self.length): a = llong_prod_mod(u,self.list.v[4*i],N) + llong_prod_mod(v,self.list.v[4*i+2], N) b = llong_prod_mod(u,self.list.v[4*i+1],N) + llong_prod_mod(v,self.list.v[4*i+3], N) export.c_p1_normalize_llong(N, a, b, &c, &d, &s, 0) X = (c,d) if X in M: M[X] = M[X] + 1 else: M[X] = 1
cdef class HeilbronnCremona(Heilbronn): cdef public int p
def __init__(self, int p): """ Create the list of Heilbronn-Cremona matrices of determinant p.
EXAMPLES::
sage: H = HeilbronnCremona(3) ; H The Cremona-Heilbronn matrices of determinant 3 sage: H.to_list() [[1, 0, 0, 3], [3, 1, 0, 1], [1, 0, 1, 3], [3, 0, 0, 1], [3, -1, 0, 1], [-1, 0, 1, -3]] """ raise ValueError("p must be >= 2 and prime")
def __repr__(self): """ Return the string representation of self.
EXAMPLES::
sage: HeilbronnCremona(691).__repr__() 'The Cremona-Heilbronn matrices of determinant 691' """
def _initialize_list(self): """ Initialize the list of matrices corresponding to self. (This function is automatically called during initialization.)
EXAMPLES::
sage: H = HeilbronnCremona.__new__(HeilbronnCremona) sage: H.p = 5 sage: H The Cremona-Heilbronn matrices of determinant 5 sage: H.to_list() [] sage: H._initialize_list() sage: H.to_list() [[1, 0, 0, 5], [5, 2, 0, 1], [2, 1, 1, 3], [1, 0, 3, 5], [5, 1, 0, 1], [1, 0, 1, 5], [5, 0, 0, 1], [5, -1, 0, 1], [-1, 0, 1, -5], [5, -2, 0, 1], [-2, 1, 1, -3], [1, 0, -3, 5]] """ cdef int r, x1, x2, y1, y2, a, b, c, x3, y3, q, n, p cdef list *L
# When p==2, then Heilbronn matrices are # [[1,0,0,2], [2,0,0,1], [2,1,0,1], [1,0,1,2]] # which are not given by the algorithm below.
# NOTE: In C, -p/2 means "round toward 0", but in Python it # means "round down."
cdef class HeilbronnMerel(Heilbronn): cdef public int n
def __init__(self, int n): r""" Initialize the list of Merel-Heilbronn matrices of determinant `n`.
EXAMPLES::
sage: H = HeilbronnMerel(3) ; H The Merel-Heilbronn matrices of determinant 3 sage: H.to_list() [[1, 0, 0, 3], [1, 0, 1, 3], [1, 0, 2, 3], [2, 1, 1, 2], [3, 0, 0, 1], [3, 1, 0, 1], [3, 2, 0, 1]] """ raise ValueError("n (=%s) must be >= 1" % n)
def __repr__(self): """ Return the string representation of self.
EXAMPLES::
sage: HeilbronnMerel(8).__repr__() 'The Merel-Heilbronn matrices of determinant 8' """
def _initialize_list(self): """ Initialize the list of matrices corresponding to self. (This function is automatically called during initialization.)
EXAMPLES::
sage: H = HeilbronnMerel.__new__(HeilbronnMerel) sage: H.n = 5 sage: H The Merel-Heilbronn matrices of determinant 5 sage: H.to_list() [] sage: H._initialize_list() sage: H.to_list() [[1, 0, 0, 5], [1, 0, 1, 5], [1, 0, 2, 5], [1, 0, 3, 5], [1, 0, 4, 5], [2, 1, 1, 3], [2, 1, 3, 4], [3, 1, 1, 2], [3, 2, 2, 3], [4, 3, 1, 2], [5, 0, 0, 1], [5, 1, 0, 1], [5, 2, 0, 1], [5, 3, 0, 1], [5, 4, 0, 1]] """ cdef int a, q, d, b, c, n cdef llong bc cdef list *L
## We have ad-bc=n so c=0 and ad=n, or b=(ad-n)/c ## Must have ad - n >= 0, so d must be >= Ceiling(n/a). ## Divisor c of bc must satisfy Floor(bc/c) lt a and c lt d. ## c ge (bc div a + 1) <=> Floor(bc/c) lt a (for integers) ## c le d - 1 <=> c lt d
############################################################################ # Fast application of Heilbronn operators to computation of # systems of Hecke eigenvalues. # GAMMA0 trivial character weight 2 case ############################################################################
def hecke_images_gamma0_weight2(int u, int v, int N, indices, R): """ INPUT:
- ``u, v, N`` - integers so that gcd(u,v,N) = 1 - ``indices`` - a list of positive integers - ``R`` - matrix over QQ that writes each elements of P1 = P1List(N) in terms of a subset of P1.
OUTPUT: a dense matrix whose columns are the images T_n(x) for n in indices and x the Manin symbol (u,v), expressed in terms of the basis.
EXAMPLES::
sage: M = ModularSymbols(23,2,1) sage: A = sage.modular.modsym.heilbronn.hecke_images_gamma0_weight2(1,0,23,[1..6],M.manin_gens_to_basis()) sage: A [ 1 0 0] [ 3 0 -1] [ 4 -2 -1] [ 7 -2 -2] [ 6 0 -2] [12 -2 -4] sage: z = M((1,0)) sage: [M.T(n)(z).element() for n in [1..6]] [(1, 0, 0), (3, 0, -1), (4, -2, -1), (7, -2, -2), (6, 0, -2), (12, -2, -4)]
TESTS::
sage: M = ModularSymbols(389,2,1,GF(7)) sage: C = M.cuspidal_subspace() sage: N = C.new_subspace() sage: D = N.decomposition() sage: D[1].q_eigenform(10, 'a') # indirect doctest q + 4*q^2 + 2*q^3 + 6*q^5 + q^6 + 5*q^7 + 6*q^8 + q^9 + O(q^10)
"""
# Create a zero dense matrix over QQ with len(indices) rows # and #P^1(N) columns. cdef Matrix_rational_dense T
cdef Py_ssize_t i, j cdef int *a cdef int *b cdef int k
cdef Heilbronn H
# List the Heilbronn matrices of determinant n defined by Cremona or Merel
# Allocate memory to hold images of (u,v) under all Heilbronn matrices
# Compute images of (u,v) under all Heilbronn matrices
# Compute the indexes of these images. # We just store them in the array a for simplicity. # Compute index of the symbol a[j], b[j] in the standard list.
# Now fill in row i of the matrix T. # The following line is just a dangerous direct way to do: T[i,k] += 1
# Free a and b
# Return the product T * R, whose rows are the image of (u,v) under # the Hecke operators T_n for n in indices. else: ans = T * R sage.misc.misc.verbose("did reduction using dense multiplication", t, level=1, caller_name='hecke_images_gamma0_weight2')
############################################################################ # Fast application of Heilbronn operators to computation of # systems of Hecke eigenvalues. # Nontrivial character but weight 2. ############################################################################
def hecke_images_nonquad_character_weight2(int u, int v, int N, indices, chi, R): """ Return images of the Hecke operators `T_n` for `n` in the list indices, where chi must be a quadratic Dirichlet character with values in QQ.
R is assumed to be the relation matrix of a weight modular symbols space over QQ with character chi.
INPUT:
- ``u, v, N`` - integers so that gcd(u,v,N) = 1 - ``indices`` - a list of positive integers - ``chi`` - a Dirichlet character that takes values in a nontrivial extension of QQ. - ``R`` - matrix over QQ that writes each elements of P1 = P1List(N) in terms of a subset of P1.
OUTPUT: a dense matrix with entries in the field QQ(chi) (the values of chi) whose columns are the images T_n(x) for n in indices and x the Manin symbol (u,v), expressed in terms of the basis.
EXAMPLES::
sage: chi = DirichletGroup(13).0^2 sage: M = ModularSymbols(chi) sage: eps = M.character() sage: R = M.manin_gens_to_basis() sage: sage.modular.modsym.heilbronn.hecke_images_nonquad_character_weight2(1,0,13,[1,2,6],eps,R) [ 1 0 0 0] [ zeta6 + 2 0 0 -1] [ 7 -2*zeta6 + 1 -zeta6 - 1 -2*zeta6] sage: x = M((1,0)); x.element() (1, 0, 0, 0) sage: M.T(2)(x).element() (zeta6 + 2, 0, 0, -1) sage: M.T(6)(x).element() (7, -2*zeta6 + 1, -zeta6 - 1, -2*zeta6) """
raise TypeError("character must not be trivial or quadratic")
R = R.change_ring(K)
# Create a zero dense matrix over K with len(indices) rows # and #P^1(N) columns. cdef Matrix_cyclo_dense T
cdef Py_ssize_t i, j cdef int *a cdef int *b cdef int k, scalar
cdef Heilbronn H
# Make a matrix over the rational numbers each of whose columns # are the values of the character chi. cdef Matrix_rational_dense chi_vals
# Allocate memory to hold images of (u,v) under all Heilbronn matrices
# Compute images of (u,v) under all Heilbronn matrices
# Compute index of the symbol a[j], b[j] in the standard list. # Now fill in row i of the matrix T. # The following line is just a dangerous direct way to do: T[i,k] += chi(scalar) # T[i,k] += chi(scalar) # This code makes assumptions about the internal structure # of matrices over cyclotomic fields. It's nasty, but it # is exactly what is needed to get a solid 100 or more # times speedup. # Note that the next line totally dominates the runtime of this whole function.
# Free a and b
def hecke_images_quad_character_weight2(int u, int v, int N, indices, chi, R): """ INPUT:
- ``u, v, N`` - integers so that gcd(u,v,N) = 1 - ``indices`` - a list of positive integers - ``chi`` - a Dirichlet character that takes values in QQ - ``R`` - matrix over QQ(chi) that writes each elements of P1 = P1List(N) in terms of a subset of P1.
OUTPUT: a dense matrix with entries in the rational field QQ (the values of chi) whose columns are the images T_n(x) for n in indices and x the Manin symbol (u,v), expressed in terms of the basis.
EXAMPLES::
sage: chi = DirichletGroup(29,QQ).0 sage: M = ModularSymbols(chi) sage: R = M.manin_gens_to_basis() sage: sage.modular.modsym.heilbronn.hecke_images_quad_character_weight2(2,1,29,[1,3,4],chi,R) [ 0 0 0 0 0 -1] [ 0 1 0 1 1 1] [ 0 -2 0 2 -2 -1] sage: x = M((2,1)) ; x.element() (0, 0, 0, 0, 0, -1) sage: M.T(3)(x).element() (0, 1, 0, 1, 1, 1) sage: M.T(4)(x).element() (0, -2, 0, 2, -2, -1) """ raise TypeError("character must takes values in QQ")
# Create a zero dense matrix over QQ with len(indices) rows # and #P^1(N) columns. cdef Matrix_rational_dense T
R = R.change_ring(QQ)
cdef Py_ssize_t i, j cdef int *a cdef int *b cdef int k, scalar cdef Heilbronn H
# Make a matrix over the rational numbers each of whose columns # are the values of the character chi.
# This is just T[i,k] += chi(scalar)
############################################################################ # Fast application of Heilbronn operators to computation of # systems of Hecke eigenvalues. # Trivial character and weight > 2. ############################################################################
def hecke_images_gamma0_weight_k(int u, int v, int i, int N, int k, indices, R): """ INPUT:
- ``u, v, N`` - integers so that gcd(u,v,N) = 1 - ``i`` - integer with 0 = i = k-2 - ``k`` - weight - ``indices`` - a list of positive integers - ``R`` - matrix over QQ that writes each elements of P1 = P1List(N) in terms of a subset of P1.
OUTPUT: a dense matrix with rational entries whose columns are the images T_n(x) for n in indices and x the Manin symbol [`X^i*Y^(k-2-i), (u,v)`], expressed in terms of the basis.
EXAMPLES::
sage: M = ModularSymbols(15,6,sign=-1) sage: R = M.manin_gens_to_basis() sage: sage.modular.modsym.heilbronn.hecke_images_gamma0_weight_k(4,1,3,15,6,[1,11,12], R) [ 0 0 1/8 -1/8 0 0 0 0] [-4435/22 -1483/22 -112 -4459/22 2151/22 -5140/11 4955/22 2340/11] [ 1253/22 1981/22 -2 3177/22 -1867/22 6560/11 -7549/22 -612/11] sage: x = M((3,4,1)) ; x.element() (0, 0, 1/8, -1/8, 0, 0, 0, 0) sage: M.T(11)(x).element() (-4435/22, -1483/22, -112, -4459/22, 2151/22, -5140/11, 4955/22, 2340/11) sage: M.T(12)(x).element() (1253/22, 1981/22, -2, 3177/22, -1867/22, 6560/11, -7549/22, -612/11) """
# The Manin symbols are enumerated as # all [0, (u,v)] for (u,v) in P^1(N) then # all [1, (u,v)] for (u,v) in P^1(N) etc. # So we create a zero dense matrix over QQ with len(indices) rows # and #P^1(N) * (k-1) columns. cdef Matrix_rational_dense T
R = R.change_ring(QQ)
cdef Py_ssize_t j, m, z, w, n, p cdef int *a cdef int *b
cdef Heilbronn H cdef fmpz_poly_t* poly
# Allocate memory to hold images of (u,v) under all Heilbronn matrices
# Compute images of (u,v) under all Heilbronn matrices
# Compute images of X^i Y^(2-k-i) under each Heilbronn matrix
# The following line dominates the runtime of this entire function:
# Compute the indexes of these images. # We just store them in the array a for simplicity. # Compute index of the symbol a[j], b[j] in the standard list. # Now fill in row z of the matrix T. # The following two lines are just a vastly faster version of: # T[z, n*w + p] += poly[j][w] # They use that we know that T has only integer entries. fmpq_mat_entry_num(T._matrix, z, n*w+p), fmpz_poly_get_coeff_ptr(poly[j], w))
# Free a and b
# Free poly part
# Return the product T * R, whose rows are the image of (u,v) under # the Hecke operators T_n for n in indices.
############################################################################ # Fast application of Heilbronn operators to computation of # systems of Hecke eigenvalues. # Nontrivial character of order > 2 and weight > 2 ############################################################################
# TODO
############################################################################ # Fast application of Heilbronn operators to computation of # systems of Hecke eigenvalues. # Nontrivial character of order = 2 and weight > 2 ############################################################################
# TODO |