Coverage for local/lib/python2.7/site-packages/sage/schemes/elliptic_curves/padics.py : 83%

# -*- coding: utf-8 -*- 

""" 

Miscellaneous `p`-adic functions 

`p`-adic functions from ell_rational_field.py, moved here to reduce 

crowding in that file. 

""" 

###################################################################### 

# Copyright (C) 2007 William Stein <wstein@gmail.com> 

# 

# Distributed under the terms of the GNU General Public License (GPL) 

# 

# This code is distributed in the hope that it will be useful, 

# but WITHOUT ANY WARRANTY; without even the implied warranty of 

# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 

# General Public License for more details. 

# 

# The full text of the GPL is available at: 

# 

# http://www.gnu.org/licenses/ 

###################################################################### 

from __future__ import absolute_import 

import sage.rings.all as rings 

from . import padic_lseries as plseries 

import sage.arith.all as arith 

from sage.rings.all import ( 

Qp, Zp, 

Integers, 

Integer, 

O, 

PowerSeriesRing, 

LaurentSeriesRing, 

RationalField) 

import math 

import sage.misc.misc as misc 

import sage.matrix.all as matrix 

sqrt = math.sqrt 

import sage.schemes.hyperelliptic_curves.monsky_washnitzer 

import sage.schemes.hyperelliptic_curves.hypellfrob 

from sage.misc.all import cached_method 

def __check_padic_hypotheses(self, p): 

r""" 

Helper function that determines if `p` 

is an odd prime of good ordinary reduction. 

EXAMPLES:: 

sage: E = EllipticCurve('11a1') 

sage: from sage.schemes.elliptic_curves.padics import __check_padic_hypotheses 

sage: __check_padic_hypotheses(E,5) 

5 

sage: __check_padic_hypotheses(E,29) 

Traceback (most recent call last): 

... 

ArithmeticError: p must be a good ordinary prime 

""" 

p = rings.Integer(p) 

if not p.is_prime(): 

raise ValueError("p = (%s) must be prime"%p) 

if p == 2: 

raise ValueError("p must be odd") 

if self.conductor() % p == 0 or self.ap(p) % p == 0: 

raise ArithmeticError("p must be a good ordinary prime") 

return p 

def _normalize_padic_lseries(self, p, normalize, use_eclib, implementation, precision): 

r""" 

Normalize parameters for :meth:`padic_lseries`. 

TESTS:: 

sage: from sage.schemes.elliptic_curves.padics import _normalize_padic_lseries 

sage: u = _normalize_padic_lseries(None, 5, None, None, 'sage', 10) 

sage: v = _normalize_padic_lseries(None, 5, "L_ratio", None, 'sage', 10) 

sage: u == v 

True 

""" 

if use_eclib is not None: 

from sage.misc.superseded import deprecation 

deprecation(812,"Use the option 'implementation' instead of 'use_eclib'") 

if implementation == 'pollackstevens': 

raise ValueError 

if use_eclib: 

implementation = 'eclib' 

else: 

implementation = 'sage' 

if implementation == 'eclib': 

if normalize is None: 

normalize = "L_ratio" 

elif implementation == 'sage': 

if normalize is None: 

normalize = "L_ratio" 

elif implementation == 'pollackstevens': 

if precision is None: 

raise ValueError("Must specify precision when using 'pollackstevens'") 

if normalize is not None: 

raise ValueError("The 'normalize' parameter is not used for Pollack-Stevens' overconvergent modular symbols") 

else: 

raise ValueError("Implementation should be one of 'sage', 'eclib' or 'pollackstevens'") 

#if precision is not None and implementation != 'pollackstevens': 

# raise ValueError("Must *not* specify precision unless using 'pollackstevens'") 

return (p, normalize, implementation, precision) 

@cached_method(key=_normalize_padic_lseries) 

def padic_lseries(self, p, normalize = None, use_eclib = None, implementation = 'eclib', precision = None): 

r""" 

Return the `p`-adic `L`-series of self at 

`p`, which is an object whose approx method computes 

approximation to the true `p`-adic `L`-series to 

any desired precision. 

INPUT: 

- ``p`` - prime 

- ``normalize`` - 'L_ratio' (default), 'period' or 'none'; 

this is describes the way the modular symbols 

are normalized. See modular_symbol for 

more details. 

- ``use_eclib`` - deprecated, use ``implementation`` instead 

- ``implementation`` - 'eclib' (default), 'sage', 'pollackstevens'; 

Whether to use John Cremona's eclib, the Sage implementation, 

or Pollack-Stevens' implementation of overconvergent 

modular symbols. 

EXAMPLES:: 

sage: E = EllipticCurve('37a') 

sage: L = E.padic_lseries(5); L 

5-adic L-series of Elliptic Curve defined by y^2 + y = x^3 - x over Rational Field 

sage: type(L) 

<class 'sage.schemes.elliptic_curves.padic_lseries.pAdicLseriesOrdinary'> 

We compute the `3`-adic `L`-series of two curves of 

rank `0` and in each case verify the interpolation property 

for their leading coefficient (i.e., value at 0):: 

sage: e = EllipticCurve('11a') 

sage: ms = e.modular_symbol() 

sage: [ms(1/11), ms(1/3), ms(0), ms(oo)] 

[0, -3/10, 1/5, 0] 

sage: ms(0) 

1/5 

sage: L = e.padic_lseries(3) 

sage: P = L.series(5) 

sage: P(0) 

2 + 3 + 3^2 + 2*3^3 + 2*3^5 + 3^6 + O(3^7) 

sage: alpha = L.alpha(9); alpha 

2 + 3^2 + 2*3^3 + 2*3^4 + 2*3^6 + 3^8 + O(3^9) 

sage: R.<x> = QQ[] 

sage: f = x^2 - e.ap(3)*x + 3 

sage: f(alpha) 

O(3^9) 

sage: r = e.lseries().L_ratio(); r 

1/5 

sage: (1 - alpha^(-1))^2 * r 

2 + 3 + 3^2 + 2*3^3 + 2*3^5 + 3^6 + 3^7 + O(3^9) 

sage: P(0) 

2 + 3 + 3^2 + 2*3^3 + 2*3^5 + 3^6 + O(3^7) 

Next consider the curve 37b:: 

sage: e = EllipticCurve('37b') 

sage: L = e.padic_lseries(3) 

sage: P = L.series(5) 

sage: alpha = L.alpha(9); alpha 

1 + 2*3 + 3^2 + 2*3^5 + 2*3^7 + 3^8 + O(3^9) 

sage: r = e.lseries().L_ratio(); r 

1/3 

sage: (1 - alpha^(-1))^2 * r 

3 + 3^2 + 2*3^4 + 2*3^5 + 2*3^6 + 3^7 + O(3^9) 

sage: P(0) 

3 + 3^2 + 2*3^4 + 2*3^5 + O(3^6) 

We can use Sage modular symbols instead to compute the `L`-series:: 

sage: e = EllipticCurve('11a') 

sage: L = e.padic_lseries(3, implementation = 'sage') 

sage: L.series(5,prec=10) 

2 + 3 + 3^2 + 2*3^3 + 2*3^5 + 3^6 + O(3^7) + (1 + 3 + 2*3^2 + 3^3 + O(3^4))*T + (1 + 2*3 + O(3^4))*T^2 + (3 + 2*3^2 + O(3^3))*T^3 + (2*3 + 3^2 + O(3^3))*T^4 + (2 + 2*3 + 2*3^2 + O(3^3))*T^5 + (1 + 3^2 + O(3^3))*T^6 + (2 + 3^2 + O(3^3))*T^7 + (2 + 2*3 + 2*3^2 + O(3^3))*T^8 + (2 + O(3^2))*T^9 + O(T^10) 

Finally, we can use the overconvergent method of Pollack-Stevens.:: 

sage: e = EllipticCurve('11a') 

sage: L = e.padic_lseries(3, implementation = 'pollackstevens', precision = 6) 

sage: L.series(5) 

2 + 3 + 3^2 + 2*3^3 + 2*3^5 + O(3^6) + (1 + 3 + 2*3^2 + 3^3 + O(3^4))*T + (1 + 2*3 + O(3^2))*T^2 + (3 + O(3^2))*T^3 + O(3^0)*T^4 + O(T^5) 

sage: L[3] 

3 + O(3^2) 

Another example with a semistable prime.:: 

sage: E = EllipticCurve("11a1") 

sage: L = E.padic_lseries(11, implementation = 'pollackstevens', precision=3) 

sage: L[1] 

10 + 3*11 + O(11^2) 

sage: L[3] 

O(11^0) 

""" 

p, normalize, implementation, precision = self._normalize_padic_lseries(p,\ 

normalize, use_eclib, implementation, precision) 

if implementation in ['sage', 'eclib']: 

if self.ap(p) % p != 0: 

Lp = plseries.pAdicLseriesOrdinary(self, p, 

normalize = normalize, implementation = implementation) 

else: 

Lp = plseries.pAdicLseriesSupersingular(self, p, 

normalize = normalize, implementation = implementation) 

else: 

phi = self.pollack_stevens_modular_symbol(sign=0) 

if phi.parent().level() % p == 0: 

Phi = phi.lift(p, precision, eigensymbol = True) 

else: 

Phi = phi.p_stabilize_and_lift(p, precision, eigensymbol = True) 

Lp = Phi.padic_lseries() #mm TODO should this pass precision on too ? 

Lp._cinf = self.real_components() 

return Lp 

def padic_regulator(self, p, prec=20, height=None, check_hypotheses=True): 

r""" 

Computes the cyclotomic `p`-adic regulator of this curve. 

INPUT: 

- ``p`` - prime = 5 

- ``prec`` - answer will be returned modulo 

`p^{\mathrm{prec}}` 

- ``height`` - precomputed height function. If not 

supplied, this function will call padic_height to compute it. 

- ``check_hypotheses`` - boolean, whether to check 

that this is a curve for which the p-adic height makes sense 

OUTPUT: The p-adic cyclotomic regulator of this curve, to the 

requested precision. 

If the rank is 0, we output 1. 

.. TODO:: 

Remove restriction that curve must be in minimal 

Weierstrass form. This is currently required for E.gens(). 

AUTHORS: 

- Liang Xiao: original implementation at the 2006 MSRI 

graduate workshop on modular forms 

- David Harvey (2006-09-13): cleaned up and integrated into Sage, 

removed some redundant height computations 

- Chris Wuthrich (2007-05-22): added multiplicative and 

supersingular cases 

- David Harvey (2007-09-20): fixed some precision loss that was 

occurring 

EXAMPLES:: 

sage: E = EllipticCurve("37a") 

sage: E.padic_regulator(5, 10) 

5 + 5^2 + 5^3 + 3*5^6 + 4*5^7 + 5^9 + O(5^10) 

An anomalous case:: 

sage: E.padic_regulator(53, 10) 

26*53^-1 + 30 + 20*53 + 47*53^2 + 10*53^3 + 32*53^4 + 9*53^5 + 22*53^6 + 35*53^7 + 30*53^8 + O(53^9) 

An anomalous case where the precision drops some:: 

sage: E = EllipticCurve("5077a") 

sage: E.padic_regulator(5, 10) 

5 + 5^2 + 4*5^3 + 2*5^4 + 2*5^5 + 2*5^6 + 4*5^7 + 2*5^8 + 5^9 + O(5^10) 

Check that answers agree over a range of precisions:: 

sage: max_prec = 30 # make sure we get past p^2 # long time 

sage: full = E.padic_regulator(5, max_prec) # long time 

sage: for prec in range(1, max_prec): # long time 

....: assert E.padic_regulator(5, prec) == full # long time 

A case where the generator belongs to the formal group already 

(:trac:`3632`):: 

sage: E = EllipticCurve([37,0]) 

sage: E.padic_regulator(5,10) 

2*5^2 + 2*5^3 + 5^4 + 5^5 + 4*5^6 + 3*5^8 + 4*5^9 + O(5^10) 

The result is not dependent on the model for the curve:: 

sage: E = EllipticCurve([0,0,0,0,2^12*17]) 

sage: Em = E.minimal_model() 

sage: E.padic_regulator(7) == Em.padic_regulator(7) 

True 

Allow a Python int as input:: 

sage: E = EllipticCurve('37a') 

sage: E.padic_regulator(int(5)) 

5 + 5^2 + 5^3 + 3*5^6 + 4*5^7 + 5^9 + 5^10 + 3*5^11 + 3*5^12 + 5^13 + 4*5^14 + 5^15 + 2*5^16 + 5^17 + 2*5^18 + 4*5^19 + O(5^20) 

""" 

p = Integer(p) # this is assumed in code below 

if check_hypotheses: 

if not p.is_prime(): 

raise ValueError("p = (%s) must be prime"%p) 

if p == 2: 

raise ValueError("p must be odd") # todo 

if self.conductor() % (p**2) == 0: 

raise ArithmeticError("p must be a semi-stable prime") 

if self.conductor() % p == 0: 

Eq = self.tate_curve(p) 

reg = Eq.padic_regulator(prec=prec) 

return reg 

elif self.ap(p) % p == 0: 

lp = self.padic_lseries(p) 

reg = lp.Dp_valued_regulator(prec=prec) 

return reg 

else: 

if self.rank() == 0: 

return Qp(p,prec)(1) 

if height is None: 

height = self.padic_height(p, prec, check_hypotheses=False) 

d = self.padic_height_pairing_matrix(p=p, prec=prec, height=height, check_hypotheses=False) 

return d.determinant() 

def padic_height_pairing_matrix(self, p, prec=20, height=None, check_hypotheses=True): 

r""" 

Computes the cyclotomic `p`-adic height pairing matrix of 

this curve with respect to the basis self.gens() for the 

Mordell-Weil group for a given odd prime p of good ordinary 

reduction. 

INPUT: 

- ``p`` - prime = 5 

- ``prec`` - answer will be returned modulo 

`p^{\mathrm{prec}}` 

- ``height`` - precomputed height function. If not 

supplied, this function will call padic_height to compute it. 

- ``check_hypotheses`` - boolean, whether to check 

that this is a curve for which the p-adic height makes sense 

OUTPUT: The p-adic cyclotomic height pairing matrix of this curve 

to the given precision. 

.. TODO:: 

remove restriction that curve must be in minimal 

Weierstrass form. This is currently required for E.gens(). 

AUTHORS: 

- David Harvey, Liang Xiao, Robert Bradshaw, Jennifer 

Balakrishnan: original implementation at the 2006 MSRI graduate 

workshop on modular forms 

- David Harvey (2006-09-13): cleaned up and integrated into Sage, 

removed some redundant height computations 

EXAMPLES:: 

sage: E = EllipticCurve("37a") 

sage: E.padic_height_pairing_matrix(5, 10) 

[5 + 5^2 + 5^3 + 3*5^6 + 4*5^7 + 5^9 + O(5^10)] 

A rank two example:: 

sage: e =EllipticCurve('389a') 

sage: e._set_gens([e(-1, 1), e(1,0)]) # avoid platform dependent gens 

sage: e.padic_height_pairing_matrix(5,10) 

[ 3*5 + 2*5^2 + 5^4 + 5^5 + 5^7 + 4*5^9 + O(5^10) 5 + 4*5^2 + 5^3 + 2*5^4 + 3*5^5 + 4*5^6 + 5^7 + 5^8 + 2*5^9 + O(5^10)] 

[5 + 4*5^2 + 5^3 + 2*5^4 + 3*5^5 + 4*5^6 + 5^7 + 5^8 + 2*5^9 + O(5^10) 4*5 + 2*5^4 + 3*5^6 + 4*5^7 + 4*5^8 + O(5^10)] 

An anomalous rank 3 example:: 

sage: e = EllipticCurve("5077a") 

sage: e._set_gens([e(-1,3), e(2,0), e(4,6)]) 

sage: e.padic_height_pairing_matrix(5,4) 

[4 + 3*5 + 4*5^2 + 4*5^3 + O(5^4) 4 + 4*5^2 + 2*5^3 + O(5^4) 3*5 + 4*5^2 + 5^3 + O(5^4)] 

[ 4 + 4*5^2 + 2*5^3 + O(5^4) 3 + 4*5 + 3*5^2 + 5^3 + O(5^4) 2 + 4*5 + O(5^4)] 

[ 3*5 + 4*5^2 + 5^3 + O(5^4) 2 + 4*5 + O(5^4) 1 + 3*5 + 5^2 + 5^3 + O(5^4)] 

""" 

if check_hypotheses: 

p = __check_padic_hypotheses(self, p) 

K = Qp(p, prec=prec) 

rank = self.rank() 

M = matrix.matrix(K, rank, rank, 0) 

if rank == 0: 

return M 

basis = self.gens() 

if height is None: 

height = self.padic_height(p, prec, check_hypotheses=False) 

# Use <P, Q> =1/2*( h(P + Q) - h(P) - h(Q) ) 

for i in range(rank): 

M[i,i] = height(basis[i]) 

for i in range(rank): 

for j in range(i+1, rank): 

M[i, j] = ( height(basis[i] + basis[j]) - M[i,i] - M[j,j] ) / 2 

M[j, i] = M[i, j] 

return M 

def _multiply_point(E, R, P, m): 

r""" 

Computes coordinates of a multiple of P with entries in a ring. 

INPUT: 

- ``E`` - elliptic curve over Q with integer 

coefficients 

- ``P`` - a rational point on P that reduces to a 

non-singular point at all primes 

- ``R`` - a ring in which 2 is invertible (typically 

`\ZZ/L\ZZ` for some positive odd integer 

`L`). 

- ``m`` - an integer, = 1 

OUTPUT: A triple `(a', b', d')` such that if the point 

`mP` has coordinates `(a/d^2, b/d^3)`, then we have 

`a' \equiv a`, `b' \equiv \pm b`, 

`d' \equiv \pm d` all in `R` (i.e. modulo 

`L`). 

Note the ambiguity of signs for `b'` and `d'`. 

There's not much one can do about this, but at least one can say 

that the sign for `b'` will match the sign for 

`d'`. 

ALGORITHM: Proposition 9 of "Efficient Computation of p-adic 

Heights" (David Harvey, to appear in LMS JCM). 

Complexity is soft-`O(\log L \log m + \log^2 m)`. 

AUTHORS: 

- David Harvey (2008-01): replaced _DivPolyContext with 

_multiply_point 

EXAMPLES: 

37a has trivial Tamagawa numbers so all points have nonsingular 

reduction at all primes:: 

sage: E = EllipticCurve("37a") 

sage: P = E([0, -1]); P 

(0 : -1 : 1) 

sage: 19*P 

(-59997896/67387681 : 88075171080/553185473329 : 1) 

sage: R = Integers(625) 

sage: from sage.schemes.elliptic_curves.padics import _multiply_point 

sage: _multiply_point(E, R, P, 19) 

(229, 170, 541) 

sage: -59997896 % 625 

229 

sage: -88075171080 % 625 # note sign is flipped 

170 

sage: -67387681.sqrt() % 625 # sign is flipped here too 

541 

Trivial cases (:trac:`3632`):: 

sage: _multiply_point(E, R, P, 1) 

(0, 624, 1) 

sage: _multiply_point(E, R, 19*P, 1) 

(229, 455, 84) 

sage: (170 + 455) % 625 # note the sign did not flip here 

0 

sage: (541 + 84) % 625 

0 

Test over a range of `n` for a single curve with fairly 

random coefficients:: 

sage: R = Integers(625) 

sage: E = EllipticCurve([4, -11, 17, -8, -10]) 

sage: P = E.gens()[0] * LCM(E.tamagawa_numbers()) 

sage: from sage.schemes.elliptic_curves.padics import _multiply_point 

sage: Q = P 

sage: for n in range(1, 25): 

....: naive = R(Q[0].numerator()), \ 

....: R(Q[1].numerator()), \ 

....: R(Q[0].denominator().sqrt()) 

....: triple = _multiply_point(E, R, P, n) 

....: assert (triple[0] == naive[0]) and ( \ 

....: (triple[1] == naive[1] and triple[2] == naive[2]) or \ 

....: (triple[1] == -naive[1] and triple[2] == -naive[2])), \ 

....: "_multiply_point() gave an incorrect answer" 

....: Q = Q + P 

""" 

assert m >= 1 

alpha = R(P[0].numerator()) 

beta = R(P[1].numerator()) 

d = R(P[0].denominator().sqrt()) 

if m == 1: 

return alpha, beta, d 

a1 = R(E.a1()) * d 

a3 = R(E.a3()) * d**3 

b2 = R(E.b2()) * d**2 

b4 = R(E.b4()) * d**4 

b6 = R(E.b6()) * d**6 

b8 = R(E.b8()) * d**8 

B4 = 6*alpha**2 + b2*alpha + b4 

B6 = 4*alpha**3 + b2*alpha**2 + 2*b4*alpha + b6 

B6_sqr = B6*B6 

B8 = 3*alpha**4 + b2*alpha**3 + 3*b4*alpha**2 + 3*b6*alpha + b8 

T = 2*beta + a1*alpha + a3 

# make a list of disjoint intervals [a[i], b[i]) such that we need to 

# compute g(k) for all a[i] <= k <= b[i] for each i 

intervals = [] 

interval = (m - 2, m + 3) 

while interval[0] < interval[1]: 

intervals.append(interval) 

interval = max((interval[0] - 3) >> 1, 0), \ 

min((interval[1] + 5) >> 1, interval[0]) 

# now walk through list and compute g(k) 

g = {0 : R(0), 1 : R(1), 2 : R(-1), 3 : B8, 4 : B6**2 - B4*B8} 

last = [0, 1, 2, 3, 4] # last few k 

for i in reversed(intervals): 

k = i[0] 

while k < i[1]: 

if k > 4: 

j = k >> 1 

if k & 1: 

t1 = g[j] 

t2 = g[j+1] 

prod1 = g[j+2] * t1*t1*t1 

prod2 = g[j-1] * t2*t2*t2 

g[k] = prod1 - B6_sqr * prod2 if j & 1 else B6_sqr * prod1 - prod2 

else: 

t1 = g[j-1] 

t2 = g[j+1] 

g[k] = g[j] * (g[j-2] * t2*t2 - g[j+2] * t1*t1) 

k = k + 1 

if m & 1: 

psi_m = g[m] 

psi_m_m1 = g[m-1] * T 

psi_m_p1 = g[m+1] * T 

else: 

psi_m = g[m] * T 

psi_m_m1 = g[m-1] 

psi_m_p1 = g[m+1] 

theta = alpha * psi_m * psi_m - psi_m_m1 * psi_m_p1 

t1 = g[m-2] * g[m+1] * g[m+1] - g[m+2] * g[m-1] * g[m-1] 

if m & 1: 

t1 = t1 * T 

omega = (t1 + (a1 * theta + a3 * psi_m * psi_m) * psi_m) / -2 

return theta, omega, psi_m * d 

def padic_height(self, p, prec=20, sigma=None, check_hypotheses=True): 

r""" 

Compute the cyclotomic p-adic height. 

The equation of the curve must be minimal at `p`. 

INPUT: 

- ``p`` - prime = 5 for which the curve has 

semi-stable reduction 

- ``prec`` - integer = 1, desired precision of result 

- ``sigma`` - precomputed value of sigma. If not 

supplied, this function will call padic_sigma to compute it. 

- ``check_hypotheses`` - boolean, whether to check 

that this is a curve for which the p-adic height makes sense 

OUTPUT: A function that accepts two parameters: 

- a Q-rational point on the curve whose height should be computed 

- optional boolean flag 'check': if False, it skips some input 

checking, and returns the p-adic height of that point to the 

desired precision. 

- The normalization (sign and a factor 1/2 with respect to some other 

normalizations that appear in the literature) is chosen in such a way 

as to make the p-adic Birch Swinnerton-Dyer conjecture hold as stated 

in [Mazur-Tate-Teitelbaum]. 

AUTHORS: 

- Jennifer Balakrishnan: original code developed at the 2006 MSRI 

graduate workshop on modular forms 

- David Harvey (2006-09-13): integrated into Sage, optimised to 

speed up repeated evaluations of the returned height function, 

addressed some thorny precision questions 

- David Harvey (2006-09-30): rewrote to use division polynomials 

for computing denominator of `nP`. 

- David Harvey (2007-02): cleaned up according to algorithms in 

"Efficient Computation of p-adic Heights" 

- Chris Wuthrich (2007-05): added supersingular and multiplicative heights 

EXAMPLES:: 

sage: E = EllipticCurve("37a") 

sage: P = E.gens()[0] 

sage: h = E.padic_height(5, 10) 

sage: h(P) 

5 + 5^2 + 5^3 + 3*5^6 + 4*5^7 + 5^9 + O(5^10) 

An anomalous case:: 

sage: h = E.padic_height(53, 10) 

sage: h(P) 

26*53^-1 + 30 + 20*53 + 47*53^2 + 10*53^3 + 32*53^4 + 9*53^5 + 22*53^6 + 35*53^7 + 30*53^8 + 17*53^9 + O(53^10) 

Boundary case:: 

sage: E.padic_height(5, 3)(P) 

5 + 5^2 + O(5^3) 

A case that works the division polynomial code a little harder:: 

sage: E.padic_height(5, 10)(5*P) 

5^3 + 5^4 + 5^5 + 3*5^8 + 4*5^9 + O(5^10) 

Check that answers agree over a range of precisions:: 

sage: max_prec = 30 # make sure we get past p^2 # long time 

sage: full = E.padic_height(5, max_prec)(P) # long time 

sage: for prec in range(1, max_prec): # long time 

....: assert E.padic_height(5, prec)(P) == full # long time 

A supersingular prime for a curve:: 

sage: E = EllipticCurve('37a') 

sage: E.is_supersingular(3) 

True 

sage: h = E.padic_height(3, 5) 

sage: h(E.gens()[0]) 

(3 + 3^3 + O(3^6), 2*3^2 + 3^3 + 3^4 + 3^5 + 2*3^6 + O(3^7)) 

sage: E.padic_regulator(5) 

5 + 5^2 + 5^3 + 3*5^6 + 4*5^7 + 5^9 + 5^10 + 3*5^11 + 3*5^12 + 5^13 + 4*5^14 + 5^15 + 2*5^16 + 5^17 + 2*5^18 + 4*5^19 + O(5^20) 

sage: E.padic_regulator(3, 5) 

(3 + 2*3^2 + 3^3 + O(3^4), 3^2 + 2*3^3 + 3^4 + O(3^5)) 

A torsion point in both the good and supersingular cases:: 

sage: E = EllipticCurve('11a') 

sage: P = E.torsion_subgroup().gen(0).element(); P 

(5 : 5 : 1) 

sage: h = E.padic_height(19, 5) 

sage: h(P) 

0 

sage: h = E.padic_height(5, 5) 

sage: h(P) 

0 

The result is not dependent on the model for the curve:: 

sage: E = EllipticCurve([0,0,0,0,2^12*17]) 

sage: Em = E.minimal_model() 

sage: P = E.gens()[0] 

sage: Pm = Em.gens()[0] 

sage: h = E.padic_height(7) 

sage: hm = Em.padic_height(7) 

sage: h(P) == hm(Pm) 

True 

TESTS: 

Check that ticket :trac:`20798` is solved:: 

sage: E = EllipticCurve("91b") 

sage: h = E.padic_height(7,10) 

sage: P = E.gen(0) 

sage: h(P) 

2*7 + 7^2 + 5*7^3 + 6*7^4 + 2*7^5 + 3*7^6 + 7^7 + O(7^9) 

sage: h(P+P) 

7 + 5*7^2 + 6*7^3 + 5*7^4 + 4*7^5 + 6*7^6 + 5*7^7 + O(7^9) 

""" 

if check_hypotheses: 

if not p.is_prime(): 

raise ValueError("p = (%s) must be prime"%p) 

if p == 2: 

raise ValueError("p must be odd") # todo 

if self.conductor() % (p**2) == 0: 

raise ArithmeticError("p must be a semi-stable prime") 

prec = int(prec) 

if prec < 1: 

raise ValueError("prec (=%s) must be at least 1" % prec) 

if self.conductor() % p == 0: 

Eq = self.tate_curve(p) 

return Eq.padic_height(prec=prec) 

elif self.ap(p) % p == 0: 

lp = self.padic_lseries(p) 

return lp.Dp_valued_height(prec=prec) 

# else good ordinary case 

# For notation and definitions, see "Efficient Computation of 

# p-adic Heights", David Harvey (unpublished) 

n1 = self.change_ring(rings.GF(p)).cardinality() 

n2 = arith.LCM(self.tamagawa_numbers()) 

n = arith.LCM(n1, n2) 

m = int(n / n2) 

adjusted_prec = prec + 2 * arith.valuation(n, p) # this is M' 

R = rings.Integers(p ** adjusted_prec) 

if sigma is None: 

sigma = self.padic_sigma(p, adjusted_prec, check_hypotheses=False) 

# K is the field for the final result 

K = Qp(p, prec=adjusted_prec-1) 

E = self 

def height(P, check=True): 

if P.is_finite_order(): 

return K(0) 

if check: 

assert P.curve() == E, "the point P must lie on the curve " \ 

"from which the height function was created" 

Q = n2 * P 

alpha, beta, d = _multiply_point(E, R, Q, m) 

assert beta.lift() % p != 0, "beta should be a unit!" 

assert d.lift() % p == 0, "d should not be a unit!" 

t = -d * alpha / beta 

total = R(1) 

t_power = t 

for k in range(2, adjusted_prec + 1): 

total = total + t_power * sigma[k].lift() 

t_power = t_power * t 

total = (-alpha / beta) * total 

L = Qp(p, prec=adjusted_prec) 

total = L(total.lift(), adjusted_prec) # yuck... get rid of this lift! 

# changed sign to make it correct for p-adic bsd 

answer = -total.log() * 2 / n**2 

if check: 

assert answer.precision_absolute() >= prec, "we should have got an " \ 

"answer with precision at least prec, but we didn't." 

return K(answer) 

# (man... I love python's local function definitions...) 

return height 

def padic_height_via_multiply(self, p, prec=20, E2=None, check_hypotheses=True): 

r""" 

Computes the cyclotomic p-adic height. 

The equation of the curve must be minimal at `p`. 

INPUT: 

- ``p`` - prime = 5 for which the curve has good 

ordinary reduction 

- ``prec`` - integer = 2, desired precision of result 

- ``E2`` - precomputed value of E2. If not supplied, 

this function will call padic_E2 to compute it. The value supplied 

must be correct mod `p^(prec-2)` (or slightly higher in the 

anomalous case; see the code for details). 

- ``check_hypotheses`` - boolean, whether to check 

that this is a curve for which the p-adic height makes sense 

OUTPUT: A function that accepts two parameters: 

- a Q-rational point on the curve whose height should be computed 

- optional boolean flag 'check': if False, it skips some input 

checking, and returns the p-adic height of that point to the 

desired precision. 

- The normalization (sign and a factor 1/2 with respect to some other 

normalizations that appear in the literature) is chosen in such a way 

as to make the p-adic Birch Swinnerton-Dyer conjecture hold as stated 

in [Mazur-Tate-Teitelbaum]. 

AUTHORS: 

- David Harvey (2008-01): based on the padic_height() function, 

using the algorithm of"Computing p-adic heights via 

point multiplication" 

EXAMPLES:: 

sage: E = EllipticCurve("37a") 

sage: P = E.gens()[0] 

sage: h = E.padic_height_via_multiply(5, 10) 

sage: h(P) 

5 + 5^2 + 5^3 + 3*5^6 + 4*5^7 + 5^9 + O(5^10) 

An anomalous case:: 

sage: h = E.padic_height_via_multiply(53, 10) 

sage: h(P) 

26*53^-1 + 30 + 20*53 + 47*53^2 + 10*53^3 + 32*53^4 + 9*53^5 + 22*53^6 + 35*53^7 + 30*53^8 + 17*53^9 + O(53^10) 

Supply the value of E2 manually:: 

sage: E2 = E.padic_E2(5, 8) 

sage: E2 

2 + 4*5 + 2*5^3 + 5^4 + 3*5^5 + 2*5^6 + O(5^8) 

sage: h = E.padic_height_via_multiply(5, 10, E2=E2) 

sage: h(P) 

5 + 5^2 + 5^3 + 3*5^6 + 4*5^7 + 5^9 + O(5^10) 

Boundary case:: 

sage: E.padic_height_via_multiply(5, 3)(P) 

5 + 5^2 + O(5^3) 

Check that answers agree over a range of precisions:: 

sage: max_prec = 30 # make sure we get past p^2 # long time 

sage: full = E.padic_height(5, max_prec)(P) # long time 

sage: for prec in range(2, max_prec): # long time 

....: assert E.padic_height_via_multiply(5, prec)(P) == full # long time 

""" 

if check_hypotheses: 

if not p.is_prime(): 

raise ValueError("p = (%s) must be prime"%p) 

if p == 2: 

raise ValueError("p must be odd") # todo 

if self.conductor() % p == 0: 

raise ArithmeticError("must have good reduction at p") 

if self.ap(p) % p == 0: 

raise ArithmeticError("must be ordinary at p") 

prec = int(prec) 

if prec < 1: 

raise ValueError("prec (=%s) must be at least 1" % prec) 

# For notation and definitions, see ``Computing p-adic heights via point 

# multiplication'' (David Harvey, still in draft form) 

n1 = self.change_ring(rings.GF(p)).cardinality() 

n2 = arith.LCM(self.tamagawa_numbers()) 

n = arith.LCM(n1, n2) 

m = int(n / n2) 

lamb = int(math.floor(math.sqrt(prec))) 

adjusted_prec = prec + 2 * arith.valuation(n, p) # this is M' 

R = rings.Integers(p ** (adjusted_prec + 2*lamb)) 

sigma = self.padic_sigma_truncated(p, N=adjusted_prec, E2=E2, lamb=lamb)