CoCalc Shared Filescode / modular_symbol_map.pyx
1#################################################################################
2#
3# (c) Copyright 2010 William Stein
4#
5#  This file is part of PSAGE
6#
7#  PSAGE is free software: you can redistribute it and/or modify
9#  the Free Software Foundation, either version 3 of the License, or
10#  (at your option) any later version.
11#
12#  PSAGE is distributed in the hope that it will be useful,
13#  but WITHOUT ANY WARRANTY; without even the implied warranty of
14#  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
15#  GNU General Public License for more details.
16#
17#  You should have received a copy of the GNU General Public License
18#  along with this program.  If not, see <http://www.gnu.org/licenses/>.
19#
20#################################################################################
21
22
23"""
24Highly optimized computation of the modular symbols map associated to
25a simple new space of modular symbols.
26
27AUTHOR:
28    - William Stein
29"""
30
31from sage.matrix.all import matrix
32from sage.rings.all import QQ, ZZ, Integer
33
34include 'stdsage.pxi'
35
36cdef enum:
37    MAX_CONTFRAC = 100
38    MAX_DEG = 10000
39
40from sage.modular.modsym.p1list cimport P1List
41
42
43def test_contfrac_q(a, b):
44    """
45    EXAMPLES::
46
47        sage: import psage.modform.rational.modular_symbol_map
48        sage: psage.modform.rational.modular_symbol_map.test_contfrac_q(7, 97)
49        [1, 13, 14, 97]
50        sage: continued_fraction(7/97).convergents()
51        [0, 1/13, 1/14, 7/97]
52        sage: psage.modform.rational.modular_symbol_map.test_contfrac_q(137, 93997)
53        [1, 686, 6175, 43911, 93997]
54        sage: continued_fraction(137/93997).convergents()
55        [0, 1/686, 9/6175, 64/43911, 137/93997]
56    """
57    cdef long qi[MAX_CONTFRAC]
58    cdef int n = contfrac_q(qi, a, b)
59    return [qi[i] for i in range(n)]
60
61cdef long contfrac_q(long qi[MAX_CONTFRAC], long a, long b) except -1:
62    """
63    Given coprime integers a, b, compute the sequence q_i of
64    denominators in the continued fraction expansion of a/b, storing
65    the results in the array q.
66
67    Returns the number of q_i.  The rest of the q array is
68    garbage.
69    """
70    cdef long c
71
72    qi[0] = 1
73    if a == 0:
74        return 1
75    if b == 0:
76        raise ZeroDivisionError
77
78    # one iteration isn't needed, since we are only computing the q_i.
79    a,b = b, a%b
80    cdef int i=1
81    while b:
82        if i >= 2:
83            qi[i] = qi[i-1]*(a//b) + qi[i-2]
84        else:
85            qi[1] = a//b
86        a,b = b, a%b
87        i += 1
88
89    return i
90
91cdef class ModularSymbolMap:
92    cdef long d, N
93    cdef public long denom
94    cdef long* X  # coefficients of linear map from P^1 to Q^d.
95    cdef public object C
96    cdef public P1List P1
97
98    def __cinit__(self):
99        self.X = NULL
100
101    def __repr__(self):
102        return "Modular symbols map for modular symbols factor of dimension %s and level %s"%(self.d, self.N)
103
104    def __init__(self, A):
105        """
106        EXAMPLES::
107
108        We illustrate a few ways to construct the modular symbol map for 389a.
109
110            sage: from psage.modform.rational.modular_symbol_map import ModularSymbolMap
111            sage: A = ModularSymbols(389,sign=1).cuspidal_subspace().new_subspace().decomposition()[0]
112            sage: f = ModularSymbolMap(A)
113            sage: f._eval1(-3,7)
114            [-2]
115            sage: f.denom
116            2
117
118            sage: E = EllipticCurve('389a'); g = E.modular_symbol()
119            sage: h = ModularSymbolMap(g)
120            sage: h._eval1(-3,7)
121            [-2]
122            sage: h.denom
123            1
124            sage: g(-3/7)
125            -2
126
127            sage: f.denom*f.C == h.denom*h.C
128            True
129        """
130        import sage.schemes.elliptic_curves.ell_modular_symbols as ell_modular_symbols
131        from sage.modular.modsym.space import ModularSymbolsSpace
132
133        if isinstance(A, ell_modular_symbols.ModularSymbol):
134            C = self._init_using_ell_modular_symbol(A)
135
136        elif isinstance(A, ModularSymbolsSpace):
137            C = self._init_using_modsym_space(A)
138
139        else:
140
141            raise NotImplementedError, "Creating modular symbols from object of type %s not implemented"%type(A)
142
143        # useful for debugging only -- otherwise is a waste of memory/space
144        self.C = C
145
146        X, self.denom = C._clear_denom()
147        # Now store in a C data structure the entries of X (as long's)
148        self.X = <long*>sage_malloc(sizeof(long*)*X.nrows()*X.ncols())
149        cdef Py_ssize_t i, j, n
150        n = 0
151        for a in X.list():
152            self.X[n] = a
153            n += 1
154
155    def _init_using_ell_modular_symbol(self, f):
156        # Input f is an elliptic curve modular symbol map
157        assert f.sign() != 0
158        self.d = 1  # the dimension
159
160        E = f.elliptic_curve()
161        self.N = E.conductor()
162        self.P1 = P1List(self.N)
163
164        # Make a matrix whose rows are the images of the Manin symbols
165        # corresponding to the elements of P^1 under f.
166        n = len(self.P1)
167        C = matrix(QQ, n, 1)
168        for i in range(n):
169            # If the ith element of P1 is (u,v) for some u,v modulo N.
170            # We need to turn this into something we can evaluate f on.
171            # 1. Lift to a 2x2 SL2Z matrix M=(a,b;c,d) with c=u, d=v (mod N).
172            # 2. The modular symbol is M{0,oo} = {b/d,a/c}.
173            # 3. So {b/d, a/c} = {b/d,oo} + {oo,a/c} = -{oo,b/d} + {oo,a/c} = f(a/c)-f(b/d).
174            # 4. Thus x |--> f(a/c)-f(b/d).
175            a,b,c,d = self.P1.lift_to_sl2z(i)  # output are Python ints, so careful!
176            C[i,0] = (f(Integer(a)/c) if c else 0) - (f(Integer(b)/d) if d else 0)
177        return C
178
179    def _init_using_modsym_space(self, A):
180        # Very slow generic setup code.  This is "O(1)" in that we
181        # care mainly about evaluation time being fast, at least in
182        # this code.
183        if A.sign() == 0:
184            raise ValueError, "A must have sign +1 or -1"
185        self.d = A.dimension()
186
187        self.N = A.level()
188        self.P1 = P1List(self.N)
189
190        # The key data we need from the modular symbols space is the
191        # map that assigns to an element of P1 the corresponding
192        # element of ZZ^n.  That's it.  We forget everything else.
193        M = A.ambient_module()
194        W = matrix([M.manin_symbol(x).element() for x in self.P1])
195        B = A.dual_free_module().basis_matrix().transpose()
196
197        # Return matrix whose rows are the images of the Manin symbols
198        # corresponding to the elements of P^1 under the modular
199        # symbol map.
200        return W*B
201
202
203    def __dealloc__(self):
204        if self.X:
205            sage_free(self.X)
206
207    cdef int evaluate(self, long v[MAX_DEG], long a, long b) except -1:
208        cdef long q[MAX_CONTFRAC]
209        cdef int i, j, k, n, sign=1
210        cdef long* x
211
212        # initialize answer vector to 0
213        for i in range(self.d):
214            v[i] = 0
215
216        # compute continued fraction
217        n = contfrac_q(q, a, b)
218
219        # compute corresponding modular symbols, mapping over...
220        for i in range(1,n):
221            j = self.P1.index((sign*q[i])%self.N, q[i-1]%self.N)
222            # map over, adding a row of the matrix self.X
223            # to the answer vector v.
224            x = self.X + j*self.d
225            for k in range(self.d):
226                v[k] += x[k]
227            # change sign, so q[i] is multiplied by (-1)^(i-1)
228            sign *= -1
229
230    def dimension(self):
231        return self.d
232
233    def _eval0(self, a, b):
234        cdef long v[MAX_DEG]
235        self.evaluate(v, a, b)
236
237    def _eval1(self, a, b):
238        """
239        EXAMPLE::
240
241            sage: from psage.modform.rational.modular_symbol_map import ModularSymbolMap
242            sage: A = ModularSymbols(188,sign=1).cuspidal_subspace().new_subspace().decomposition()[-1]
243            sage: f = ModularSymbolMap(A)
244            sage: f._eval1(-3,7)
245            [-3, 0]
246
247        """
248        cdef long v[MAX_DEG]
249        self.evaluate(v, a, b)
250        cdef int i
251        return [v[i] for i in range(self.d)]
252
253
254
255
256
257