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

Fix some doctests and bugs.

Line
1"""
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#
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#
22######################################################################
23
24
25from sage.rings.integer_ring import   ZZ
26from sage.rings.rational_field import QQ
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
48    """
49    The p-adic L-series of an elliptic curve.
50
51    EXAMPLES:
52    An ordinary example:
53        sage: e = EllipticCurve('389a')
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:
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:
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:
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')
127            3-adic L-series of Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field
129            3-adic L-series of Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field (not normalized)
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:
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')
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:
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.
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:
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):
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:
268            sage: L.order_of_vanishing()
269            0
271            sage: L.order_of_vanishing()
272            0
274            sage: L.order_of_vanishing()
275            1
277            sage: L.order_of_vanishing()
278            1
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):
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:
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:
360        return None
361
362    def _set_series_in_cache(self, n, prec, f):
363        self.__series[(n,prec)] = f
364
365
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
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
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
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
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)
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:
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:
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:
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
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
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
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
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:
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')
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([+ (E.ap(p))*, - 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')
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()
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')
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')
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')
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.