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
8
#  it under the terms of the GNU General Public License as published by
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
"""
24
Highly optimized computation of the modular symbols map associated to
25
a simple new space of modular symbols.
26

27
AUTHOR:
28
    - William Stein
29
"""
30

31
from sage.matrix.all import matrix
32
from sage.rings.all import QQ, ZZ, Integer
33

34
include 'stdsage.pxi'
35

36
cdef enum:
37
    MAX_CONTFRAC = 100
38
    MAX_DEG = 10000
39

40
from sage.modular.modsym.p1list cimport P1List
41

42

43
def 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

61
cdef 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

91
cdef 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

258