Ticket #23342: wrapper.spyx

File wrapper.spyx, 4.9 KB (added by kedlaya, 2 years ago)

Sage+Cython wrapper to Hein's C++ code

Line 
1# Copyright (C) 2017 by Kiran S. Kedlaya <kskedl@gmail.com>
2#
3# distutils: language = c++
4#encoding=utf8
5
6#clang C++
7#clib gmpxx gmp m pthread
8#cargs -g -Wall -O3 -Wextra -Werror -pedantic -std=c++11
9#cinclude /usr/local/include
10#cfile AutomorphismZZ.cpp CharacterZZ.cpp Eigenvector.cpp GenusZZ.cpp IsometryZZ.cpp MathZZ.cpp NeighborIteratorZZ.cpp QuadFormZZ.cpp SparseMatrix.cpp
11
12import sys
13reload(sys)
14sys.setdefaultencoding('utf8')
15
16from cython.operator cimport dereference as deref
17
18from libcpp.map cimport map
19from libcpp.unordered_map cimport unordered_map
20from libcpp.memory cimport shared_ptr, make_shared
21
22from sage.rings.integer cimport Integer
23from sage.libs.gmp.types cimport mpz_t
24
25from sage.quadratic_forms.ternary_qf import find_a_ternary_qf_by_level_disc
26from sage.arith.misc import squarefree_divisors
27from sage.matrix.constructor import matrix
28   
29cdef extern from "gmpxx.h":
30    cdef cppclass mpz_class:
31        mpz_class(mpz_t a)
32       
33    cdef cppclass mpq_class:
34        pass
35
36cdef extern from "QuadForm.h":
37    cdef cppclass QuadForm[R,F]:
38        QuadForm(const R& a, const R& b, const R& c,
39            const R& f, const R& g, const R& h, long reduced=false)
40
41cdef extern from "SparseMatrix.h":
42    cdef cppclass SparseMatrix:
43        ctypedef map[long, map[long, long]] dataMap
44        const dataMap& data() const;
45        long num_rows() const;
46
47cdef extern from "Character.h":
48    cdef cppclass Character[R,F]:
49        Character(const R& cond)
50
51cdef extern from "Genus.h":
52    cdef cppclass Genus[R,F]:
53        ctypedef shared_ptr[QuadForm[R,F]] QuadFormPtr
54        ctypedef shared_ptr[HeckeOperator[R,F]] HeckePtr
55
56        Genus(const QuadForm[R,F]& q, long numThreads=0)
57        void add_character(const Character[R,F]& chi)
58        void compute_hecke_operators(const R& p, long numThreads)
59        HeckePtr hecke_operator(const R& p, const R& cond)
60
61cdef extern from "HeckeOperator.h":
62    cdef cppclass HeckeOperator[R,F]:
63        SparseMatrix& matrix() const;
64         
65ctypedef Genus[mpz_class, mpq_class] GenusZZ
66ctypedef QuadForm[mpz_class, mpq_class] QuadFormZZ
67ctypedef Character[mpz_class, mpq_class] CharacterZZ
68
69ctypedef shared_ptr[QuadFormZZ] QuadFormZZPtr
70ctypedef shared_ptr[GenusZZ] GenusZZPtr
71
72cdef GenusZZPtr create_genus_from_coeffs(Integer a, Integer b,
73                                          Integer c, Integer f,
74                                          Integer g, Integer h):
75    cdef QuadFormZZPtr q
76    cdef GenusZZPtr genus
77    q = make_shared[QuadFormZZ](mpz_class(a.value), mpz_class(b.value),
78                                mpz_class(c.value), mpz_class(f.value),
79                                mpz_class(g.value), mpz_class(h.value), False)
80    genus = make_shared[GenusZZ](deref(q))
81    return(genus)
82
83cpdef hecke_birch(N, l):
84    r"""
85    Compute Hecke operators at level N for each prime in l via Birch's method.
86
87    INPUT:
88 
89    - ``N``: a squarefree positive integer
90
91    - ``l``: a list of prime positive integers not dividing N
92
93    OUTPUT:
94
95    A dict with one entry for each pair (d, p) where d is a squarefree
96    divisor of N and p is an entry of l. This entry is a sparse integer matrix
97    representing the action of the Hecke operator T_p on the subspace of the
98    space S_2(Gamma_0(N), QQ)^{new} on which the Atkin-Lehner involutions act
99    via the Kronecker character of level d.
100
101    EXAMPLES::
102
103        sage: h = hecke_birch(91, [2, 3])
104        sage: h[(13, 3)]
105        [-2]
106
107    TESTS:
108
109    Compare with Brandt symbols:
110
111        sage: B = BrandtModule(151)
112        sage: h = hecke_birch(151, [2])
113        sage: P = prod(h[i].charpoly() for i in h)
114        sage: P == B.hecke_polynomial(2)
115        True
116
117    """
118    cdef Integer d, p
119    cdef SparseMatrix mat
120    cdef SparseMatrix.dataMap dat
121
122    if not N.is_squarefree():
123        raise ValueError, "level " + str(N) + \
124            " is not squarefree"
125    for p in l:
126        if not p.is_prime:
127            raise ValueError, "index " + str(p) + " is not prime"
128        if N % p == 0:
129            raise ValueError, "prime " + str(p) + \
130                " divides level " + str(N)
131    ans = {}
132    Q = find_a_ternary_qf_by_level_disc(4*N, N)
133    (a,b,c,f,g,h) = Q.coefficients()
134    cdef GenusZZPtr genus = create_genus_from_coeffs(a,b,c,f,g,h)
135    for d in squarefree_divisors(N):
136        deref(genus).add_character(CharacterZZ(mpz_class(d.value)))
137           
138    for p in l:
139        deref(genus).compute_hecke_operators(mpz_class(p.value), 0)
140        for d in squarefree_divisors(N):
141            hecke = deref(genus).hecke_operator(mpz_class(p.value),
142                                                      mpz_class(d.value))
143            mat = deref(hecke).matrix()
144            dat = mat.data()
145            datdict = {}
146            for (i,dict1) in dat:
147                for j in dict1:
148                    datdict[(i,j)] = Integer(dat[i][j])
149            n = mat.num_rows()
150            ans[(d, p)] = matrix(n,n,datdict)
151    return(ans)
152
153
154