source: sage/schemes/elliptic_curves/padic_lseries.py @ 5456:b89fb9f92c7d

Revision 5456:b89fb9f92c7d, 35.0 KB checked in by William Stein <wstein@…>, 6 years ago (diff)

Fix some doctests and bugs.
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, Zp
28from sage.rings.infinity import infinity
29from sage.rings.all import LaurentSeriesRing, 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, get_verbose
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        q = E0.period_lattice()[0]/E.period_lattice()[0]*E0.real_components()/E.real_components()
97        self._quotient_of_periods = QQ(gp.bestappr(q,200))
98       
99        self._modular_symbol = E.modular_symbol(sign=1, normalize=normalize)
100
101    def elliptic_curve(self):
102        """
103        Return the elliptic curve to which this p-adic L-series is associated.
104       
105        EXAMPLES:
106            sage: L = EllipticCurve('11a').padic_lseries(5)
107            sage: L.elliptic_curve()
108            Elliptic Curve defined by y^2 + y = x^3 - x^2 - 10*x - 20 over Rational Field
109        """
110        return self._E
111
112    def prime(self):
113        """
114        EXAMPLES:
115            sage: L = EllipticCurve('11a').padic_lseries(5)
116            sage: L.prime()
117            5
118        """
119        return self._p
120
121    def _repr_(self):
122        """
123        Return print representation.
124
125            sage: e = EllipticCurve('37a')
126            sage: e.padic_lseries(3)
127            3-adic L-series of Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field
128            sage: e.padic_lseries(3,normalize=False)
129            3-adic L-series of Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field (not normalized)
130            sage: L = e.padic_lseries(3,normalize=False)
131            sage: L.rename('(factor)*L_3(T)')
132            sage: L
133            (factor)*L_3(T)
134        """
135        s = "%s-adic L-series of %s"%(self._p, self._E)
136        if not self._normalize:
137            s += ' (not normalized)'
138        return s
139
140    def modular_symbol(self, r):
141        """
142        Return the modular symbol used to compute this p-adic
143        L-series evaluated at r.
144
145        EXAMPLES:
146            sage: L = EllipticCurve('11a').padic_lseries(5)
147            sage: [L.modular_symbol(r) for r in [0,1/5,oo,1/11]]
148            [1/5, 6/5, 0, 0]       
149        """
150        return self._modular_symbol(r)
151
152    def measure(self, a, n, prec):
153        r"""
154        Return the measure on $\ZZ_p^*$ defined by
155           $$
156             \mu_{E,\alpha}^+ ( a + p^n \ZZ_p  ) =
157                   \frac{1}{\alpha^n} \modsym{a}{p^n} - \frac{1}{\alpha^{n+1}} \modsym{a}{p^{n-1}}
158           $$
159        that is used to define this $p$-adic $L$-function.
160
161        INPUT:
162            a -- an integer
163            n -- a non-negative integer
164            prec -- an integer
165
166        EXAMPLES:
167            sage: E = EllipticCurve('37a')
168            sage: L = E.padic_lseries(5)
169            sage: L.measure(1,2, prec=9)
170            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)
171        """
172        try:
173            p, alpha, z, w, f = self.__measure_data[(n,prec)]
174        except (KeyError, AttributeError):
175            if not hasattr(self, '__measure_data'):
176                self.__measure_data = {}
177            p = self._p
178            alpha = self.alpha(prec=prec)
179            z = 1/(alpha**n)
180            w = p**(n-1)
181            f = self._modular_symbol
182            self.__measure_data[(n,prec)] = (p,alpha,z,w,f)
183           
184        if self._E.conductor() % p == 0:
185            return z * f(a/(p*w))
186        return z * f(a/(p*w)) - (z/alpha) * f(a/w)
187       
188    def alpha(self, prec=20):
189        r"""
190        Return a p-adic root $\alpha$ of the polynomial $x^2 - a_p x
191        + p$ with $\ord_p(\alpha) < 1$.  In the ordinary case this is
192        just the unit root.
193
194        INPUT:
195            prec -- positive integer, the p-adic precision of the root.
196
197        EXAMPLES:
198        Consider the elliptic curve 37a:
199            sage: E = EllipticCurve('37a')
200
201        An ordinary prime:
202            sage: L = E.padic_lseries(5)
203            sage: alpha = L.alpha(10); alpha
204            3 + 2*5 + 4*5^2 + 2*5^3 + 5^4 + 4*5^5 + 2*5^7 + 5^8 + 5^9 + O(5^10)
205            sage: alpha^2 - E.ap(5)*alpha + 5
206            O(5^10)
207
208        A supersingular prime.
209            sage: L = E.padic_lseries(3)
210            sage: alpha = L.alpha(10); alpha
211            (1 + O(3^10))*alpha
212            sage: alpha^2 - E.ap(3)*alpha + 3
213            (O(3^10))*alpha^2 + (O(3^11))*alpha + (O(3^11))
214
215        A reducible prime:
216            sage: L = EllipticCurve('11a').padic_lseries(5)
217            sage: L.alpha(5)
218            1 + 4*5 + 3*5^2 + 2*5^3 + 4*5^4 + O(5^5)
219        """
220        try:
221            return self._alpha[prec]
222        except AttributeError:
223            self._alpha = {}
224        except KeyError:
225            pass
226        E = self._E
227        p = self._p
228        a_p = E.ap(p)
229        K = Qp(p, prec, print_mode='series')
230 
231        if E.conductor() % p == 0:
232            self._alpha[prec] = K(a_p)
233            return K(a_p) 
234
235        R = ZZ['x']
236        f = R([p, -a_p, 1])
237        if E.is_ordinary(p):
238            G = f.factor_padic(p, prec+5)
239            for pr, e in G:
240                a = -pr[0]
241                if a.valuation() < 1:
242                    self._alpha[prec] = K(a)
243                    return K(a)
244            raise ValueError, "bug in p-adic L-function alpha"
245        else: # supersingular case
246            f = f.change_ring(Qp(p, prec, print_mode='series'))
247            a = f.root_field('alpha', check_irreducible=False).gen()
248            self._alpha[prec] = a
249            return a
250
251    def order_of_vanishing(self, proof=True):
252        """
253        Return the order of vanishing of this $p$-adic $L$-series.
254
255        The output of this function is provably correct, due to a
256        theorem of Kato.  This function will terminate if and only if
257        the Mazur-Tate-Teitelbaum analogue of the BSD conjecture about
258        the rank of the curve is true and the subgroup of elements of
259        p-power order in the Shafarevich-Tate group of this curve is
260        finite.  I.e., if this function terminates (with no errors!),
261        then you may conclude that the p-adic BSD rank conjecture is
262        true and that the p-part of Sha is finite.
263
264        NOTE: currently $p$ must be a prime of good ordinary reduction.
265
266        EXAMPLES:
267            sage: L = EllipticCurve('11a').padic_lseries(3)
268            sage: L.order_of_vanishing()
269            0
270            sage: L = EllipticCurve('11a').padic_lseries(5)                       
271            sage: L.order_of_vanishing()
272            0       
273            sage: L = EllipticCurve('37a').padic_lseries(5)
274            sage: L.order_of_vanishing()
275            1
276            sage: L = EllipticCurve('43a').padic_lseries(3)
277            sage: L.order_of_vanishing()
278            1
279            sage: L = EllipticCurve('37b').padic_lseries(3)
280            sage: L.order_of_vanishing()
281            0
282
283        We verify that Sha(E)(p) is finite for p=3,5,7 for the
284        first curve of rank 2:
285            sage: e = EllipticCurve('389a')
286            sage: for p in primes(3,10):
287            ...    print p, e.padic_lseries(p).order_of_vanishing()
288            3 2
289            5 2
290            7 2           
291        """
292        try:
293            return self.__ord
294        except AttributeError:
295            pass
296       
297        if not self.is_ordinary():
298            raise NotImplementedError
299        E = self.elliptic_curve()
300        if not E.is_good(self.prime()):
301            raise ValueError, "prime must be of good reduction"
302        r = E.rank()
303        n = 1
304        while True:
305            f = self.series(n)
306            v = f.valuation()
307            if v < r:
308                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)
309            if v == r:
310                self.__ord = v
311                return v
312            n += 1
313
314
315    def _c_bounds(self, n):
316        raise NotImplementedError
317
318    def _prec_bounds(self, n,prec):
319        raise NotImplementedError
320
321    def teichmuller(self, prec):
322        r"""
323        Return Teichmuller lifts to the given precision.
324       
325        INPUT:
326            prec -- a positive integer.
327           
328        OUTPUT:
329            the cached Teichmuller lifts
330
331        EXAMPLES:
332            sage: L = EllipticCurve('11a').padic_lseries(7)
333            sage: L.teichmuller(1)
334            [0, 1, 2, 3, 4, 5, 6]
335            sage: L.teichmuller(2)
336            [0, 1, 30, 31, 18, 19, 48]
337        """
338        p = self._p
339        K = Qp(p, prec, print_mode='series')       
340        return [Integer(0)] + \
341               [a.residue(prec).lift() for a in K.teichmuller_system()]
342
343    def _e_bounds(self, n, prec):
344        p = self._p
345        prec = max(2,prec)
346        R = PowerSeriesRing(ZZ,'T',prec+1)
347        T = R(R.gen(),prec +1)
348        w = (1+T)**(p**n) - 1
349        return [infinity] + [valuation(w[j],p) for j in range(1,min(w.degree()+1,prec))]
350       
351    def _get_series_from_cache(self, n, prec):
352        try:
353            return self.__series[(n,prec)]
354        except AttributeError:
355            self.__series = {}
356        except KeyError:
357            for _n, _prec in self.__series.keys():
358                if _n == n and _prec <= prec:
359                    return self.__series[(_n,_prec)].add_bigoh(prec)
360        return None
361
362    def _set_series_in_cache(self, n, prec, f):
363        self.__series[(n,prec)] = f
364       
365
366class pAdicLseriesOrdinary(pAdicLseries):
367    def series(self, n=2, prec=5):
368        r"""
369        Return the $n$-th approximation to the $p$-adic $L$-series as
370        a power series in $T$ (corresponding to $\gamma-1$ with
371        $\gamma=1+p$ as a generator of $1+p\mathbb{Z}_p$).  Each
372        coefficient is a $p$-adic number whose precision is provably
373        correct.
374
375        INPUT:
376            n -- (default: 2) a positive integer
377            prec -- (default: 5) maxima number of terms of the series
378                    to compute; to compute as many as possible just
379                    give a very large number for prec; the result will
380                    still be correct.
381
382        ALIAS:
383            power_series is identical to series.
384
385        EXAMPLES:
386        We compute some $p$-adic $L$-functions associated to the elliptic
387        curve 11a.
388       
389            sage: E = EllipticCurve('11a')
390            sage: p = 3
391            sage: E.is_ordinary(p)
392            True
393            sage: L = E.padic_lseries(p)
394            sage: L.series(3)
395            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)
396
397
398        Another example at a prime of bad reduction, where the
399        $p$-adic $L$-function has an extra 0 (compared to the non
400        $p$-adic $L$-function).
401       
402            sage: E = EllipticCurve('11a')
403            sage: p = 11
404            sage: E.is_ordinary(p)
405            True
406            sage: L = E.padic_lseries(p)
407            sage: L.series(2)
408            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)
409           
410        We compute a $p$-adic $L$-function that vanishes to order $2$.
411            sage: E = EllipticCurve('389a')
412            sage: p = 3
413            sage: E.is_ordinary(p)
414            True
415            sage: L = E.padic_lseries(p)
416            sage: L.series(1)
417            O(T^1)
418            sage: L.series(2)
419            O(3^4) + O(3^1)*T + (2 + O(3))*T^2 + O(T^3)
420            sage: L.series(3)
421            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)
422        """
423        n = ZZ(n)
424        if n < 1:
425            raise ValueError, "n (=%s) must be a positive integer"%n
426       
427
428        p = self._p
429        #verbose("computing L-series for p=%s, n=%s, and prec=%s"%(p,n,prec))
430       
431        bounds = self._prec_bounds(n,prec)
432        padic_prec = max(bounds[1:]) + 5
433        verbose("using p-adic precision of %s"%padic_prec)
434       
435        res_series_prec = min(p**(n-1), prec)
436        verbose("using series precision of %s"%res_series_prec)
437       
438        ans = self._get_series_from_cache(n, res_series_prec)
439        if not ans is None:
440            verbose("found series in cache")
441            return ans
442 
443        K = QQ
444        gamma = K(1 + p)
445        R = PowerSeriesRing(K,'T',res_series_prec)
446        T = R(R.gen(),res_series_prec )
447        L = R(0) 
448        one_plus_T_factor = R(1) 
449        gamma_power = K(1)
450        teich = self.teichmuller(padic_prec)
451        p_power = p**(n-1)
452
453        verbose("Now iterating over %s summands"%((p-1)*p_power))
454        verbose_level = get_verbose()
455        count_verb = 0
456        for j in range(p_power):
457            s = K(0)
458            if verbose_level >= 2 and j/p_power*100 > count_verb + 3:
459                verbose("%.2f percent done"%(float(j)/p_power*100))
460                count_verb += 3
461            for a in range(1,p):
462                b = teich[a] * gamma_power
463                s += self.measure(b, n, padic_prec).lift()
464            L += s * one_plus_T_factor
465            one_plus_T_factor *= 1+T
466            gamma_power *= gamma
467
468        # Now create series but with each coefficient truncated
469        # so it is proven correct:
470        K = Qp(p, padic_prec, print_mode='series')
471        R = PowerSeriesRing(K,'T',res_series_prec)
472        L = R(L,res_series_prec)
473        aj = L.list()
474        if len(aj) > 0:
475            aj = [aj[0].add_bigoh(padic_prec-2)] + [aj[j].add_bigoh(bounds[j]) for j in range(1,len(aj))]
476        L = R(aj,res_series_prec ) * self._quotient_of_periods
477       
478        self._set_series_in_cache(n, res_series_prec, L)
479       
480        return L
481
482    power_series = series
483
484    def is_ordinary(self):
485        """
486        Return True if the elliptic that this $L$-function is attached
487        to is ordinary.
488       
489        EXAMPLES:
490            sage: L = EllipticCurve('11a').padic_lseries(5)           
491            sage: L.is_ordinary()
492            True
493        """
494        return True
495
496    def is_supersingular(self):
497        """
498        Return True if the elliptic that this L function is attached
499        to is supersingular.
500               
501        EXAMPLES:
502            sage: L = EllipticCurve('11a').padic_lseries(5)           
503            sage: L.is_supersingular()
504            False
505        """
506        return False
507
508    def _c_bound(self):
509        try:
510            return self.__c_bound
511        except AttributeError:
512            pass
513        E = self._E
514        p = self._p
515        if E.is_irreducible(p):
516            ans = 0
517        else:
518            m = E.modular_symbol_space(sign=1)
519            b = m.boundary_map().codomain()
520            C = b._known_cusps()  # all known, since computed the boundary map
521            ans = max([valuation(self.modular_symbol(a).denominator(), p) for a in C]) 
522        self.__c_bound = ans
523        return ans
524               
525    def _prec_bounds(self, n, prec):
526        p = self._p
527        e = self._e_bounds(n-1, prec)
528        c = self._c_bound()
529        return [e[j] - c for j in range(len(e))]
530
531
532class pAdicLseriesSupersingular(pAdicLseries):
533    def series(self, n=3, prec=5):
534        r"""
535        Return the $n$-th approximation to the $p$-adic $L$-series as a
536        power series in $T$ (corresponding to $\gamma-1$ with
537        $\gamma=1+p$ as a generator of $1+p\mathbb{Z}_p$).  Each
538        coefficient is a $p$-adic number whose precision is probably
539        correct.
540
541        INPUT:
542            n -- (default: 3) a positive integer
543            prec -- (default: 5) maxima number of terms of the series
544                    to compute; to compute as many as possible just
545                    give a very large number for prec; the result will
546                    still be correct.
547
548        ALIAS:
549            power_series is identical to series.
550
551        EXAMPLES:
552        A superingular example, where we must compute to higher precision to see anything.
553            sage: e = EllipticCurve('37a')
554            sage: L = e.padic_lseries(3); L
555            3-adic L-series of Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field
556            sage: L.series(2)
557            O(T^3)
558            sage: L.series(4)         # takes a long time (several seconds)
559            (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)
560            sage: L.alpha(2).parent()
561            Univariate Quotient Polynomial Ring in alpha over 3-adic Field with capped
562            relative precision 2 with modulus (1 + O(3^2))*x^2 + (3 + O(3^3))*x + (3 + O(3^3))
563        """
564        n = ZZ(n)
565        if n < 1:
566            raise ValueError, "n (=%s) must be a positive integer"%n
567       
568        p = self._p
569        prec = min(p**(n-1), prec)
570        bounds = self._prec_bounds(n,prec)
571        padic_prec = max(sum(bounds[1:],[])) + 5
572        verbose("using p-adic precision of %s"%padic_prec) 
573        ans = self._get_series_from_cache(n, prec)
574        if not ans is None:
575            verbose("found series in cache")
576            return ans
577
578        alpha = self.alpha(prec=padic_prec)
579        K = alpha.parent()
580        gamma = 1 + p
581        R = PowerSeriesRing(K,'T',prec)
582        T = R(R.gen(), prec)
583        L = R(0)
584        one_plus_T_factor = R(1)
585        gamma_power = 1
586        teich = self.teichmuller(padic_prec)
587
588        verbose("Now iterating over %s summands"%((p-1)*p**(n-1)))
589        verbose_level = get_verbose()
590        count_verb = 0
591        for j in range(p**(n-1)):
592            s = K(0)
593            if verbose_level >= 2 and j/p**(n-1)*100 > count_verb + 3:
594                verbose("%.2f percent done"%(float(j)/p**(n-1)*100))
595                count_verb += 3
596            for a in range(1,p):
597                b = teich[a] * gamma_power
598                s += self.measure(b, n, padic_prec)
599            L += s * one_plus_T_factor
600            one_plus_T_factor *= 1+T
601            gamma_power *= gamma
602
603        # Now create series but with each coefficient truncated
604        # so it is proven correct:
605        L = R(L,prec) 
606        aj = L.list()
607        if len(aj) > 0:
608            bj = [aj[0][0].add_bigoh(padic_prec-2) + alpha * aj[0][1].add_bigoh(padic_prec-2)]
609            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))]
610            L = R(bj, prec)
611        L = L * self._quotient_of_periods
612        self._set_series_in_cache(n, prec, L)
613        return L
614
615    power_series = series
616
617    def is_ordinary(self):
618        return False
619
620    def is_supersingular(self):
621        return True
622
623    def _prec_bounds(self, n,prec):
624        p = self._p
625        e = self._e_bounds(n-1,prec)
626        c0 = ZZ(n+2)/2
627        c1 = ZZ(n+3)/2
628        return [[infinity,infinity]] + [[(e[j] - c0).floor(), (e[j] - c1).floor()] for j in range(1,len(e))]
629
630
631    def Dp_valued_series(self, n=3, prec=5):
632        r"""
633        Returns a vector of two components which are p-adic power series.
634        The answer v is such that
635            $$(1-\varphi)^(-2)* L_p(E,T) = v[1] * \omega + v[2] * \varphi(\omega)$$
636        as an element of the Dieudonne module $D_p(E) = H^1_{dR}(E/\QQ_p)$ where
637        $\omega$ is the invariant differential and $\varphi$ is the Frobenius on $D_p(E)$.   
638        According to the p-adic BSD this function has a zero of order
639        rank(E(Q)) and it's leading term is
640        \begin{verbatim}
641           +- #Sha(E/Q) * Tamagawa product / Torsion^2 * padic height regulator with values in D_p(E).
642        \end{verbatim}
643   
644       
645        INPUT:
646            n -- (default: 3) a positive integer
647            prec -- (default: 5) a positive integer
648
649        EXAMPLES:
650            sage: E = EllipticCurve('14a')
651            sage: L = E.padic_lseries(5)
652            sage: L.Dp_valued_series(4)
653            (4 + 4*5^2 + O(5^4) + (1 + O(5))*T + (4 + O(5))*T^2 + (1 + O(5))*T^3 + (3 + O(5))*T^4 + O(T^5), O(5^4) + O(5^1)*T + O(5^1)*T^2 + O(5^1)*T^3 + (3 + O(5))*T^4 + O(T^5))
654        """
655        E = self._E
656        p = self._p
657        lps = self.series(n, prec=prec)
658   
659        # now split up the series in two lps = G + H * alpha
660        R = lps.base_ring().base_ring() # Qp
661        QpT , T = PowerSeriesRing(R,'T',prec).objgen()
662        G = QpT([lps[n][0] for n in range(0,lps.prec())], prec)
663        H = QpT([lps[n][1] for n in range(0,lps.prec())], prec) 
664         
665        # now compute phi
666        phi = matrix.matrix([[0,-1/p],[1,E.ap(p)/p]])
667        lpv = vector([G  + (E.ap(p))*H  , - R(p) * H ])  # this is L_p
668        eps = (1-phi)**(-2) 
669        resu = lpv*eps.transpose()
670        return resu
671
672
673    def frobenius(self, prec=20, method = "mw"):
674        r"""
675        This returns a geometric Frobenius $\varphi$ on the Diedonne module $D_p(E)$
676        with respect to the basis $\omega$, the invariant differential and $\eta=x\omega$.
677        It satisfies  $phi^2 - a_p/p*phi + 1/p = 0$.
678       
679        INPUT:
680            prec -- (default: 20) a positive integer
681            method -- either "mw" (default) for Monsky-Washintzer
682                      or "approx" for the method describedby Bernardi and Perrin-Riou
683                      (much slower)
684                         
685
686        EXAMPLES:
687            sage: E = EllipticCurve('14a')
688            sage: L = E.padic_lseries(5)
689            sage: phi = L.frobenius(5)
690            sage: phi
691            [                  2 + 5^2 + 5^4 + O(5^5)    3*5^-1 + 3 + 5 + 4*5^2 + 5^3 + O(5^4)]
692            [      3 + 3*5^2 + 4*5^3 + 3*5^4 + O(5^5) 3 + 4*5 + 3*5^2 + 4*5^3 + 3*5^4 + O(5^5)]
693            sage: -phi^2
694            [5^-1 + O(5^4)        O(5^4)]
695            [       O(5^5) 5^-1 + O(5^4)]
696        """
697        E = self._E
698        p = self._p
699        if method != "mw" and method !="approx":
700            raise ValueError, "Unknown method %s."%method
701        if method == "approx":
702            return self.__phi_bpr(prec=prec)
703        if p < 4 and method == "mw":
704            print "Warning: If this fails try again using method=\"approx\""
705        Ew = E.integral_weierstrass_model()
706        adjusted_prec = monsky_washnitzer.adjusted_prec(p, prec)
707        modprecring = Integers(p**adjusted_prec)
708        output_ring = Qp(p, prec)
709        R, x = PolynomialRing(modprecring, 'x').objgen()
710        Q = x**3 + modprecring(Ew.a4()) * x + modprecring(Ew.a6())
711        trace = Ew.ap(p)
712        fr = monsky_washnitzer.matrix_of_frobenius(Q, p, adjusted_prec, trace)
713        fr = matrix.matrix(output_ring,2,2,fr)
714       
715        # return a vector for pari's ellchangecurve to pass from e1 to e2
716        def isom(e1,e2):
717            if not e1.is_isomorphic(e2):
718                raise ValueError, "Curves must be isomorphic."
719            usq = (e1.discriminant()/e2.discriminant()).nth_root(6)
720            u = usq.sqrt()
721            s = (u   *  e2.a1() - e1.a1() )/ZZ(2)
722            r = (usq *  e2.a2() - e1.a2() + s**2 + e1.a1()*s)/ZZ(3)
723            t = (u**3 * e2.a3() - e1.a3() - e1.a1()*r)/ZZ(2)
724            return [u,r,s,t]
725
726        v = isom(E,Ew)
727        u = v[0]
728        r = v[1]
729       
730        # change basis
731        A = matrix.matrix([[u,-r/u],[0,1/u]])
732        frn = A * fr * A**(-1) 
733        return 1/p*frn
734
735
736    # returns the phi using the definition of bernardi-perrin-riou on page 232.
737   
738    def __phi_bpr(self, prec=0):
739        E = self._E
740        p = self._p
741        if prec > 10:
742            print "Warning: Very large value for the precision." 
743        if prec == 0:
744            prec = floor((log(10000)/log(p)))
745            verbose("prec set to %s"%prec)
746        eh = E.formal()
747        om = eh.differential(prec = p**prec+3)
748        verbose("differential computed")
749        xt = eh.x(prec=p**prec + 3)
750        et = xt*om
751        # c_(p^k) = cs[k] d...   
752        cs = [om[p**k-1] for k in range(0,prec+1)]
753        ds = [et[p**k-1] for k in range(0,prec+1)]
754        delta = 0
755        dpr = 0
756        gamma = 0
757        dga = 0
758        for k in range(1,prec+1):
759            # this is the equation eq[0]*x+eq[1]*y+eq[2] == 0
760            # such that delta_ = delta + d^dpr*x ...
761            eq = [(p**dpr*cs[k]) % p**k,(-p**dga*ds[k]) % p**k , (delta*cs[k]-gamma*ds[k]-cs[k-1]) % p**k ] 
762            verbose("valuations : %s"%([x.valuation(p) for x in eq]))
763            v = min([x.valuation(p) for x in eq])
764            if v == infinity:
765                verbose("no new information at step k=%s"%k)
766            else:
767                eq = [ZZ(x/p**v) for x in eq]
768                verbose("renormalised eq mod p^%s is now %s"%(k-v,eq))
769                if eq[0].valuation(p) == 0:
770                    l = min(eq[1].valuation(p),k-v)
771                    if l == 0:
772                        verbose("not uniquely determined at step k=%s"%k)
773                    else: 
774                        ainv = eq[0].inverse_mod(p**l)
775                        delta = delta - eq[2]*ainv*p**dpr
776                        dpr = dpr + l
777                        delta = delta % p**dpr
778                        verbose("delta_prec increased to %s\n delta is now %s"%(dpr,delta))
779                elif eq[1].valuation(p) == 0:
780                    l = min(eq[0].valuation(p),k-v)
781                    ainv = eq[1].inverse_mod(p**l)
782                    gamma = gamma - eq[2]*ainv*p**dga
783                    dga = dga + l
784                    gamma = gamma % p**dga
785                    verbose("gamma_prec increased to %s\n gamma is now %s"%(dga,gamma))
786                else:
787                    raise RuntimeError,  "Bug: no delta or gamma can exist"
788   
789        # end of approximation of delta and gamma
790        delta = Qp(p,dpr)(delta)
791        gamma = Qp(p,dga)(gamma)
792        verbose("result delta = %s\n      gamma = %s\n check : %s"%(delta,gamma, [Qp(3,k)(delta * cs[k] - gamma * ds[k] - cs[k-1]) for k in range(1,prec+1)] ))
793        a = delta
794        c = -gamma
795        d = E.ap(p) - a
796        b = (-1/p+a*d)/c
797        phi = matrix.matrix([[a,b],[c,d]])
798        return phi
799
800
801    def bernardi_sigma_function(self, prec=20):
802        r"""
803        Return the  p-adic sigma function of Bernardi in terms of $z = log(t)$.
804        This is the same as padic_sigma with E2 = 0.
805
806        EXAMPLES:
807            sage: E = EllipticCurve('14a')
808            sage: L = E.padic_lseries(5)
809            sage: L.bernardi_sigma_function(5) # Todo: some sort of consistency check!?
810            z + 1/24*z^3 + 29/384*z^5 - 8399/322560*z^7 - 291743/92897280*z^9 - 4364831/5225472*z^10 + 2172371753/955514880*z^11 - 17875714529/6897623040*z^12 + 2839176621047/1605264998400*z^13 + 32012675789849/10042939146240*z^14 - 367444910151047/89894839910400*z^15 + 973773806885959/241030539509760*z^16 - 33997971208432501/17259809262796800*z^17 - 10331978660756704339/842918229599846400*z^18 + 18601407947897364480389/950670294194847744000*z^19 - 118837570440101901119321/8071784966648129126400*z^20 + O(z^21)
811        """
812        E = self._E
813        p = self._p
814       
815        Eh = E.formal()
816        lo = Eh.log(prec + 5)
817        F = lo.reversion()
818       
819        S = LaurentSeriesRing(QQ,'z')
820        z = S.gen()
821        F = F(z)
822        xofF = Eh.x(prec + 2)(F)
823        #r =  ( E.a1()**2 + 4*E.a2() ) / ZZ(12)
824        g = (1/z**2 - xofF ).power_series()
825        h = g.integral().integral()
826        sigma_of_z = z.power_series() * h.exp()
827
828        return sigma_of_z
829 
830
831    def Dp_valued_height(self,prec=20):
832        r"""
833        Returns the canonical $p$-adic height with values in the Dieudonne module $D_p(E)$.
834        It is defined to be
835            $$h_{\eta} \cdot \omega - h_{\omega} \cdot \eta$$
836        where $h_{\eta}$ is made out of the sigma function of Bernardi and
837        $h_{\omega}$ is $-log^2$.
838        The answer $v$ is given as $v[1]*omega + v[2]*eta$.
839        The coordinates of $v$ are dependent of the
840        Weierstrass equation.
841 
842        EXAMPLES:
843            sage: E = EllipticCurve('53a')
844            sage: L = E.padic_lseries(5)
845            sage: h = L.Dp_valued_height(7)
846            sage: h(E.gens()[0])
847            (2*5 + 3*5^2 + 2*5^3 + 5^4 + 3*5^6 + 3*5^7 + O(5^8),
848            4*5^2 + 4*5^3 + 4*5^5 + 4*5^6 + 2*5^7 + 5^8 + O(5^9))
849        """
850        E = self._E
851        p = self._p
852        Ehat = E.formal()
853        elog = Ehat.log(prec + Integer(3))
854   
855        # we will have to do it properly with David Harvey's _DivPolyContext(E, R, Q)
856        n = arith.LCM(E.tamagawa_numbers())
857        n = arith.LCM(n, E.Np(p)) # allowed here because E has good reduction at p
858   
859        if p < 5:
860            phi = self.frobenius(min(6,prec),method="approx")
861        else:
862            phi = self.frobenius(prec+2,method="mw")
863           
864        def height(P,check=True):
865            if P.is_finite_order():
866                return Qp(p,prec)(0)
867            if check:
868                assert P.curve() == E, "the point P must lie on the curve from which the height function was created"
869            Q = n * P
870            tt = - Q[0]/Q[1]
871            R = Qp(p,prec+5)
872            tt = R(tt)
873            zz = elog(tt)
874
875            homega = -zz**2/n**2
876       
877            eQ = denominator(Q[1])/denominator(Q[0])
878            si = self.bernardi_sigma_function(prec=prec+4)
879            heta =  2 * log(si(zz)/eQ) / n**2 
880       
881            R = Qp(p,prec)
882               
883            return vector([R(heta),-R(homega)])
884
885        return height
886   
887   
888   
889    def Dp_valued_regulator(self,prec=20,v1=0,v2=0):
890        """
891        Returns the canonical $p$-adic regulator with values in the Dieudonne module $D_p(E)$
892        as defined by Perrin-Riou using the $p$-adic height with values in $D_p(E)$.
893        The result is written in the basis $\omega$, $\varphi(\omega)$, and hence the
894        coordinates of the result are independent of the chosen Weierstrass equation.
895       
896        NOTE:
897            The definition here is corrected with repect to Perrin-Riou's article
898            'Arithm\'etique des courbes elliptiques \`a r\'eduction supersinguli\`ere en $p$'.
899       
900   
901        EXAMPLES:
902            sage: E = EllipticCurve('43a')
903            sage: L = E.padic_lseries(7)
904            sage: L.Dp_valued_regulator(7)
905            (2*7 + 2*7^3 + 2*7^4 + 5*7^5 + 6*7^6 + 2*7^7 + O(7^8), 3*7^2 + 4*7^3 + 3*7^4 + 5*7^5 + 2*7^7 + O(7^8))
906        """
907   
908        p = self._p
909        E = self._E
910       
911        h =  self.Dp_valued_height(prec=prec)
912   
913        # this is the height_{v} (P) for a v in D_p
914        def hv(vec,P):
915            hP = h(P)
916            return - vec[0]*hP[1] +vec[1]*hP[0]
917       
918        #    def hvpairing(vec,P,Q):
919        #        return (hv(vec,    P+Q) - hv(vec,P)-hv(vec,Q))/2
920        K = Qp(p, prec)   
921       
922        if v1 ==0 and v2 ==0 :
923            v1 = vector([K(0),K(1)])  # that is eta
924            v2 = vector([K(-1),K(1)])  # and this is eta-omega.
925        #                      the rest should not depend on this choice
926        #                      as long as it is outside Q_p * omega                   
927   
928        rk = E.rank()
929        if rk == 0:
930            return vector([K(1),K(0)])
931
932
933        basis = E.gens()
934   
935        def regv(vec):
936            M = matrix.matrix(K,rk,rk,0)
937            point_height = [hv(vec,P) for P in basis]
938            for i in range(rk):
939                for j in range(i+1, rk):
940                    M[i, j] = M[j, i] = (hv(vec,basis[i] + basis[j])- point_height[i] - point_height[j] )/2
941            for i in range(rk):
942                M[i,i] = point_height[i]
943           
944            return M.determinant()
945       
946       
947        def Dp_pairing(vec1,vec2):
948            return (vec1[0]*vec2[1]-vec1[1]*vec2[0])
949   
950        omega_vec = vector([K(1),K(0)])
951       
952        # note the correction here with respect to Perrin-Riou's definition.
953        # only this way the result will be indep of the choice of v1 and v2.
954        reg1 = regv(v1)/Dp_pairing(omega_vec,v1)**(rk-1)
955       
956        reg2 = regv(v2)/Dp_pairing(omega_vec,v2)**(rk-1)
957   
958       
959        # the regulator in the basis omega,eta
960        reg_oe = (reg1 * v2 - reg2 * v1 ) / Dp_pairing(v2,v1)
961 
962        if p < 5:
963            phi = self.frobenius(min(6,prec),method="approx")
964        else:
965            phi = self.frobenius(prec+2,method="mw")
966           
967        c = phi[1,0]  # this is the 'period' [omega,phi(omega)]
968        a = phi[0,0]
969 
970        return vector([reg_oe[0] - a/c*reg_oe[1],reg_oe[1]/c])
Note: See TracBrowser for help on using the repository browser.