source: sage/schemes/elliptic_curves/padic_lseries.py @ 4784:b37da68b32d4

Revision 4784:b37da68b32d4, 27.5 KB checked in by William Stein <wstein@…>, 6 years ago (diff)

Fix doctests.
Line 
1"""
2p-adic L-functions of elliptic curves
3
4AUTHORS:
5   -- William Stein (2007-01-01): first version
6   -- chris wuthrich (22/05/2007): changed minor issued and added supersingular things
7"""
8
9######################################################################
10#       Copyright (C) 2007 William Stein <wstein@gmail.com>
11#
12#  Distributed under the terms of the GNU General Public License (GPL)
13#
14#    This code is distributed in the hope that it will be useful,
15#    but WITHOUT ANY WARRANTY; without even the implied warranty of
16#    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
17#    General Public License for more details.
18#
19#  The full text of the GPL is available at:
20#
21#                  http://www.gnu.org/licenses/
22######################################################################
23
24
25from sage.rings.integer_ring import   ZZ
26from sage.rings.rational_field import QQ
27from sage.rings.padics.factory import Qp
28from sage.rings.infinity import infinity
29from sage.rings.all import PowerSeriesRing, PolynomialRing, Integers
30
31from sage.rings.integer import Integer
32from sage.rings.arith import valuation, binomial
33
34from sage.structure.sage_object import SageObject
35
36from sage.misc.all import verbose, denominator
37from sage.databases.cremona import parse_cremona_label
38from sage.schemes.elliptic_curves.constructor import EllipticCurve
39import sage.rings.arith as arith
40
41from sage.modules.free_module_element import vector
42import sage.matrix.all as matrix
43import monsky_washnitzer
44from sage.interfaces.all import gp
45from sage.misc.functional import log
46
47class pAdicLseries(SageObject):
48    """
49    The p-adic L-series of an elliptic curve.
50
51    EXAMPLES:
52    An ordinary example:
53        sage: e = EllipticCurve('389a')
54        sage: L = e.padic_lseries(5)
55        sage: L.series(0)
56        Traceback (most recent call last):
57        ...
58        ValueError: n (=0) must be a positive integer
59        sage: L.series(1)
60        O(T^1)
61        sage: L.series(2)
62        O(5^4) + O(5^1)*T + (4 + O(5))*T^2 + (2 + O(5))*T^3 + (3 + O(5))*T^4 + O(T^5)
63        sage: L.series(3, prec=10)
64        O(5^5) + O(5^2)*T + (4 + 4*5 + O(5^2))*T^2 + (2 + 4*5 + O(5^2))*T^3 + (3 + O(5^2))*T^4 + (1 + O(5))*T^5 + (3*5 + O(5^2))*T^6 + (4 + 5 + O(5^2))*T^7 + (2 + 5 + O(5^2))*T^8 + O(5^2)*T^9 + O(T^10)
65
66    A prime p such that E[p] is reducible:
67        sage: L = EllipticCurve('11a').padic_lseries(5)
68        sage: L.series(1)
69        5 + O(5^2) + O(T)
70        sage: L.series(2)
71        5 + 4*5^2 + O(5^3) + O(5^0)*T + O(5^0)*T^2 + O(5^0)*T^3 + O(5^0)*T^4 + O(T^5)
72        sage: L.series(3)
73        5 + 4*5^2 + 4*5^3 + O(5^4) + O(5^1)*T + O(5^1)*T^2 + O(5^1)*T^3 + O(5^1)*T^4 + O(T^5)       
74    """
75    def __init__(self, E, p, normalize):
76        """
77        INPUT:
78            E -- an elliptic curve
79            p -- a prime of good reduction
80            normalize -- (bool, default: True); whether or not to correctly
81                 normalize the L-series, up to a power of -1 and 2.
82                 If False computations may be faster.
83        """
84        self._E = E
85        self._p = ZZ(p)
86        self._normalize = normalize
87        if not self._p.is_prime():
88            raise ValueError, "p (=%s) must be a prime"%p
89        if E.conductor() % (self._p)**2 == 0:
90            raise NotImplementedError, "p (=%s) must be a prime of semi-stable reduction"%p
91       
92        # this factor adjusts the p-adic L-series so that it is correct for any element of the isogeny class
93        crla = parse_cremona_label(E.label()) 
94        cr0 = Integer(crla[0]).str() + crla[1] + '1'
95        E0 = EllipticCurve(cr0)
96        self._quotient_of_periods = QQ(E0.period_lattice()[0]/E.period_lattice()[0])
97        self._modular_symbol = E.modular_symbol(sign=1, normalize=normalize)
98
99    def elliptic_curve(self):
100        """
101        Return the elliptic curve to which this p-adic L-series is associated.
102       
103        EXAMPLES:
104            sage: L = EllipticCurve('11a').padic_lseries(5)
105            sage: L.elliptic_curve()
106            Elliptic Curve defined by y^2 + y = x^3 - x^2 - 10*x - 20 over Rational Field
107        """
108        return self._E
109
110    def prime(self):
111        """
112        EXAMPLES:
113            sage: L = EllipticCurve('11a').padic_lseries(5)
114            sage: L.prime()
115            5
116        """
117        return self._p
118
119    def _repr_(self):
120        """
121        Return print representation.
122
123            sage: e = EllipticCurve('37a')
124            sage: e.padic_lseries(3)
125            3-adic L-series of Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field
126            sage: e.padic_lseries(3,normalize=False)
127            3-adic L-series of Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field (not normalized)
128            sage: L = e.padic_lseries(3,normalize=False)
129            sage: L.rename('(factor)*L_3(T)')
130            sage: L
131            (factor)*L_3(T)
132        """
133        s = "%s-adic L-series of %s"%(self._p, self._E)
134        if not self._normalize:
135            s += ' (not normalized)'
136        return s
137
138    def modular_symbol(self, r):
139        """
140        Return the modular symbol used to compute this p-adic
141        L-series evaluated at r.
142
143        EXAMPLES:
144            sage: L = EllipticCurve('11a').padic_lseries(5)
145            sage: [L.modular_symbol(r) for r in [0,1/5,oo,1/11]]
146            [1/5, 6/5, 0, 0]       
147        """
148        return self._modular_symbol(r)
149
150    def measure(self, a, n, prec):
151        r"""
152        Return the measure on $\ZZ_p^*$ defined by
153           $$
154             \mu_{E,\alpha}^+ ( a + p^n \ZZ_p  ) =
155                   \frac{1}{\alpha^n} \modsym{a}{p^n} - \frac{1}{\alpha^{n+1}} \modsym{a}{p^{n-1}}
156           $$
157        that is used to define this $p$-adic $L$-function.
158
159        INPUT:
160            a -- an integer
161            n -- a non-negative integer
162            prec -- an integer
163
164        EXAMPLES:
165            sage: E = EllipticCurve('37a')
166            sage: L = E.padic_lseries(5)
167            sage: L.measure(1,2, prec=9)
168            1 + 4*5 + 2*5^2 + 4*5^3 + 3*5^4 + 5^5 + 4*5^6 + 4*5^7 + 4*5^8 + O(5^9)
169        """
170        p = self._p
171        alpha = self.alpha(prec=prec)
172        z = 1/(alpha**n)
173        w = p**(n-1)
174        f = self._modular_symbol
175        if self._E.conductor() % p == 0:
176            return z * f(a/(p*w))
177        return z * f(a/(p*w)) - (z/alpha) * f(a/w)
178       
179    def alpha(self, prec=20):
180        r"""
181        Return a p-adic root $\alpha$ of the polynomial $x^2 - a_p x
182        + p$ with $\ord_p(\alpha) < 1$.  In the ordinary case this is
183        just the unit root.
184
185        INPUT:
186            prec -- positive integer, the p-adic precision of the root.
187
188        EXAMPLES:
189        Consider the elliptic curve 37a:
190            sage: E = EllipticCurve('37a')
191
192        An ordinary prime:
193            sage: L = E.padic_lseries(5)
194            sage: alpha = L.alpha(10); alpha
195            3 + 2*5 + 4*5^2 + 2*5^3 + 5^4 + 4*5^5 + 2*5^7 + 5^8 + 5^9 + O(5^10)
196            sage: alpha^2 - E.ap(5)*alpha + 5
197            O(5^10)
198
199        A supersingular prime.
200            sage: L = E.padic_lseries(3)
201            sage: alpha = L.alpha(10); alpha
202            (1 + O(3^10))*alpha
203            sage: alpha^2 - E.ap(3)*alpha + 3
204            (O(3^10))*alpha^2 + (O(3^11))*alpha + (O(3^11))
205
206        A reducible prime:
207            sage: L = EllipticCurve('11a').padic_lseries(5)
208            sage: L.alpha(5)
209            1 + 4*5 + 3*5^2 + 2*5^3 + 4*5^4 + O(5^5)
210        """
211        try:
212            return self._alpha[prec]
213        except AttributeError:
214            self._alpha = {}
215        except KeyError:
216            pass
217        E = self._E
218        p = self._p
219        a_p = E.ap(p)
220        K = Qp(p, prec, print_mode='series')
221 
222        if E.conductor() % p == 0:
223            self._alpha[prec] = K(a_p)
224            return K(a_p) 
225
226        R = ZZ['x']
227        f = R([p, -a_p, 1])
228        if E.is_ordinary(p):
229            G = f.factor_padic(p, prec+5)
230            for pr, e in G:
231                a = -pr[0]
232                if a.valuation() < 1:
233                    self._alpha[prec] = K(a)
234                    return K(a)
235            raise ValueError, "bug in p-adic L-function alpha"
236        else: # supersingular case
237            f = f.change_ring(Qp(p, prec, print_mode='series'))
238            a = f.root_field('alpha', check_irreducible=False).gen()
239            self._alpha[prec] = a
240            return a
241
242    def order_of_vanishing(self, proof=True):
243        """
244        Return the order of vanishing of this $p$-adic $L$-series.
245
246        The output of this function is provably correct, due to a
247        theorem of Kato.  This function will terminate if and only if
248        the Mazur-Tate-Teitelbaum analogue of the BSD conjecture about
249        the rank of the curve is true and the subgroup of elements of
250        p-power order in the Shafarevich-Tate group of this curve is
251        finite.  I.e., if this function terminates (with no errors!),
252        then you may conclude that the p-adic BSD rank conjecture is
253        true and that the p-part of Sha is finite.
254
255        NOTE: currently $p$ must be a prime of good ordinary reduction.
256
257        EXAMPLES:
258            sage: L = EllipticCurve('11a').padic_lseries(3)
259            sage: L.order_of_vanishing()
260            0
261            sage: L = EllipticCurve('11a').padic_lseries(5)                       
262            sage: L.order_of_vanishing()
263            0       
264            sage: L = EllipticCurve('37a').padic_lseries(5)
265            sage: L.order_of_vanishing()
266            1
267            sage: L = EllipticCurve('43a').padic_lseries(3)
268            sage: L.order_of_vanishing()
269            1
270            sage: L = EllipticCurve('37b').padic_lseries(3)
271            sage: L.order_of_vanishing()
272            0
273
274        We verify that Sha(E)(p) is finite for p=3,5,7 for the
275        first curve of rank 2:
276            sage: e = EllipticCurve('389a')
277            sage: for p in primes(3,10):
278            ...    print p, e.padic_lseries(p).order_of_vanishing()
279            3 2
280            5 2
281            7 2           
282        """
283        try:
284            return self.__ord
285        except AttributeError:
286            pass
287       
288        if not self.is_ordinary():
289            raise NotImplementedError
290        E = self.elliptic_curve()
291        if not E.is_good(self.prime()):
292            raise ValueError, "prime must be of good reduction"
293        r = E.rank()
294        n = 1
295        while True:
296            f = self.series(n)
297            v = f.valuation()
298            if v < r:
299                raise RuntimeError, "while computing p-adic order of vanishing, got a contradiction: the curve is %s, the curve has rank %s, but the p-adic L-series vanishes to order <= %s"%(E, r, v)
300            if v == r:
301                self.__ord = v
302                return v
303            n += 1
304
305
306    def _c_bounds(self, n):
307        raise NotImplementedError
308
309    def _prec_bounds(self, n):
310        raise NotImplementedError
311
312    def teichmuller(self, prec):
313        r"""
314        Return Teichmuller lifts to the given precision.
315       
316        INPUT:
317            prec -- a positive integer.
318           
319        OUTPUT:
320            the cached Teichmuller lifts
321
322        EXAMPLES:
323            sage: L = EllipticCurve('11a').padic_lseries(7)
324            sage: L.teichmuller(1)
325            [0, 1, 2, 3, 4, 5, 6]
326            sage: L.teichmuller(2)
327            [0, 1, 30, 31, 18, 19, 48]
328        """
329        p = self._p
330        K = Qp(p, prec, print_mode='series')       
331        return [Integer(0)] + \
332               [a.residue(prec).lift() for a in K.teichmuller_system()]
333
334    def _e_bounds(self, n):
335        p = self._p
336        T = (ZZ['T']).gen()
337        w = (1+T)**(p**n) - 1
338        return [infinity] + [valuation(w[j],p) for j in range(1,w.degree()+1)]
339       
340    def _get_series_from_cache(self, n, prec):
341        try:
342            return self.__series[(n,prec)]
343        except AttributeError:
344            self.__series = {}
345        except KeyError:
346            for _n, _prec in self.__series.keys():
347                if _n == n and _prec <= prec:
348                    return self.__series[(_n,_prec)].add_bigoh(prec)
349        return None
350
351    def _set_series_in_cache(self, n, prec, f):
352        self.__series[(n,prec)] = f
353       
354
355class pAdicLseriesOrdinary(pAdicLseries):
356    def series(self, n=2, prec=5):
357        r"""
358        Return the $n$-th approximation to the $p$-adic $L$-series as
359        a power series in $T$ (corresponding to $\gamma-1$ with
360        $\gamma=1+p$ as a generator of $1+p\mathbb{Z}_p$).  Each
361        coefficient is a $p$-adic number whose precision is provably
362        correct.
363
364        INPUT:
365            n -- (default: 2) a positive integer
366            prec -- (default: 5) maxima number of terms of the series
367                    to compute; to compute as many as possible just
368                    give a very large number for prec; the result will
369                    still be correct.
370
371        ALIAS:
372            power_series is identical to series.
373
374        EXAMPLES:
375        We compute some $p$-adic $L$-functions associated to the elliptic
376        curve 11a.
377       
378            sage: E = EllipticCurve('11a')
379            sage: p = 3
380            sage: E.is_ordinary(p)
381            True
382            sage: L = E.padic_lseries(p)
383            sage: L.series(3)
384            2 + 3 + 3^2 + 2*3^3 + O(3^5) + (1 + 3 + O(3^2))*T + (1 + 2*3 + O(3^2))*T^2 + O(3^1)*T^3 + (2*3 + O(3^2))*T^4 + O(T^5)
385
386
387        Another example at a prime of bad reduction, where the
388        $p$-adic $L$-function has an extra 0 (compared to the non
389        $p$-adic $L$-function).
390       
391            sage: E = EllipticCurve('11a')
392            sage: p = 11
393            sage: E.is_ordinary(p)
394            True
395            sage: L = E.padic_lseries(p)
396            sage: L.series(2)
397            O(11^4) + (10 + O(11))*T + (6 + O(11))*T^2 + (2 + O(11))*T^3 + (5 + O(11))*T^4 + O(T^5)
398           
399        We compute a $p$-adic $L$-function that vanishes to order $2$.
400            sage: E = EllipticCurve('389a')
401            sage: p = 3
402            sage: E.is_ordinary(p)
403            True
404            sage: L = E.padic_lseries(p)
405            sage: L.series(1)
406            O(T^1)
407            sage: L.series(2)
408            O(3^4) + O(3^1)*T + (2 + O(3))*T^2 + O(T^3)
409            sage: L.series(3)
410            O(3^5) + O(3^2)*T + (2 + 2*3 + O(3^2))*T^2 + (2 + O(3))*T^3 + (1 + 3 + O(3^2))*T^4 + O(T^5)
411        """
412        n = ZZ(n)
413        if n < 1:
414            raise ValueError, "n (=%s) must be a positive integer"%n
415       
416
417        p = self._p
418       
419        bounds = self._prec_bounds(n)
420        padic_prec = max(bounds[1:]) + 5
421        res_series_prec = min(p**(n-1), prec)
422       
423        ans = self._get_series_from_cache(n, res_series_prec)
424        if not ans is None:
425            return ans
426 
427        K = QQ
428        gamma = K(1 + p)
429        R = PowerSeriesRing(K,'T',res_series_prec)
430        T = R(R.gen(),res_series_prec )
431        L = R(0) 
432        one_plus_T_factor = R(1) 
433        gamma_power = 1
434        teich = self.teichmuller(padic_prec)
435        for j in range(p**(n-1)):
436            s = K(0) 
437            for a in range(1,p):
438                b = teich[a] * gamma_power
439                s += self.measure(b, n, padic_prec).lift()
440            L += s * one_plus_T_factor
441            one_plus_T_factor *= 1+T
442            gamma_power *= gamma
443
444        # Now create series but with each coefficient truncated
445        # so it is proven correct:
446        K = Qp(p, padic_prec, print_mode='series')
447        R = PowerSeriesRing(K,'T',res_series_prec)
448        L = R(L,res_series_prec)
449        aj = L.list()
450        if len(aj) > 0:
451            aj = [aj[0].add_bigoh(padic_prec-2)] + [aj[j].add_bigoh(bounds[j]) for j in range(1,len(aj))]
452        L = R(aj,res_series_prec ) * self._quotient_of_periods
453       
454        self._set_series_in_cache(n, res_series_prec, L)
455       
456        return L
457
458    power_series = series
459
460    def is_ordinary(self):
461        """
462        Return True if the elliptic that this $L$-function is attached
463        to is ordinary.
464       
465        EXAMPLES:
466            sage: L = EllipticCurve('11a').padic_lseries(5)           
467            sage: L.is_ordinary()
468            True
469        """
470        return True
471
472    def is_supersingular(self):
473        """
474        Return True if the elliptic that this L function is attached
475        to is supersingular.
476               
477        EXAMPLES:
478            sage: L = EllipticCurve('11a').padic_lseries(5)           
479            sage: L.is_supersingular()
480            False
481        """
482        return False
483
484    def _c_bound(self):
485        try:
486            return self.__c_bound
487        except AttributeError:
488            pass
489        E = self._E
490        p = self._p
491        if E.is_irreducible(p):
492            ans = 0
493        else:
494            m = E.modular_symbol_space(sign=1)
495            b = m.boundary_map().codomain()
496            C = b._known_cusps()  # all known, since computed the boundary map
497            ans = max([valuation(self.modular_symbol(a).denominator(), p) for a in C]) 
498        self.__c_bound = ans
499        return ans
500               
501    def _prec_bounds(self, n):
502        p = self._p
503        e = self._e_bounds(n-1)
504        c = self._c_bound()
505        return [e[j] - c for j in range(len(e))]
506
507
508class pAdicLseriesSupersingular(pAdicLseries):
509    def series(self, n=2, prec=5):
510        r"""
511        Return the $n$-th approximation to the $p$-adic $L$-series as a
512        power series in $T$ (corresponding to $\gamma-1$ with
513        $\gamma=1+p$ as a generator of $1+p\mathbb{Z}_p$).  Each
514        coefficient is a $p$-adic number whose precision is probably
515        {\em not} correct.
516
517        WARNING: ** The output precision is not as high as claimed! **
518
519        INPUT:
520            n -- (default: 2) a positive integer
521            prec -- (default: 5) maxima number of terms of the series
522                    to compute; to compute as many as possible just
523                    give a very large number for prec; the result will
524                    still be correct.
525
526        ALIAS:
527            power_series is identical to series.
528
529        EXAMPLES:
530        A superingular example, where we must compute to higher precision to see anything.
531            sage: e = EllipticCurve('37a')
532            sage: L = e.padic_lseries(3); L
533            3-adic L-series of Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field
534            sage: L.series(2)
535            O(T^3)
536            sage: L.series(4)         # takes a long time (several seconds)
537    (O(3^1))*alpha + (O(3^2)) + ((O(3^-1))*alpha + (2*3^-1 + O(3^0)))*T + ((O(3^-1))*alpha + (2*3^-1 + O(3^0)))*T^2 + ((O(3^-2))*alpha + (O(3^-1)))*T^3 + ((O(3^-1))*alpha + (3^-1 + O(3^0)))*T^4 + O(T^5)
538            sage: L.alpha(2).parent()
539            Univariate Quotient Polynomial Ring in alpha over 3-adic Field with capped relative precision 2 with modulus (1 + O(3^2))*x^2 + (3 + O(3^3))*x + (3 + O(3^3))
540        """
541        n = ZZ(n)
542        if n < 1:
543            raise ValueError, "n (=%s) must be a positive integer"%n
544       
545       
546        bounds = self._prec_bounds(n)
547        padic_prec = max(sum(bounds[1:],[])) + 5
548        p = self._p
549
550        prec = min(p**(n-1), prec)
551        ans = self._get_series_from_cache(n, prec)
552        if not ans is None:
553            return ans
554
555        alpha = self.alpha(prec=padic_prec)
556        K = alpha.parent()
557        gamma = 1 + p
558        R = PowerSeriesRing(K,'T',prec)
559        T = R(R.gen(), prec)
560        L = R(0)
561        one_plus_T_factor = R(1)
562        gamma_power = 1
563        teich = self.teichmuller(padic_prec)
564        for j in range(p**(n-1)):
565            s = K(0)
566            for a in range(1,p):
567                b = teich[a] * gamma_power
568                s += self.measure(b, n, padic_prec)
569            L += s * one_plus_T_factor
570            one_plus_T_factor *= 1+T
571            gamma_power *= gamma
572
573        # Now create series but with each coefficient truncated
574        # so it is proven correct:
575        L = R(L,prec) 
576        aj = L.list()
577        if len(aj) > 0:
578            bj = [aj[0][0].add_bigoh(padic_prec-2) + alpha * aj[0][1].add_bigoh(padic_prec-2)]
579            bj += [aj[j][0].add_bigoh(bounds[j][0]) + alpha * aj[j][1].add_bigoh(bounds[j][1]) for j in range(1,len(aj))]
580            L = R(bj, prec)
581        L = L * self._quotient_of_periods
582        self._set_series_in_cache(n, prec, L)
583        return L
584
585    power_series = series
586
587    def is_ordinary(self):
588        return False
589
590    def is_supersingular(self):
591        return True
592
593    def _prec_bounds(self, n):
594        p = self._p
595        e = self._e_bounds(n-1)
596        c0 = ZZ(n+2)/2
597        c1 = ZZ(n+3)/2
598        return [[infinity,infinity]] + [[(e[j] - c0).floor(), (e[j] - c1).floor()] for j in range(1,len(e))]
599
600
601    def Dp_valued_series(self, n=2, prec=5):
602        r"""
603        Returns a vector of two components which are p-adic power series.
604        The answer v is such that
605            $$(1-\varphi)^(-2)* L_p(E,T) = v[1] * \omega + v[2] * \eta$$
606        as an element of the Dieudonne module $D_p(E) = H^1_{dR}(E/\QQ_p)$ where
607        $\omega$ is the invariant differential and $\eta$ is $x\cdot\omega$.   
608        According to the p-adic BSD this function has a zero of order
609        rank(E(Q)) and it's leading term is
610        \begin{verbatim}
611           +- #Sha(E/Q) * Tamagawa product / Torsion^2 * padic height regulator with values in D_p(E).
612        \end{verbatim}
613
614        WARNING: ** The output precision is not as high as claimed! **       
615       
616        INPUT:
617            n -- (default: 2) a positive integer
618            prec -- (default: 5) a positive integer
619
620        EXAMPLES:
621            sage: E = EllipticCurve('14a')
622            sage: L = E.padic_lseries(5)
623            sage: L.Dp_valued_series(2)
624            (4 + 4*5^2 + O(5^3), 0)
625        """
626        E = self._E
627        p = self._p
628        lps = self.series(n, prec=prec)
629   
630        # now split up the series in two lps = G + H * alpha
631        R = lps.base_ring().base_ring() # Qp
632        QpT , T = PowerSeriesRing(R,'T',prec).objgen()
633        G = sum([R(lps[n][0])*T**n for n in range(0,lps.prec())])
634        H = sum([R(lps[n][1])*T**n for n in range(0,lps.prec())]) 
635   
636        # now compute phi
637        phi =  self.geometric_frob_on_Dp(p) 
638        phi_omega_0 = phi[0,0]
639        phi_omega_1 = phi[1,0]
640        R = phi_omega_0.parent()
641        lpv = vector([G  + R(E.ap(p))*H - R(p) * phi_omega_0* H , - R(p)*phi_omega_1*H])  # this is L_p
642        eps = (1-phi)**(-2) 
643        resu = lpv*eps.transpose()
644        return resu
645
646
647    def geometric_frob_on_Dp(self, prec=20):
648        r"""
649        This returns the geometric Frobenius $\varphi$ on the Diedonne module $D_p(E)$
650        with respect to the basis $\omega$, the invariant differential and $\eta=x\omega$.
651        It satisfies  $phi^2 - a_p/p*phi + 1/p = 0$.
652
653        EXAMPLES:
654            sage: E = EllipticCurve('14a')
655            sage: L = E.padic_lseries(5)
656            sage: F = L.geometric_frob_on_Dp(5)
657            sage: F
658            [ 3 + 4*5 + 3*5^2 + 4*5^3 + O(5^4) 2*5^-1 + 1 + 3*5 + 3*5^3 + O(5^4)]
659            [     2 + 4*5 + 5^2 + 5^4 + O(5^5)            2 + 5^2 + 5^4 + O(5^5)]
660            sage: F^2 - E.ap(5)/5 * F + 1/5
661            [O(5^4) O(5^3)]
662            [O(5^4) O(5^4)]
663        """
664        E = self._E
665        p = self._p
666        Ew = E.weierstrass_model()
667        adjusted_prec = monsky_washnitzer.adjusted_prec(p, prec)
668        modprecring = Integers(p**adjusted_prec)
669        output_ring = Qp(p, prec)
670        R, x = PolynomialRing(modprecring, 'x').objgen()
671        Q = x**3 + modprecring(Ew.a4()) * x + modprecring(Ew.a6())
672        trace = Ew.ap(p)
673        fr = monsky_washnitzer.matrix_of_frobenius(Q, p, adjusted_prec, trace)
674        fr = matrix.matrix(output_ring,2,2,fr)
675        a=fr[0,0]
676        b=fr[0,1]
677        c=fr[1,0]
678        d=fr[1,1]
679        usq = (Ew.discriminant()/E.discriminant()).nth_root(6)
680        r = (4*E.a2() + E.a1())/12*usq;
681        frn = matrix.matrix([[a+c*r,(b-a*r+d*r-r**2*c)/usq],[usq*c,d-r*c]])
682        return frn**(-1)
683
684
685
686    def bernardi_sigma_function(self, prec=20):
687        r"""
688        Return the  p-adic sigma function of Bernardi in terms of $z = log(t)$.
689        This is the same as padic_sigma with E2 = 0.
690
691        EXAMPLES:
692            sage: E = EllipticCurve('14a')
693            sage: L = E.padic_lseries(5)
694            sage: L.bernardi_sigma_function(5) # Todo: some sort of consistency check!?
695            z + 1/24*z^3 + 29/384*z^5 - 8399/322560*z^7 - 291743/92897280*z^9 + O(z^11)
696        """
697        E = self._E
698        p = self._p
699        wp_in_pari = gp(E).ellwp('z',prec + 5)
700   
701        # transform the series from pari to sage
702        Qz , z = PowerSeriesRing(QQ,'z',prec+5).objgen()
703        # we dropped the 1/z^2
704        wp = sum([ QQ(wp_in_pari.polcoeff(k)) * z**k for k in range(1,prec + 5)]) 
705        minusx = (E.a1()**2+4*E.a2())/12 - wp
706   
707        si = z * minusx.integral().integral().exp() 
708   
709        return si
710   
711
712    def Dp_valued_height(self,prec=20):
713        """
714        Returns the canonical $p$-adic height with values in the Dieudonne module $D_p(E)$.
715        It is defined to be
716            $$h_{\eta} \cdot \omega - h_{\omega} \cdot \eta$$
717        where $h_{\eta}$ is made out of the sigma function of Bernardi and
718        $h_{\omega}$ is $-log^2$.
719        """
720        E = self._E
721        p = self._p
722        Ehat = E.formal()
723        elog = Ehat.log(prec + Integer(3))
724   
725        # we will have to do it properly with David Harvey's _DivPolyContext(E, R, Q)
726        n = arith.LCM(E.tamagawa_numbers())
727        n = arith.LCM(n, E.Np(p)) # allowed here because E has good reduction at p
728   
729        def height(P,check=True): 
730            if check:
731                assert P.curve() == E, "the point P must lie on the curve from which the height function was created"
732            Q = n * P
733            tt = - Q[0]/Q[1]
734            R = Qp(p,prec+5)
735            tt = R(tt)
736            zz = elog(tt)
737
738            homega = -zz**2/n**2
739       
740            eQ = denominator(Q[1])/denominator(Q[0])
741            si = self.bernardi_sigma_function(prec=prec+4)
742            heta =  2 * log(si(zz)/eQ) / n**2 
743       
744            R = Qp(p,prec)
745   
746            vec = vector([R(heta),-R(homega)]) 
747            return vec
748
749        return height
750   
751   
752   
753    def Dp_valued_regulator(self,prec=20):
754        """
755        Returns the canonical $p$-adic regulator with values in the Dieudonne module $D_p(E)$
756        as defined by Perrin-Riou using the canonical $p$-adic height.
757        """
758   
759        p = self._p
760        E = self._E
761       
762        h =  self.Dp_valued_height(prec=prec)
763   
764        # this is the height_{v} (P) for a v in D_p
765        def hv(vec,P):
766            hP = h(P)
767            return - vec[0]*hP[1] +vec[1]*hP[0]
768       
769        #    def hvpairing(vec,P,Q):
770        #        return (hv(vec,    P+Q) - hv(vec,P)-hv(vec,Q))/2
771        K = Qp(p, prec)   
772       
773        v1 = vector([K(0),K(1)])  # that is eta
774        v2 = vector([K(-1),K(1)])  # and this is eta-omega.
775        #                      the rest should not depend on this choice
776        #                      as long as it is outside Q_p * omega                   
777   
778        rk = E.rank()
779        if rk == 0:
780            return vector([K(1),K(0)])
781
782
783        basis = E.gens()
784   
785        def regv(vec):
786            M = matrix.matrix(K,rk,rk,0)
787            point_height = [hv(vec,P) for P in basis]
788            for i in range(rk):
789                for j in range(i+1, rk):
790                    M[i, j] = M[j, i] = (hv(vec,basis[i] + basis[j])- point_height[i] - point_height[j] )/2
791            for i in range(rk):
792                M[i,i] = point_height[i]
793           
794            return M.determinant()
795       
796        reg1 = regv(v1)
797        reg2 = regv(v2)
798   
799        def Dp_pairing(vec1,vec2):
800            return vec1[0]*vec2[1]-vec1[1]*vec2[0]
801   
802       
803        return (reg1 * v2 - reg2 * v1 ) / Dp_pairing(v2,v1)
Note: See TracBrowser for help on using the repository browser.