r"""
Computation of Frobenius matrix on Monsky-Washnitzer cohomology
The most interesting functions to be exported here are
:func:`matrix_of_frobenius` and :func:`adjusted_prec`.
Currently this code is limited to the case `p \geq 5` (no
`GF(p^n)` for `n > 1`), and only handles the
elliptic curve case (not more general hyperelliptic curves).
REFERENCES:
- Kedlaya, K., "Counting points on hyperelliptic curves using
Monsky-Washnitzer cohomology", J. Ramanujan Math. Soc. 16 (2001) no
4, 323-338
- Edixhoven, B., "Point counting after Kedlaya", EIDMA-Stieltjes
graduate course, Lieden (lecture notes?).
AUTHORS:
- David Harvey and Robert Bradshaw: initial code developed at the 2006
MSRI graduate workshop, working with Jennifer Balakrishnan and Liang
Xiao
- David Harvey (2006-08): cleaned up, rewrote some chunks, lots more
documentation, added Newton iteration method, added more complete
'trace trick', integrated better into Sage.
- David Harvey (2007-02): added algorithm with sqrt(p) complexity
(removed in May 2007 due to better C++ implementation)
- Robert Bradshaw (2007-03): keep track of exact form in reduction
algorithms
- Robert Bradshaw (2007-04): generalization to hyperelliptic curves
"""
#*****************************************************************************
# Copyright (C) 2006 William Stein <wstein@gmail.com>
# 2006 Robert Bradshaw <robertwb@math.washington.edu>
# 2006 David Harvey <dmharvey@math.harvard.edu>
#
# Distributed under the terms of the GNU General Public License (GPL)
# http://www.gnu.org/licenses/
#*****************************************************************************
from sage.rings.all import Integers, Integer, PolynomialRing, is_Polynomial, PowerSeriesRing, Rationals, Rational, LaurentSeriesRing
from sage.modules.module import Module
from sage.structure.element import ModuleElement
from sage.matrix.all import matrix
from sage.modules.all import vector
from sage.rings.ring import CommutativeAlgebra
from sage.structure.element import CommutativeAlgebraElement
from sage.rings.arith import binomial, integer_ceil as ceil
from sage.misc.functional import log
from sage.misc.misc import newton_method_sizes
from ell_generic import is_EllipticCurve
from constructor import EllipticCurve
class SpecialCubicQuotientRing(CommutativeAlgebra):
r"""
Specialised class for representing the quotient ring
`R[x,T]/(T - x^3 - ax - b)`, where `R` is an
arbitrary commutative base ring (in which 2 and 3 are invertible),
`a` and `b` are elements of that ring.
Polynomials are represented internally in the form
`p_0 + p_1 x + p_2 x^2` where the `p_i` are
polynomials in `T`. Multiplication of polynomials always
reduces high powers of `x` (i.e. beyond `x^2`) to
powers of `T`.
Hopefully this ring is faster than a general quotient ring because
it uses the special structure of this ring to speed multiplication
(which is the dominant operation in the frobenius matrix
calculation). I haven't actually tested this theory though...
TODO: - Eventually we will want to run this in characteristic 3, so
we need to: (a) Allow `Q(x)` to contain an `x^2` term, and
(b) Remove the requirement that 3 be invertible. Currently this is
used in the Toom-Cook algorithm to speed multiplication.
EXAMPLES::
sage: B.<t> = PolynomialRing(Integers(125))
sage: R = monsky_washnitzer.SpecialCubicQuotientRing(t^3 - t + B(1/4))
sage: R
SpecialCubicQuotientRing over Ring of integers modulo 125 with polynomial T = x^3 + 124*x + 94
Get generators::
sage: x, T = R.gens()
sage: x
(0) + (1)*x + (0)*x^2
sage: T
(T) + (0)*x + (0)*x^2
Coercions::
sage: R(7)
(7) + (0)*x + (0)*x^2
Create elements directly from polynomials::
sage: A, z = R.poly_ring().objgen()
sage: A
Univariate Polynomial Ring in T over Ring of integers modulo 125
sage: R.create_element(z^2, z+1, 3)
(T^2) + (T + 1)*x + (3)*x^2
Some arithmetic::
sage: x^3
(T + 31) + (1)*x + (0)*x^2
sage: 3 * x**15 * T**2 + x - T
(3*T^7 + 90*T^6 + 110*T^5 + 20*T^4 + 58*T^3 + 26*T^2 + 124*T) + (15*T^6 + 110*T^5 + 35*T^4 + 63*T^2 + 1)*x + (30*T^5 + 40*T^4 + 8*T^3 + 38*T^2)*x^2
Retrieve coefficients (output is zero-padded)::
sage: x^10
(3*T^2 + 61*T + 8) + (T^3 + 93*T^2 + 12*T + 40)*x + (3*T^2 + 61*T + 9)*x^2
sage: (x^10).coeffs()
[[8, 61, 3, 0], [40, 12, 93, 1], [9, 61, 3, 0]]
TODO: write an example checking multiplication of these polynomials
against Sage's ordinary quotient ring arithmetic. I can't seem to
get the quotient ring stuff happening right now...
"""
def __init__(self, Q, laurent_series = False):
"""
Constructor.
INPUT:
- ``Q`` - a polynomial of the form
`Q(x) = x^3 + ax + b`, where `a`, `b` belong to a ring in which
2, 3 are invertible.
- ``laurent_series`` - whether or not to allow
negative powers of `T` (default=False)
EXAMPLES::
sage: B.<t> = PolynomialRing(Integers(125))
sage: R = monsky_washnitzer.SpecialCubicQuotientRing(t^3 - t + B(1/4))
sage: R
SpecialCubicQuotientRing over Ring of integers modulo 125 with polynomial T = x^3 + 124*x + 94
::
sage: R = monsky_washnitzer.SpecialCubicQuotientRing(t^3 + 2*t^2 - t + B(1/4))
Traceback (most recent call last):
...
ValueError: Q (=t^3 + 2*t^2 + 124*t + 94) must be of the form x^3 + ax + b
::
sage: B.<t> = PolynomialRing(Integers(10))
sage: R = monsky_washnitzer.SpecialCubicQuotientRing(t^3 - t + 1)
Traceback (most recent call last):
...
ArithmeticError: 2 and 3 must be invertible in the coefficient ring (=Ring of integers modulo 10) of Q
"""
if not is_Polynomial(Q):
raise TypeError, "Q (=%s) must be a polynomial" % Q
if Q.degree() != 3 or not Q[2].is_zero():
raise ValueError, "Q (=%s) must be of the form x^3 + ax + b" % Q
base_ring = Q.parent().base_ring()
if not base_ring(6).is_unit():
raise ArithmeticError, \
"2 and 3 must be invertible in the coefficient ring (=%s) of Q" % \
base_ring
# CommutativeAlgebra.__init__ tries to establish a coercion
# from the base ring, by trac ticket #9138. The corresponding
# hom set is cached. In order to use self as cache key, its
# string representation is used. In otder to get the string
# representation, we need to know the attributes _a and
# _b. Hence, in #9138, we have to move CommutativeAlgebra.__init__
# further down:
self._a = Q[1]
self._b = Q[0]
if laurent_series:
self._poly_ring = LaurentSeriesRing(base_ring, 'T') # R[T]
else:
self._poly_ring = PolynomialRing(base_ring, 'T') # R[T]
self._poly_generator = self._poly_ring.gen(0) # the generator T
CommutativeAlgebra.__init__(self, base_ring)
# Precompute a matrix that is used in the Toom-Cook multiplication.
# This is where we need 2 and 3 invertible.
# (a good description of Toom-Cook is online at:
# http://www.gnu.org/software/gmp/manual/html_node/Toom-Cook-3-Way-Multiplication.html )
self._speedup_matrix = \
(matrix(Integers(), 3, 3, [2, 4, 8,
1, 1, 1,
8, 4, 2])**(-1)
).change_ring(base_ring).list()
def __repr__(self):
"""
EXAMPLES::
sage: B.<t> = PolynomialRing(Integers(125))
sage: R = monsky_washnitzer.SpecialCubicQuotientRing(t^3 - t + B(1/4))
sage: print R
SpecialCubicQuotientRing over Ring of integers modulo 125 with polynomial T = x^3 + 124*x + 94
"""
return "SpecialCubicQuotientRing over %s with polynomial T = %s" % \
(self.base_ring(), PolynomialRing(self.base_ring(), 'x')(
[self._b, self._a, 0, 1]))
def poly_ring(self):
"""
Return the underlying polynomial ring in `T`.
EXAMPLES::
sage: B.<t> = PolynomialRing(Integers(125))
sage: R = monsky_washnitzer.SpecialCubicQuotientRing(t^3 - t + B(1/4))
sage: R.poly_ring()
Univariate Polynomial Ring in T over Ring of integers modulo 125
"""
return self._poly_ring
def gens(self):
"""
Return a list [x, T] where x and T are the generators of the ring
(as element *of this ring*).
.. note::
I have no idea if this is compatible with the usual Sage
'gens' interface.
EXAMPLES::
sage: B.<t> = PolynomialRing(Integers(125))
sage: R = monsky_washnitzer.SpecialCubicQuotientRing(t^3 - t + B(1/4))
sage: x, T = R.gens()
sage: x
(0) + (1)*x + (0)*x^2
sage: T
(T) + (0)*x + (0)*x^2
"""
return [SpecialCubicQuotientRingElement(self,
self._poly_ring(0), self._poly_ring(1), self._poly_ring(0),
check=False),
SpecialCubicQuotientRingElement(self,
self._poly_generator, self._poly_ring(0), self._poly_ring(0),
check=False)]
def create_element(self, p0, p1, p2, check=True):
"""
Creates the element `p_0 + p_1*x + p_2*x^2`, where the `p_i`
are polynomials in `T`.
INPUT:
- ``p0, p1, p2`` - coefficients; must be coercible
into poly_ring()
- ``check`` - bool (default True): whether to carry
out coercion
EXAMPLES::
sage: B.<t> = PolynomialRing(Integers(125))
sage: R = monsky_washnitzer.SpecialCubicQuotientRing(t^3 - t + B(1/4))
sage: A, z = R.poly_ring().objgen()
sage: R.create_element(z^2, z+1, 3)
(T^2) + (T + 1)*x + (3)*x^2
"""
return SpecialCubicQuotientRingElement(self, p0, p1, p2, check)
def __call__(self, value):
"""
EXAMPLES::
sage: B.<t> = PolynomialRing(Integers(125))
sage: R = monsky_washnitzer.SpecialCubicQuotientRing(t^3 - t + B(1/4))
sage: R(3)
(3) + (0)*x + (0)*x^2
"""
return self._coerce_(value)
def _coerce_impl(self, value):
"""
EXAMPLES::
sage: B.<t> = PolynomialRing(Integers(125))
sage: R = monsky_washnitzer.SpecialCubicQuotientRing(t^3 - t + B(1/4))
sage: R._coerce_impl(3)
(3) + (0)*x + (0)*x^2
"""
# coerce to underlying polynomial ring (possibly via base ring):
value = self._poly_ring._coerce_(value)
return SpecialCubicQuotientRingElement(self,
value, self._poly_ring(0), self._poly_ring(0), check=False)
class SpecialCubicQuotientRingElement(CommutativeAlgebraElement):
"""
An element of a SpecialCubicQuotientRing.
"""
def __init__(self, parent, p0, p1, p2, check=True):
"""
Constructs the element `p_0 + p_1*x + p_2*x^2`, where
the `p_i` are polynomials in `T`.
INPUT:
- ``parent`` - a SpecialCubicQuotientRing
- ``p0, p1, p2`` - coefficients; must be coercible
into parent.poly_ring()
- ``check`` - bool (default True): whether to carry
out coercion
EXAMPLES::
sage: B.<t> = PolynomialRing(Integers(125))
sage: R = monsky_washnitzer.SpecialCubicQuotientRing(t^3 - t + B(1/4))
sage: from sage.schemes.elliptic_curves.monsky_washnitzer import SpecialCubicQuotientRingElement
sage: SpecialCubicQuotientRingElement(R, 2, 3, 4)
(2) + (3)*x + (4)*x^2
"""
if not isinstance(parent, SpecialCubicQuotientRing):
raise TypeError, \
"parent (=%s) must be a SpecialCubicQuotientRing" % parent
CommutativeAlgebraElement.__init__(self, parent)
if check:
poly_ring = parent.poly_ring()
p0 = poly_ring(p0)
p1 = poly_ring(p1)
p2 = poly_ring(p2)
self._triple = (p0, p1, p2)
def coeffs(self):
"""
Returns list of three lists of coefficients, corresponding to the
`x^0`, `x^1`, `x^2` coefficients. The lists
are zero padded to the same length. The list entries belong to the
base ring.
EXAMPLES::
sage: B.<t> = PolynomialRing(Integers(125))
sage: R = monsky_washnitzer.SpecialCubicQuotientRing(t^3 - t + B(1/4))
sage: p = R.create_element(t, t^2 - 2, 3)
sage: p.coeffs()
[[0, 1, 0], [123, 0, 1], [3, 0, 0]]
"""
coeffs = [column.coeffs() for column in self._triple]
degree = max([len(x) for x in coeffs])
base_ring = self.parent().base_ring()
for column in coeffs:
column.extend([base_ring(0)] * (degree - len(column)))
return coeffs
def __nonzero__(self):
"""
EXAMPLES::
sage: B.<t> = PolynomialRing(Integers(125))
sage: R = monsky_washnitzer.SpecialCubicQuotientRing(t^3 - t + B(1/4))
sage: x, T = R.gens()
sage: not x
False
sage: not T
False
sage: not R.create_element(0, 0, 0)
True
"""
return not not self._triple[0] or not not self._triple[1] or not not self._triple[2]
def __cmp__(self, other):
"""
EXAMPLES::
sage: B.<t> = PolynomialRing(Integers(125))
sage: x, t = monsky_washnitzer.SpecialCubicQuotientRing(t^3 - t + B(1/4)).gens()
sage: x == t
False
sage: x == x
True
sage: x == x + x - x
True
"""
return cmp(self._triple, other._triple)
def _repr_(self):
"""
EXAMPLES::
sage: B.<t> = PolynomialRing(Integers(125))
sage: R = monsky_washnitzer.SpecialCubicQuotientRing(t^3 - t + B(1/4))
sage: x, T = R.gens()
sage: x + T*x - 2*T^2
(123*T^2) + (T + 1)*x + (0)*x^2
"""
return "(%s) + (%s)*x + (%s)*x^2" % self._triple
def _latex_(self):
"""
EXAMPLES::
sage: B.<t> = PolynomialRing(Integers(125))
sage: R = monsky_washnitzer.SpecialCubicQuotientRing(t^3 - t + B(1/4))
sage: x, T = R.gens()
sage: f = x + T*x - 2*T^2
sage: latex(f)
(123 T^{2}) + (T + 1)x + (0)x^2
"""
return "(%s) + (%s)x + (%s)x^2" % \
tuple([column._latex_() for column in self._triple])
def _add_(self, other):
"""
EXAMPLES::
sage: B.<t> = PolynomialRing(Integers(125))
sage: R = monsky_washnitzer.SpecialCubicQuotientRing(t^3 - t + B(1/4))
sage: f = R.create_element(2, t, t^2 - 3)
sage: g = R.create_element(3 + t, -t, t)
sage: f + g
(T + 5) + (0)*x + (T^2 + T + 122)*x^2
"""
return SpecialCubicQuotientRingElement(self.parent(),
self._triple[0] + other._triple[0],
self._triple[1] + other._triple[1],
self._triple[2] + other._triple[2],
check=False)
def _sub_(self, other):
"""
EXAMPLES::
sage: B.<t> = PolynomialRing(Integers(125))
sage: R = monsky_washnitzer.SpecialCubicQuotientRing(t^3 - t + B(1/4))
sage: f = R.create_element(2, t, t^2 - 3)
sage: g = R.create_element(3 + t, -t, t)
sage: f - g
(124*T + 124) + (2*T)*x + (T^2 + 124*T + 122)*x^2
"""
return SpecialCubicQuotientRingElement(self.parent(),
self._triple[0] - other._triple[0],
self._triple[1] - other._triple[1],
self._triple[2] - other._triple[2],
check=False)
def shift(self, n):
"""
Returns this element multiplied by `T^n`.
EXAMPLES::
sage: B.<t> = PolynomialRing(Integers(125))
sage: R = monsky_washnitzer.SpecialCubicQuotientRing(t^3 - t + B(1/4))
sage: f = R.create_element(2, t, t^2 - 3)
sage: f
(2) + (T)*x + (T^2 + 122)*x^2
sage: f.shift(1)
(2*T) + (T^2)*x + (T^3 + 122*T)*x^2
sage: f.shift(2)
(2*T^2) + (T^3)*x + (T^4 + 122*T^2)*x^2
"""
return SpecialCubicQuotientRingElement(self.parent(),
self._triple[0].shift(n),
self._triple[1].shift(n),
self._triple[2].shift(n),
check=False)
def scalar_multiply(self, scalar):
"""
Multiplies this element by a scalar, i.e. just multiply each
coefficient of `x^j` by the scalar.
INPUT:
- ``scalar`` - either an element of base_ring, or an
element of poly_ring.
EXAMPLES::
sage: B.<t> = PolynomialRing(Integers(125))
sage: R = monsky_washnitzer.SpecialCubicQuotientRing(t^3 - t + B(1/4))
sage: x, T = R.gens()
sage: f = R.create_element(2, t, t^2 - 3)
sage: f
(2) + (T)*x + (T^2 + 122)*x^2
sage: f.scalar_multiply(2)
(4) + (2*T)*x + (2*T^2 + 119)*x^2
sage: f.scalar_multiply(t)
(2*T) + (T^2)*x + (T^3 + 122*T)*x^2
"""
scalar = self.parent()._poly_ring(scalar)
return SpecialCubicQuotientRingElement(self.parent(),
scalar * self._triple[0],
scalar * self._triple[1],
scalar * self._triple[2],
check=False)
def square(self):
"""
EXAMPLES::
sage: B.<t> = PolynomialRing(Integers(125))
sage: R = monsky_washnitzer.SpecialCubicQuotientRing(t^3 - t + B(1/4))
sage: x, T = R.gens()
::
sage: f = R.create_element(1 + 2*t + 3*t^2, 4 + 7*t + 9*t^2, 3 + 5*t + 11*t^2)
sage: f.square()
(73*T^5 + 16*T^4 + 38*T^3 + 39*T^2 + 70*T + 120) + (121*T^5 + 113*T^4 + 73*T^3 + 8*T^2 + 51*T + 61)*x + (18*T^4 + 60*T^3 + 22*T^2 + 108*T + 31)*x^2
"""
return self * self
def _mul_(self, other):
"""
EXAMPLES::
sage: B.<t> = PolynomialRing(Integers(125))
sage: R = monsky_washnitzer.SpecialCubicQuotientRing(t^3 - t + B(1/4))
sage: x, T = R.gens()
::
sage: f = R.create_element(1 + 2*t + 3*t^2, 4 + 7*t + 9*t^2, 3 + 5*t + 11*t^2)
sage: g = R.create_element(4 + 3*t + 7*t^2, 2 + 3*t + t^2, 8 + 4*t + 6*t^2)
sage: f * g
(65*T^5 + 27*T^4 + 33*T^3 + 75*T^2 + 120*T + 57) + (66*T^5 + T^4 + 123*T^3 + 95*T^2 + 24*T + 50)*x + (45*T^4 + 75*T^3 + 37*T^2 + 2*T + 52)*x^2
"""
if not isinstance(other, SpecialCubicQuotientRingElement):
return self.scalar_multiply(other)
# Here we do Toom-Cook three-way multiplication, which reduces the
# naive 9 polynomial multiplications to only 5 polynomial multiplications.
a0, a1, a2 = self._triple
b0, b1, b2 = other._triple
M = self.parent()._speedup_matrix
if self is other:
# faster method if we're squaring
p0 = a0 * a0
temp = a0 + 2*a1 + 4*a2
p1 = temp * temp
temp = a0 + a1 + a2
p2 = temp * temp
temp = 4*a0 + 2*a1 + a2
p3 = temp * temp
p4 = a2 * a2
else:
p0 = a0 * b0
p1 = (a0 + 2*a1 + 4*a2) * (b0 + 2*b1 + 4*b2)
p2 = (a0 + a1 + a2) * (b0 + b1 + b2)
p3 = (4*a0 + 2*a1 + a2) * (4*b0 + 2*b1 + b2)
p4 = a2 * b2
q1 = p1 - p0 - 16*p4
q2 = p2 - p0 - p4
q3 = p3 - 16*p0 - p4
c0 = p0
c1 = M[0]*q1 + M[1]*q2 + M[2]*q3
c2 = M[3]*q1 + M[4]*q2 + M[5]*q3
c3 = M[6]*q1 + M[7]*q2 + M[8]*q3
c4 = p4
# Now the product is c0 + c1 x + c2 x^2 + c3 x^3 + c4 x^4.
# We need to reduce mod y = x^3 + ax + b and return result.
parent = self.parent()
T = parent._poly_generator
b = parent._b
a = parent._a
# todo: These lines are necessary to get binop stuff working
# for certain base rings, e.g. when we compute b*c3 in the
# final line. They shouldn't be necessary. Need to fix this
# somewhere else in Sage.
a = parent._poly_ring(a)
b = parent._poly_ring(b)
return SpecialCubicQuotientRingElement(parent,
-b*c3 + c0 + c3*T,
-b*c4 - a*c3 + c1 + c4*T,
-a*c4 + c2,
check=False)
def transpose_list(input):
"""
INPUT:
- ``input`` - a list of lists, each list of the same
length
OUTPUT:
- ``output`` - a list of lists such that output[i][j]
= input[j][i]
EXAMPLES::
sage: from sage.schemes.elliptic_curves.monsky_washnitzer import transpose_list
sage: L = [[1, 2], [3, 4], [5, 6]]
sage: transpose_list(L)
[[1, 3, 5], [2, 4, 6]]
"""
h = len(input)
w = len(input[0])
output = []
for i in range(w):
row = []
for j in range(h):
row.append(input[j][i])
output.append(row)
return output
def helper_matrix(Q):
"""
Computes the (constant) matrix used to calculate the linear
combinations of the `d(x^i y^j)` needed to eliminate the
negative powers of `y` in the cohomology (i.e. in
reduce_negative()).
INPUT:
- ``Q`` - cubic polynomial
"""
a = Q[1]
b = Q[0]
# Discriminant (should be invertible for a curve of good reduction)
D = 4*a**3 + 27*b**2
Dinv = D**(-1) # NB do not use 1/D
# This is the inverse of the matrix
# [ a, -3b, 0 ]
# [ 0, -2a, -3b ]
# [ 3, 0, -2a ]
M = Dinv * matrix( [ [4*a**2 , -6*b*a , 9*b**2 ],
[-9*b , -2*a**2 , 3*b*a ],
[ 6*a , -9*b , -2*a**2 ] ])
return M
def lift(x):
r"""
Tries to call x.lift(), presumably from the `p`-adics to ZZ.
If this fails, it assumes the input is a power series, and tries to
lift it to a power series over QQ.
This function is just a very kludgy solution to the problem of
trying to make the reduction code (below) work over both Zp and
Zp[[t]].
"""
try:
return x.lift()
except AttributeError:
return PowerSeriesRing(Rationals(), "t")(x.list(), x.prec())
def reduce_negative(Q, p, coeffs, offset, exact_form=None):
"""
Applies cohomology relations to incorporate negative powers of
`y` into the `y^0` term.
INPUT:
- ``p`` - prime
- ``Q`` - cubic polynomial
- ``coeffs`` - list of length 3 lists. The
`i^{th}` list [a, b, c] represents
`y^{2(i - offset)} (a + bx + cx^2) dx/y`.
- ``offset`` - nonnegative integer
OUTPUT: The reduction is performed in-place. The output is placed
in coeffs[offset]. Note that coeffs[i] will be meaningless for i
offset after this function is finished.
EXAMPLE::
sage: R.<x> = Integers(5^3)['x']
sage: Q = x^3 - x + R(1/4)
sage: coeffs = [[10, 15, 20], [1, 2, 3], [4, 5, 6], [7, 8, 9]]
sage: coeffs = [[R.base_ring()(a) for a in row] for row in coeffs]
sage: monsky_washnitzer.reduce_negative(Q, 5, coeffs, 3)
sage: coeffs[3]
[28, 52, 9]
::
sage: R.<x> = Integers(7^3)['x']
sage: Q = x^3 - x + R(1/4)
sage: coeffs = [[7, 14, 21], [1, 2, 3], [4, 5, 6], [7, 8, 9]]
sage: coeffs = [[R.base_ring()(a) for a in row] for row in coeffs]
sage: monsky_washnitzer.reduce_negative(Q, 7, coeffs, 3)
sage: coeffs[3]
[245, 332, 9]
"""
m = helper_matrix(Q).list()
base_ring = Q.base_ring()
next_a = coeffs[0]
if exact_form is not None:
x = exact_form.parent().gen(0)
y = exact_form.parent()(exact_form.parent().base_ring().gen(0))
try:
three_j_plus_5 = 5 - base_ring(6*offset)
three_j_plus_7 = 7 - base_ring(6*offset)
six = base_ring(6)
for i in range(0, offset):
j = 2*(i-offset)
a = next_a
next_a = coeffs[i+1]
# todo: the following divisions will sometimes involve
# a division by (a power of) p. In all cases, we know (from
# Kedlaya's estimates) that the answer should be p-integral.
# However, since we're working over $Z/p^k Z$, we're not allowed
# to "divide by p". So currently we lift to Q, divide, and coerce
# back. Eventually, when pAdicInteger is implemented, and plays
# nicely with pAdicField, we should reimplement this stuff
# using pAdicInteger.
if (p.divides(j+1)):
# need to lift here to perform the division
a[0] = base_ring(lift(a[0]) / (j+1))
a[1] = base_ring(lift(a[1]) / (j+1))
a[2] = base_ring(lift(a[2]) / (j+1))
else:
j_plus_1_inv = ~base_ring(j+1)
a[0] = a[0] * j_plus_1_inv
a[1] = a[1] * j_plus_1_inv
a[2] = a[2] * j_plus_1_inv
c1 = m[3]*a[0] + m[4]*a[1] + m[5]*a[2]
c2 = m[6]*a[0] + m[7]*a[1] + m[8]*a[2]
next_a[0] = next_a[0] - three_j_plus_5 * c1
next_a[1] = next_a[1] - three_j_plus_7 * c2
three_j_plus_7 = three_j_plus_7 + six
three_j_plus_5 = three_j_plus_5 + six
if exact_form is not None:
c0 = m[0]*a[0] + m[1]*a[1] + m[2]*a[2]
exact_form += (c0 + c1*x + c2 * x**2) * y**(j+1)
except NotImplementedError:
raise NotImplementedError, \
"It looks like you've found a non-integral matrix of Frobenius! " \
"(Q=%s, p=%s)\nTime to write a paper." % (Q, p)
coeffs[int(offset)] = next_a
return exact_form
def reduce_positive(Q, p, coeffs, offset, exact_form=None):
"""
Applies cohomology relations to incorporate positive powers of
`y` into the `y^0` term.
INPUT:
- ``Q`` - cubic polynomial
- ``coeffs`` - list of length 3 lists. The
`i^{th}` list [a, b, c] represents
`y^{2(i - offset)} (a + bx + cx^2) dx/y`.
- ``offset`` - nonnegative integer
OUTPUT: The reduction is performed in-place. The output is placed
in coeffs[offset]. Note that coeffs[i] will be meaningless for i
offset after this function is finished.
EXAMPLE::
sage: R.<x> = Integers(5^3)['x']
sage: Q = x^3 - x + R(1/4)
::
sage: coeffs = [[1, 2, 3], [10, 15, 20]]
sage: coeffs = [[R.base_ring()(a) for a in row] for row in coeffs]
sage: monsky_washnitzer.reduce_positive(Q, 5, coeffs, 0)
sage: coeffs[0]
[16, 102, 88]
::
sage: coeffs = [[9, 8, 7], [10, 15, 20]]
sage: coeffs = [[R.base_ring()(a) for a in row] for row in coeffs]
sage: monsky_washnitzer.reduce_positive(Q, 5, coeffs, 0)
sage: coeffs[0]
[24, 108, 92]
"""
base_ring = Q.base_ring()
next_a = coeffs[len(coeffs) - 1]
Qa = Q[1]
Qb = Q[0]
A = 2*Qa
B = 3*Qb
offset = Integer(offset)
if exact_form is not None:
x = exact_form.parent().gen(0)
y = exact_form.parent().base_ring().gen(0)
# y = exact_form.parent()(exact_form.parent().base_ring().gen(0))
for i in range(len(coeffs)-1, offset, -1):
j = 2*(i-offset) - 2
a = next_a
next_a = coeffs[i-1]
a[0] = a[0] - Qa*a[2]/3 # subtract d(y^j + 3)
if exact_form is not None:
exact_form += Q.base_ring()(a[2].lift() / (3*j+9)) * y**(j+3)
# todo: see comments about pAdicInteger in reduceNegative()
# subtract off c1 of d(x y^j + 1), and
if p.divides(3*j + 5):
c1 = base_ring(lift(a[0]) / (3*j + 5))
else:
c1 = a[0] / (3*j + 5)
# subtract off c2 of d(x^2 y^j + 1)
if p.divides(3*j + 7):
c2 = base_ring(lift(a[1]) / (3*j + 7))
else:
c2 = a[1] / (3*j + 7)
next_a[0] = next_a[0] + B*c1*(j+1)
next_a[1] = next_a[1] + A*c1*(j+1) + B*c2*(j+1)
next_a[2] = next_a[2] + A*c2*(j+1)
if exact_form is not None:
exact_form += (c1*x + c2 * x**2) * y**(j+1)
coeffs[int(offset)] = next_a
return exact_form
def reduce_zero(Q, coeffs, offset, exact_form=None):
"""
Applies cohomology relation to incorporate `x^2 y^0` term
into `x^0 y^0` and `x^1 y^0` terms.
INPUT:
- ``Q`` - cubic polynomial
- ``coeffs`` - list of length 3 lists. The
`i^{th}` list [a, b, c] represents
`y^{2(i - offset)} (a + bx + cx^2) dx/y`.
- ``offset`` - nonnegative integer
OUTPUT: The reduction is performed in-place. The output is placed
in coeffs[offset]. This method completely ignores coeffs[i] for i
!= offset.
EXAMPLE::
sage: R.<x> = Integers(5^3)['x']
sage: Q = x^3 - x + R(1/4)
sage: coeffs = [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
sage: coeffs = [[R.base_ring()(a) for a in row] for row in coeffs]
sage: monsky_washnitzer.reduce_zero(Q, coeffs, 1)
sage: coeffs[1]
[6, 5, 0]
"""
a = coeffs[int(offset)]
if a[2] == 0:
return exact_form
Qa = Q[1]
a[0] = a[0] - a[2]*Qa/3 # $3x^2 dx/y = -a dx/y$
coeffs[int(offset)] = a
if exact_form is not None:
x = exact_form.parent().gen(0)
y = exact_form.parent()(exact_form.parent().base_ring().gen(0))
exact_form += Q.base_ring()(a[2] / 3) * y
a[2] = 0
coeffs[int(offset)] = a
return exact_form
def reduce_all(Q, p, coeffs, offset, compute_exact_form=False):
"""
Applies cohomology relations to reduce all terms to a linear
combination of `dx/y` and `x dx/y`.
INPUT:
- ``Q`` - cubic polynomial
- ``coeffs`` - list of length 3 lists. The
`i^{th}` list [a, b, c] represents
`y^{2(i - offset)} (a + bx + cx^2) dx/y`.
- ``offset`` - nonnegative integer
OUTPUT:
- ``A, B`` - pair such that the input differential is
cohomologous to (A + Bx) dx/y.
.. note::
The algorithm operates in-place, so the data in coeffs is
destroyed.
EXAMPLE::
sage: R.<x> = Integers(5^3)['x']
sage: Q = x^3 - x + R(1/4)
sage: coeffs = [[1, 2, 3], [4, 5, 6], [7, 8, 9]]
sage: coeffs = [[R.base_ring()(a) for a in row] for row in coeffs]
sage: monsky_washnitzer.reduce_all(Q, 5, coeffs, 1)
(21, 106)
"""
R = Q.base_ring()
if compute_exact_form:
# exact_form = SpecialCubicQuotientRing(Q, laurent_series=True)(0)
exact_form = PolynomialRing(LaurentSeriesRing(Q.base_ring(), 'y'), 'x')(0)
# t = (Q.base_ring().order().factor())[0]
# from sage.rings.padics.qp import pAdicField
# exact_form = PolynomialRing(LaurentSeriesRing(pAdicField(p, t[1]), 'y'), 'x')(0)
else:
exact_form = None
while len(coeffs) <= offset:
coeffs.append([R(0), R(0), R(0)])
exact_form = reduce_negative(Q, p, coeffs, offset, exact_form)
exact_form = reduce_positive(Q, p, coeffs, offset, exact_form)
exact_form = reduce_zero(Q, coeffs, offset, exact_form)
if exact_form is None:
return coeffs[int(offset)][0], coeffs[int(offset)][1]
else:
return (coeffs[int(offset)][0], coeffs[int(offset)][1]), exact_form
def frobenius_expansion_by_newton(Q, p, M):
r"""
Computes the action of Frobenius on `dx/y` and on
`x dx/y`, using Newton's method (as suggested in Kedlaya's
paper).
(This function does *not* yet use the cohomology relations - that
happens afterwards in the "reduction" step.)
More specifically, it finds `F_0` and `F_1` in
the quotient ring `R[x, T]/(T - Q(x))`, such that
.. math::
F( dx/y) = T^{-r} F0 dx/y, \text{\ and\ } F(x dx/y) = T^{-r} F1 dx/y
where
.. math::
r = ( (2M-3)p - 1 )/2.
(Here `T` is `y^2 = z^{-2}`, and `R` is the
coefficient ring of `Q`.)
`F_0` and `F_1` are computed in the
SpecialCubicQuotientRing associated to `Q`, so all powers
of `x^j` for `j \geq 3` are reduced to powers of
`T`.
INPUT:
- ``Q`` - cubic polynomial of the form
`Q(x) = x^3 + ax + b`, whose coefficient ring is a
`Z/(p^M)Z`-algebra
- ``p`` - residue characteristic of the p-adic field
- ``M`` - p-adic precision of the coefficient ring
(this will be used to determine the number of Newton iterations)
OUTPUT:
- ``F0, F1`` - elements of
SpecialCubicQuotientRing(Q), as described above
- ``r`` - non-negative integer, as described above
"""
S = SpecialCubicQuotientRing(Q)
x, _ = S.gens() # T = y^2
base_ring = S.base_ring()
# When we compute Frob(1/y) we actually only need precision M-1, since
# we're going to multiply by p at the end anyway.
M = float(M - 1)
# Kedlaya sets s = Q(x^p)/T^p = 1 + p T^{-p} E, where
# E = (Q(x^p) - Q(x)^p) / p (has integral coefficients).
# Then he computes s^{-1/2} in S, using Newton's method to find
# successive approximations. We follow this plan, but we normalise our
# approximations so that we only ever need positive powers of T.
# Start by setting r = Q(x^p)/2 = 1/2 T^p s.
# (The 1/2 is for convenience later on.)
x_to_p_less_one = x**(p-1)
x_to_p = x_to_p_less_one * x
x_to_p_cubed = x_to_p.square() * x_to_p
r = (base_ring(1) / base_ring(2)) * (x_to_p_cubed + Q[1]*x_to_p + S(Q[0]))
# todo: this next loop would be clearer if it used the newton_method_sizes()
# function
# We will start with a hard-coded initial approximation, which we provide
# up to precision 3. First work out what precision is best to start with.
if M <= 3:
initial_precision = M
elif ceil(log(M/2, 2)) == ceil(log(M/3, 2)):
# In this case there's no advantage to starting with precision three,
# because we'll overshoot at the end. E.g. suppose the final precision
# is 8. If we start with precision 2, we need two iterations to get us
# to 8. If we start at precision 3, we will still need two iterations,
# but we do more work along the way. So may as well start with only 2.
initial_precision = 2
else:
initial_precision = 3
# Now compute the first approximation. In the main loop below, X is the
# normalised approximation, and k is the precision. More specifically,
# X = T^{p(k-1)} x_i, where x_i is an approximation to s^{-1/2}, and the
# approximation is correct mod p^k.
if initial_precision == 1:
k = 1
X = S(1)
elif initial_precision == 2:
# approximation is 3/2 - 1/2 s
k = 2
X = S(base_ring(3) / base_ring(2)).shift(p) - r
elif initial_precision == 3:
# approximation is (15 - 10 s + 3 s^2) / 8
k = 3
X = (base_ring(1) / base_ring(8)) * (
S(15).shift(2*p) - \
(base_ring(20) * r).shift(p) + \
(base_ring(12) * r.square()) \
)
# The key to the following calculation is that the T^{-m} coefficient
# of every x_i is divisible by p^(ceil(m/p)) (for m >= 0). Therefore if
# we are only expecting an answer correct mod p^k, we can truncate
# beyond the T^{-(k-1)p} term without any problems.
# todo: what would be really nice is to be able to work in a lower
# precision *coefficient ring* when we start the iteration, and move up to
# higher precision rings as the iteration proceeds. This would be feasible
# over Integers(p**n), but quite complicated (maybe impossible) over a more
# general base ring. This might give a decent constant factor speedup;
# or it might not, depending on how much the last iteration dominates the
# whole runtime. My guess is that it isn't worth the effort.
three_halves = base_ring(3) / base_ring(2)
# Newton iteration loop
while k < M:
# target_k = k' = precision we want our answer to be after this iteration
target_k = 2*k
# This prevents us overshooting. For example if the current precision
# is 3 and we want to get to 10, we're better off going up to 5
# instead of 6, because it's less work to get from 5 to 10 than it
# is to get from 6 to 10.
if ceil(log(M/target_k, 2)) == ceil(log(M/(target_k-1), 2)):
target_k -= 1
# temp = T^{p(3k-2)} 1/2 s x_i^3
temp = X.square() * (X * r)
# We know that the final result is only going to be correct mod
# p^(target_k), so we might as well truncate the extraneous terms now.
# temp = T^{p(k'-1)} 1/2 s x_i^3
temp = temp.shift(-p*(3*k - target_k - 1))
# X = T^{p(k'-1)} (3/2 x_i - 1/2 s x_i^3)
# = T^{p(k'-1)} x_{i+1}
X = (three_halves * X).shift(p*(target_k - k)) - temp
k = target_k
# Now k should equal M, since we're up to the correct precision
assert k == M, "Oops, something went wrong in the iteration"
# We should have s^{-1/2} correct to precision M.
# The following line can be uncommented to verify this.
# (It's a slow verification though, can double the whole computation time.)
#assert (p * X.square() * r * base_ring(2)).coeffs() == \
# R(p).shift(p*(2*M - 1)).coeffs()
# Finally incorporate frobenius of dx and x dx, and choose offset that
# compensates for our normalisations by powers of T.
F0 = base_ring(p) * x_to_p_less_one * X
F1 = F0 * x_to_p
offset = ((2*k-1)*p - 1)/2
return F0, F1, offset
def frobenius_expansion_by_series(Q, p, M):
r"""
Computes the action of Frobenius on `dx/y` and on `x dx/y`, using a
series expansion.
(This function computes the same thing as
frobenius_expansion_by_newton(), using a different method.
Theoretically the Newton method should be asymptotically faster,
when the precision gets large. However, in practice, this functions
seems to be marginally faster for moderate precision, so I'm
keeping it here until I figure out exactly why it's faster.)
(This function does *not* yet use the cohomology relations - that
happens afterwards in the "reduction" step.)
More specifically, it finds F0 and F1 in the quotient ring
`R[x, T]/(T - Q(x))`, such that
`F( dx/y) = T^{-r} F0 dx/y`, and
`F(x dx/y) = T^{-r} F1 dx/y` where
`r = ( (2M-3)p - 1 )/2`. (Here `T` is `y^2 = z^{-2}`,
and `R` is the coefficient ring of `Q`.)
`F_0` and `F_1` are computed in the
SpecialCubicQuotientRing associated to `Q`, so all powers
of `x^j` for `j \geq 3` are reduced to powers of
`T`.
It uses the sum
.. math::
F0 = \sum_{k=0}^{M-2} {-1/2 \choose k} p x^{p-1} E^k T^{(M-2-k)p}
and
.. math::
F1 = x^p F0,
where `E = Q(x^p) - Q(x)^p`.
INPUT:
- ``Q`` - cubic polynomial of the form
`Q(x) = x^3 + ax + b`, whose coefficient ring is a
`\ZZ/(p^M)\ZZ` -algebra
- ``p`` - residue characteristic of the `p`-adic field
- ``M`` - `p`-adic precision of the coefficient ring
(this will be used to determine the number of terms in the
series)
OUTPUT:
- ``F0, F1`` - elements of
SpecialCubicQuotientRing(Q), as described above
- ``r`` - non-negative integer, as described above
"""
S = SpecialCubicQuotientRing(Q)
x, _ = S.gens()
base_ring = S.base_ring()
x_to_p_less_1 = x**(p-1)
x_to_p = x_to_p_less_1 * x
# compute frobQ = Q(x^p)
x_to_p_squared = x_to_p * x_to_p
x_to_p_cubed = x_to_p_squared * x_to_p
frobQ = x_to_p_cubed + Q[1]*x_to_p + Q[0]*S(1)
# anticipating the day when p = 3 is supported:
# frobQ = x_to_p_cubed + Q[2]*x_to_p_squared + Q[1]*x_to_p + Q[0]*S(1)
E = frobQ - S(1).shift(p) # E = Q(x^p) - Q(x)^p
offset = int( ((2*M-3)*p-1)/2 )
term = p * x_to_p_less_1
F0 = term.shift((M-2)*p)
# todo: Possible speedup idea, perhaps by a factor of 2, but
# it requires a lot of work:
# Note that p divides E, so p^k divides E^k. So when we are
# working with high powers of E, we're doing a lot more work
# in the multiplications than we need to. To take advantage of
# this we would need some protocol for "lowering the precision"
# of a SpecialCubicQuotientRing. This would be quite messy to
# do properly over an arbitrary base ring. Perhaps it is
# feasible to do for the most common case (i.e. Z/p^nZ).
# (but it probably won't save much time unless p^n is very
# large, because the machine word size is probably pretty
# big anyway.)
for k in range(int(1), int(M-1)):
term = term * E
c = base_ring(binomial(-Integer(1)/2, k))
F0 += (term * c).shift((M-k-2)*p)
return F0, F0 * x_to_p, offset
def adjusted_prec(p, prec):
r"""
Computes how much precision is required in matrix_of_frobenius to
get an answer correct to prec `p`-adic digits.
The issue is that the algorithm used in matrix_of_frobenius
sometimes performs divisions by `p`, so precision is lost
during the algorithm.
The estimate returned by this function is based on Kedlaya's result
(Lemmas 2 and 3 of "Counting Points on Hyperelliptic Curves..."),
which implies that if we start with `M` `p`-adic
digits, the total precision loss is at most
`1 + \lfloor \log_p(2M-3) \rfloor` `p`-adic
digits. (This estimate is somewhat less than the amount you would
expect by naively counting the number of divisions by
`p`.)
INPUT:
- ``p`` - a prime = 5
- ``prec`` - integer, desired output precision, = 1
OUTPUT: adjusted precision (usually slightly more than prec)
"""
# initial estimate:
if prec <= 2:
adjusted = 2
else:
adjusted = prec + int(log(2*prec - 3, p)) - 1
# increase it until we have enough
while adjusted - int(log(2*adjusted - 3, p)) - 1 < prec:
adjusted += 1
return adjusted
def matrix_of_frobenius(Q, p, M, trace=None, compute_exact_forms=False):
"""
Computes the matrix of Frobenius on Monsky-Washnitzer cohomology,
with respect to the basis `(dx/y, x dx/y)`.
INPUT:
- ``Q`` - cubic polynomial `Q(x) = x^3 + ax + b`
defining an elliptic curve `E` by
`y^2 = Q(x)`. The coefficient ring of `Q` should be a
`\ZZ/(p^M)\ZZ`-algebra in which the matrix of
frobenius will be constructed.
- ``p`` - prime = 5 for which E has good reduction
- ``M`` - integer = 2; `p` -adic precision of
the coefficient ring
- ``trace`` - (optional) the trace of the matrix, if
known in advance. This is easy to compute because it's just the
`a_p` of the curve. If the trace is supplied,
matrix_of_frobenius will use it to speed the computation (i.e. we
know the determinant is `p`, so we have two conditions, so
really only column of the matrix needs to be computed. It's
actually a little more complicated than that, but that's the basic
idea.) If trace=None, then both columns will be computed
independently, and you can get a strong indication of correctness
by verifying the trace afterwards.
.. warning::
THE RESULT WILL NOT NECESSARILY BE CORRECT TO M p-ADIC
DIGITS. If you want prec digits of precision, you need to use
the function adjusted_prec(), and then you need to reduce the
answer mod `p^{\mathrm{prec}}` at the end.
OUTPUT:
2x2 matrix of frobenius on Monsky-Washnitzer cohomology,
with entries in the coefficient ring of Q.
EXAMPLES:
A simple example::
sage: p = 5
sage: prec = 3
sage: M = monsky_washnitzer.adjusted_prec(p, prec)
sage: M
5
sage: R.<x> = PolynomialRing(Integers(p**M))
sage: A = monsky_washnitzer.matrix_of_frobenius(x^3 - x + R(1/4), p, M)
sage: A
[3090 187]
[2945 408]
But the result is only accurate to prec digits::
sage: B = A.change_ring(Integers(p**prec))
sage: B
[90 62]
[70 33]
Check trace (123 = -2 mod 125) and determinant::
sage: B.det()
5
sage: B.trace()
123
sage: EllipticCurve([-1, 1/4]).ap(5)
-2
Try using the trace to speed up the calculation::
sage: A = monsky_washnitzer.matrix_of_frobenius(x^3 - x + R(1/4),
... p, M, -2)
sage: A
[2715 187]
[1445 408]
Hmmm... it looks different, but that's because the trace of our
first answer was only -2 modulo `5^3`, not -2 modulo
`5^5`. So the right answer is::
sage: A.change_ring(Integers(p**prec))
[90 62]
[70 33]
Check it works with only one digit of precision::
sage: p = 5
sage: prec = 1
sage: M = monsky_washnitzer.adjusted_prec(p, prec)
sage: R.<x> = PolynomialRing(Integers(p**M))
sage: A = monsky_washnitzer.matrix_of_frobenius(x^3 - x + R(1/4), p, M)
sage: A.change_ring(Integers(p))
[0 2]
[0 3]
Here's an example that's particularly badly conditioned for using
the trace trick::
sage: p = 11
sage: prec = 3
sage: M = monsky_washnitzer.adjusted_prec(p, prec)
sage: R.<x> = PolynomialRing(Integers(p**M))
sage: A = monsky_washnitzer.matrix_of_frobenius(x^3 + 7*x + 8, p, M)
sage: A.change_ring(Integers(p**prec))
[1144 176]
[ 847 185]
The problem here is that the top-right entry is divisible by 11,
and the bottom-left entry is divisible by `11^2`. So when
you apply the trace trick, neither `F(dx/y)` nor
`F(x dx/y)` is enough to compute the whole matrix to the
desired precision, even if you try increasing the target precision
by one. Nevertheless, ``matrix_of_frobenius`` knows
how to get the right answer by evaluating `F((x+1) dx/y)`
instead::
sage: A = monsky_washnitzer.matrix_of_frobenius(x^3 + 7*x + 8, p, M, -2)
sage: A.change_ring(Integers(p**prec))
[1144 176]
[ 847 185]
The running time is about ``O(p*prec**2)`` (times
some logarithmic factors), so it's feasible to run on fairly large
primes, or precision (or both?!?!)::
sage: p = 10007
sage: prec = 2
sage: M = monsky_washnitzer.adjusted_prec(p, prec)
sage: R.<x> = PolynomialRing(Integers(p**M))
sage: A = monsky_washnitzer.matrix_of_frobenius( # long time
... x^3 - x + R(1/4), p, M) # long time
sage: B = A.change_ring(Integers(p**prec)); B # long time
[74311982 57996908]
[95877067 25828133]
sage: B.det() # long time
10007
sage: B.trace() # long time
66
sage: EllipticCurve([-1, 1/4]).ap(10007) # long time
66
::
sage: p = 5
sage: prec = 300
sage: M = monsky_washnitzer.adjusted_prec(p, prec)
sage: R.<x> = PolynomialRing(Integers(p**M))
sage: A = monsky_washnitzer.matrix_of_frobenius( # long time
... x^3 - x + R(1/4), p, M) # long time
sage: B = A.change_ring(Integers(p**prec)) # long time
sage: B.det() # long time
5
sage: -B.trace() # long time
2
sage: EllipticCurve([-1, 1/4]).ap(5) # long time
-2
Let's check consistency of the results for a range of precisions::
sage: p = 5
sage: max_prec = 60
sage: M = monsky_washnitzer.adjusted_prec(p, max_prec)
sage: R.<x> = PolynomialRing(Integers(p**M))
sage: A = monsky_washnitzer.matrix_of_frobenius(x^3 - x + R(1/4), p, M) # long time
sage: A = A.change_ring(Integers(p**max_prec)) # long time
sage: result = [] # long time
sage: for prec in range(1, max_prec): # long time
... M = monsky_washnitzer.adjusted_prec(p, prec) # long time
... R.<x> = PolynomialRing(Integers(p^M),'x') # long time
... B = monsky_washnitzer.matrix_of_frobenius( # long time
... x^3 - x + R(1/4), p, M) # long time
... B = B.change_ring(Integers(p**prec)) # long time
... result.append(B == A.change_ring( # long time
... Integers(p**prec))) # long time
sage: result == [True] * (max_prec - 1) # long time
True
The remaining examples discuss what happens when you take the
coefficient ring to be a power series ring; i.e. in effect you're
looking at a family of curves.
The code does in fact work...
::
sage: p = 11
sage: prec = 3
sage: M = monsky_washnitzer.adjusted_prec(p, prec)
sage: S.<t> = PowerSeriesRing(Integers(p**M), default_prec=4)
sage: a = 7 + t + 3*t^2
sage: b = 8 - 6*t + 17*t^2
sage: R.<x> = PolynomialRing(S)
sage: Q = x**3 + a*x + b
sage: A = monsky_washnitzer.matrix_of_frobenius(Q, p, M) # long time
sage: B = A.change_ring(PowerSeriesRing(Integers(p**prec), 't', default_prec=4)) # long time
sage: B # long time
[1144 + 264*t + 841*t^2 + 1025*t^3 + O(t^4) 176 + 1052*t + 216*t^2 + 523*t^3 + O(t^4)]
[ 847 + 668*t + 81*t^2 + 424*t^3 + O(t^4) 185 + 341*t + 171*t^2 + 642*t^3 + O(t^4)]
The trace trick should work for power series rings too, even in the
badly- conditioned case. Unfortunately I don't know how to compute
the trace in advance, so I'm not sure exactly how this would help.
Also, I suspect the running time will be dominated by the
expansion, so the trace trick won't really speed things up anyway.
Another problem is that the determinant is not always p::
sage: B.det() # long time
11 + 484*t^2 + 451*t^3 + O(t^4)
However, it appears that the determinant always has the property
that if you substitute t - 11t, you do get the constant series p
(mod p\*\*prec). Similarly for the trace. And since the parameter
only really makes sense when it's divisible by p anyway, perhaps
this isn't a problem after all.
"""
M = int(M)
if M < 2:
raise ValueError, "M (=%s) must be at least 2" % M
base_ring = Q.base_ring()
# Expand out frobenius of dx/y and x dx/y.
# (You can substitute frobenius_expansion_by_series here, that will work
# as well. See its docstring for some performance notes.)
F0, F1, offset = frobenius_expansion_by_newton(Q, p, M)
#F0, F1, offset = frobenius_expansion_by_series(Q, p, M)
if compute_exact_forms:
# we need to do all the work to get the exact expressions f such that F(x^i dx/y) = df + \sum a_i x^i dx/y
F0_coeffs = transpose_list(F0.coeffs())
F0_reduced, f_0 = reduce_all(Q, p, F0_coeffs, offset, True)
F1_coeffs = transpose_list(F1.coeffs())
F1_reduced, f_1 = reduce_all(Q, p, F1_coeffs, offset, True)
elif M == 2:
# This implies that only one digit of precision is valid, so we only need
# to reduce the second column. Also, the trace doesn't help at all.
F0_reduced = [ base_ring(0), base_ring(0) ]
F1_coeffs = transpose_list(F1.coeffs())
F1_reduced = reduce_all(Q, p, F1_coeffs, offset)
elif trace is None:
# No trace provided, just reduce F(dx/y) and F(x dx/y) separately.
F0_coeffs = transpose_list(F0.coeffs())
F0_reduced = reduce_all(Q, p, F0_coeffs, offset)
F1_coeffs = transpose_list(F1.coeffs())
F1_reduced = reduce_all(Q, p, F1_coeffs, offset)
else:
# Trace has been provided.
# In most cases this can be used to quickly compute F(dx/y) from
# F(x dx/y). However, if we're unlucky, the (dx/y)-component of
# F(x dx/y) (i.e. the top-right corner of the matrix) may be divisible
# by p, in which case there isn't enough information to get the
# (x dx/y)-component of F(dx/y) to the desired precision. When this
# happens, it turns out that F((x+1) dx/y) always *does* give enough
# information (together with the trace) to get both columns to the
# desired precision.
# First however we need a quick way of telling whether the top-right
# corner is divisible by p, i.e. we want to compute the second column
# of the matrix mod p. We could do this by just running the entire
# algorithm with M = 2 (which assures precision 1). Luckily, we've
# already done most of the work by computing F1 to high precision; so
# all we need to do is extract the coefficients that would correspond
# to the first term of the series, and run the reduction on them.
# todo: actually we only need to do this reduction step mod p^2, not
# mod p^M, which is what the code currently does. If the base ring
# is Integers(p^M), then it's easy. Otherwise it's tricky to construct
# the right ring, I don't know how to do it.
F1_coeffs = transpose_list(F1.coeffs())
F1_modp_coeffs = F1_coeffs[int((M-2)*p):]
# make a copy, because reduce_all will destroy the coefficients:
F1_modp_coeffs = [[cell for cell in row] for row in F1_modp_coeffs]
F1_modp_offset = offset - (M-2)*p
F1_modp_reduced = reduce_all(Q, p, F1_modp_coeffs, F1_modp_offset)
if F1_modp_reduced[0].is_unit():
# If the first entry is invertible mod p, then F(x dx/y) is sufficient
# to get the whole matrix.
F1_reduced = reduce_all(Q, p, F1_coeffs, offset)
F0_reduced = [ base_ring(trace) - F1_reduced[1], None ]
# using that the determinant is p:
F0_reduced[1] = (F0_reduced[0] * F1_reduced[1] - base_ring(p)) \
/ F1_reduced[0]
else:
# If the first entry is zero mod p, then F((x+1) dx/y) will be sufficient
# to get the whole matrix. (Here we are using the fact that the second
# entry *cannot* be zero mod p. This is guaranteed by some results in
# section 3.2 of ``Computation of p-adic Heights and Log Convergence''
# by Mazur, Stein, Tate. But let's quickly check it anyway :-))
assert F1_modp_reduced[1].is_unit(), \
"Hey that's impossible! The second entry in the second column " \
"should be invertible mod p!"
G0_coeffs = transpose_list( (F0 + F1).coeffs())
G0_reduced = reduce_all(Q, p, G0_coeffs, offset)
# Now G0_reduced expresses F((x+1) dx/y) in terms of dx/y and x dx/y.
# Re-express this in terms of (x+1) dx/y and x dx/y.
H0_reduced = [ G0_reduced[0], G0_reduced[1] - G0_reduced[0] ]
# The thing we're about to divide by better be a unit.
assert H0_reduced[1].is_unit(), \
"Hey that's impossible! The second entry in this column " \
"should be invertible mod p!"
# Figure out the second column using the trace...
H1_reduced = [ None, base_ring(trace) - H0_reduced[0] ]
# ... and using that the determinant is p:
H1_reduced[0] = (H0_reduced[0] * H1_reduced[1] - base_ring(p)) \
/ H0_reduced[1]
# Finally, change back to the usual basis (dx/y, x dx/y)
F1_reduced = [ H1_reduced[0], \
H1_reduced[0] + H1_reduced[1] ]
F0_reduced = [ H0_reduced[0] - F1_reduced[0],
H0_reduced[0] + H0_reduced[1] - F1_reduced[1] ]
# One more sanity check: our final result should be congruent mod p
# to the approximation we used earlier.
assert not (
(F1_reduced[0] - F1_modp_reduced[0]).is_unit() or \
(F1_reduced[1] - F1_modp_reduced[1]).is_unit() or \
F0_reduced[0].is_unit() or F0_reduced[1].is_unit()), \
"Hey that's impossible! The output matrix is not congruent mod p " \
"to the approximation found earlier!"
if compute_exact_forms:
return matrix(base_ring, 2, 2, [F0_reduced[0], F1_reduced[0],
F0_reduced[1], F1_reduced[1]]), f_0, f_1
else:
return matrix(base_ring, 2, 2, [F0_reduced[0], F1_reduced[0],
F0_reduced[1], F1_reduced[1]])
#*****************************************************************************
# This is a generalization of the above functionality for hyperelliptic curves.
#
# THIS IS A WORK IN PROGRESS.
#
# I tried to embed must stuff into the rings themselves rather than
# just extract and manipulate lists of coefficients. Hence the implementations
# below are much less optimized, so are much slower, but should hopefully be
# easier to follow. (E.g. one can print/make sense of intermediate results.)
#
# AUTHOR:
# -- Robert Bradshaw (2007-04)
#
#*****************************************************************************
import weakref
from sage.schemes.hyperelliptic_curves.all import is_HyperellipticCurve, HyperellipticCurve
from sage.rings.padics.all import pAdicField
from sage.rings.all import QQ, is_LaurentSeriesRing, is_IntegralDomain
from sage.modules.all import FreeModule, is_FreeModuleElement
from sage.misc.profiler import Profiler
from sage.misc.misc import repr_lincomb
def matrix_of_frobenius_hyperelliptic(Q, p=None, prec=None, M=None):
"""
Computes the matrix of Frobenius on Monsky-Washnitzer cohomology,
with respect to the basis `(dx/2y, x dx/2y, ...x^{d-2} dx/2y)`, where
`d` is the degree of `Q`.
INPUT:
- ``Q`` - monic polynomial `Q(x)`
- ``p`` - prime `\geq 5` for which `E` has good reduction
- ``prec`` - (optional) `p`-adic precision of the coefficient ring
- ``M`` - (optional) adjusted `p`-adic precision of the coefficint ring
OUTPUT:
`(d-1)` x `(d-1)` matrix `M` of Frobenius on Monsky-Washnitzer cohomology,
and list of differentials \{f_i \} such that
.. math::
\phi^* (x^i dx/2y) = df_i + M[i]*vec(dx/2y, ..., x^{d-2} dx/2y)
EXAMPLES::
sage: p = 5
sage: prec = 3
sage: R.<x> = QQ['x']
sage: A,f = monsky_washnitzer.matrix_of_frobenius_hyperelliptic(x^5 - 2*x + 3, p, prec)
sage: A
[ 4*5 + O(5^3) 5 + 2*5^2 + O(5^3) 2 + 3*5 + 2*5^2 + O(5^3) 2 + 5 + 5^2 + O(5^3)]
[ 3*5 + 5^2 + O(5^3) 3*5 + O(5^3) 4*5 + O(5^3) 2 + 5^2 + O(5^3)]
[ 4*5 + 4*5^2 + O(5^3) 3*5 + 2*5^2 + O(5^3) 5 + 3*5^2 + O(5^3) 2*5 + 2*5^2 + O(5^3)]
[ 5^2 + O(5^3) 5 + 4*5^2 + O(5^3) 4*5 + 3*5^2 + O(5^3) 2*5 + O(5^3)]
"""
prof = Profiler()
prof("setup")
if p is None:
try:
K = Q.base_ring()
p = K.prime()
prec = K.precision_cap()
except AttributeError:
raise ValueError, "p and prec must be specified if Q is not defined over a p-adic ring"
if M is None:
M = adjusted_prec(p, prec)
extra_prec_ring = Integers(p**M)
# extra_prec_ring = pAdicField(p, M) # SLOW!
real_prec_ring = pAdicField(p, prec) # pAdicField(p, prec) # To capped absolute?
S = SpecialHyperellipticQuotientRing(Q, extra_prec_ring, True)
MW = S.monsky_washnitzer()
prof("frob basis elements")
F = MW.frob_basis_elements(M, p)
prof("rationalize")
# do reduction over Q in case we have non-integral entries (and it's so much faster than padics)
rational_S = S.change_ring(QQ)
# this is a hack until pAdics are fast
# (They are in the latest development bundle, but its not standard and I'd need to merge.
# (it will periodically cast into this ring to reduce coefficient size)
rational_S._prec_cap = p**M
rational_S._p = p
# S._p = p
# rational_S(F[0]).reduce_fast()
# prof("reduce others")
# rational_S = S.change_ring(pAdicField(p, M))
F = [rational_S(F_i) for F_i in F]
prof("reduce")
reduced = [F_i.reduce_fast(True) for F_i in F]
# reduced = [F_i.reduce() for F_i in F]
#print reduced[0][0].diff() - F[0]
# but the coeffs are WAY more precision than they need to be
# print reduced[0][1]
prof("make matrix")
# now take care of precision capping
M = matrix(real_prec_ring, [a for f, a in reduced])
for i in range(M.ncols()):
for j in range(M.nrows()):
M[i,j] = M[i,j].add_bigoh(prec)
# print prof
return M.transpose(), [f for f, a in reduced]
# For uniqueness (as many of the non-trivial calculations are cached along the way).
_special_ring_cache = {}
_mw_cache = {}
def SpecialHyperellipticQuotientRing(*args):
if _special_ring_cache.has_key(args):
R = _special_ring_cache[args]()
if R is not None:
return R
R = SpecialHyperellipticQuotientRing_class(*args)
_special_ring_cache[args] = weakref.ref(R)
return R
def MonskyWashnitzerDifferentialRing(base_ring):
if _mw_cache.has_key(base_ring):
R = _mw_cache[base_ring]()
if R is not None:
return R
R = MonskyWashnitzerDifferentialRing_class(base_ring)
_mw_cache[base_ring] = weakref.ref(R)
return R
class SpecialHyperellipticQuotientRing_class(CommutativeAlgebra):
_p = None
def __init__(self, Q, R=None, invert_y=True):
if R is None:
R = Q.base_ring()
# Trac ticket #9138: CommutativeAlgebra.__init__ must not be
# done so early. It tries to register a coercion, but that
# requires the hash bein available. But the hash, in its
# default implementation, relies on the string representation,
# which is not available at this point.
#CommutativeAlgebra.__init__(self, R) # moved to below.
x = PolynomialRing(R, 'xx').gen(0)
if is_EllipticCurve(Q):
E = Q
if E.a1() != 0 or E.a2() != 0:
raise NotImplementedError, "Curve must be in Weierstrass normal form."
Q = -E.change_ring(R).defining_polynomial()(x,0,1)
self._curve = E
elif is_HyperellipticCurve(Q):
C = Q
if C.hyperelliptic_polynomials()[1] != 0:
raise NotImplementedError, "Curve must be of form y^2 = Q(x)."
Q = C.hyperelliptic_polynomials()[0].change_ring(R)
self._curve = C
if is_Polynomial(Q):
self._Q = Q.change_ring(R)
self._coeffs = self._Q.coeffs()
if self._coeffs.pop() != 1:
raise NotImplementedError, "Polynomial must be monic."
if not hasattr(self, '_curve'):
if self._Q.degree() == 3:
ainvs = [0, self._Q[2], 0, self._Q[1], self._Q[0]]
self._curve = EllipticCurve(ainvs)
else:
self._curve = HyperellipticCurve(self._Q)
else:
raise NotImplementedError, "Must be an elliptic curve or polynomial Q for y^2 = Q(x)\n(Got element of %s)" % Q.parent()
self._n = degree = int(Q.degree())
self._series_ring = (LaurentSeriesRing if invert_y else PolynomialRing)(R, 'y')
self._series_ring_y = self._series_ring.gen(0)
self._series_ring_0 = self._series_ring(0)
# Trac ticket #9138: Initialise the commutative algebra here!
# Below, we do self(self._poly_ring.gen(0)), which requires
# the initialisation being finished.
CommutativeAlgebra.__init__(self, R)
self._poly_ring = PolynomialRing(self._series_ring, 'x')
self._x = self(self._poly_ring.gen(0))
self._y = self(self._series_ring.gen(0))
self._Q_coeffs = Q.change_ring(self._series_ring).list()
self._dQ = Q.derivative().change_ring(self)(self._x)
self._monsky_washnitzer = MonskyWashnitzerDifferentialRing(self)
self._monomial_diffs = {}
self._monomial_diff_coeffs = {}
def _repr_(self):
y_inverse = ",y^-1" if is_LaurentSeriesRing(self._series_ring) else ""
return "SpecialHyperellipticQuotientRing K[x,y%s] / (y^2 = %s) over %s"%(y_inverse, self._Q, self.base_ring())
def base_extend(self, R):
if R.has_coerce_map_from(self.base_ring()):
return self.change_ring(R)
else:
raise TypeError, "no such base extension"
def change_ring(self, R):
return SpecialHyperellipticQuotientRing(self._Q, R, is_LaurentSeriesRing(self._series_ring))
def __call__(self, val, offset=0, check=True):
if isinstance(val, SpecialHyperellipticQuotientElement) and val.parent() is self:
if offset == 0:
return val
else:
return val << offset
elif isinstance(val, MonskyWashnitzerDifferential):
return self._monsky_washnitzer(val)
return SpecialHyperellipticQuotientElement(self, val, offset, check)
def gens(self):
return self._x, self._y
def x(self):
return self._x
def y(self):
return self._y
def monomial(self, i, j, b=None):
"""
Returns `b y^j x^i`, computed quickly.
"""
i = int(i)
j = int(j)
if 0 < i and i < self._n:
if b is None:
by_to_j = self._series_ring_y << (j-1)
else:
by_to_j = self._series_ring(b) << j
v = [self._series_ring_0] * self._n
v[i] = by_to_j
return self(v)
else:
return (self._x ** i) << j if b is None else self.base_ring()(b) * (self._x ** i) << j
def monomial_diff_coeffs(self, i, j):
r"""
The key here is that the formula for `d(x^iy^j)` is messy
in terms of `i`, but varies nicely with `j`.
.. math::
d(x^iy^j) = y^{j-1} (2ix^{i-1}y^2 + j (A_i(x) + B_i(x)y^2)) \frac{dx}{2y}
Where `A,B` have degree at most `n-1` for each
`i`. Pre-compute `A_i, B_i` for each `i`
the "hard" way, and the rest are easy.
"""
try:
return self._monomial_diff_coeffs[i,j]
except KeyError:
pass
if i < self._n:
try:
A, B, two_i_x_to_i = self._precomputed_diff_coeffs[i]
except AttributeError:
self._precomputed_diff_coeffs = self._precompute_monomial_diffs()
A, B, two_i_x_to_i = self._precomputed_diff_coeffs[i]
if i == 0:
return j*A, j*B
else:
return j*A, j*B + two_i_x_to_i
else:
dg = self.monomial(i, j).diff()
coeffs = [dg.extract_pow_y(j-1), dg.extract_pow_y(j+1)]
self._monomial_diff_coeffs[i,j] = coeffs
return coeffs
def monomial_diff_coeffs_matrices(self):
self.monomial_diff_coeffs(0, 0) # precompute stuff
R = self.base_ring()
mat_1 = matrix(R, self._n, self._n)
mat_2 = matrix(R, self._n, self._n)
for i in range(self._n):
mat_1[i] = self._precomputed_diff_coeffs[i][1]
mat_2[i] = self._precomputed_diff_coeffs[i][2]
return mat_1.transpose(), mat_2.transpose()
def _precompute_monomial_diffs(self):
x, y = self.gens()
R = self.base_ring()
V = FreeModule(R, self.degree())
As = []
for i in range(self.degree()):
dg = self.monomial(i, 1).diff()
two_i_x_to_i = R(2*i) * x**(i-1) * y*y if i > 0 else self(0)
A = dg - self._monsky_washnitzer(two_i_x_to_i)
As.append( (V(A.extract_pow_y(0)), V(A.extract_pow_y(2)), V(two_i_x_to_i.extract_pow_y(2))) )
return As
def Q(self):
return self._Q
def curve(self):
return self._curve
def degree(self):
return self._n
def prime(self):
return self._p
def monsky_washnitzer(self):
return self._monsky_washnitzer
def is_field(self, proof = True):
return False
class SpecialHyperellipticQuotientElement(CommutativeAlgebraElement):
def __init__(self, parent, val=0, offset=0, check=True):
CommutativeAlgebraElement.__init__(self, parent)
if not check:
self._f = parent._poly_ring(val, check=False)
return
if isinstance(val, SpecialHyperellipticQuotientElement):
R = parent.base_ring()
self._f = parent._poly_ring([a.change_ring(R) for a in val._f])
return
if isinstance(val, tuple):
val, offset = val
if isinstance(val, list) and len(val) > 0 and is_FreeModuleElement(val[0]):
val = transpose_list(val)
self._f = parent._poly_ring(val)
if offset != 0:
self._f = self._f.parent()([a << offset for a in self._f], check=False)
def __cmp__(self, other):
"""
EXAMPLES:
"""
return cmp(self._f, other._f)
def change_ring(self, R):
return self.parent().change_ring(R)(self)
def __call__(self, *x):
return self._f(*x)
def __invert__(self):
"""
The general element in our ring is not invertible, but `y` may be. We
do not want to pass to the fraction field.
"""
if self._f.degree() == 0 and self._f[0].is_unit():
return SpecialHyperellipticQuotientElement(self.parent(), ~self._f[0])
else:
raise ZeroDivisionError, "Element not invertible"
def __nonzero__(self):
return not not self._f
def __eq__(self, other):
if not isinstance(other, SpecialHyperellipticQuotientElement):
other = self.parent()(other)
return self._f == other._f
def _add_(self, other):
return SpecialHyperellipticQuotientElement(self.parent(), self._f + other._f)
def _sub_(self, other):
return SpecialHyperellipticQuotientElement(self.parent(), self._f - other._f)
def _mul_(self, other):
# over laurent series, addition and subtraction can be expensive,
# and the degree of this poly is small enough that Karatsuba actually hurts
# significantly in some cases
if self._f[0].valuation() + other._f[0].valuation() > -200:
prod = self._f._mul_generic(other._f)
else:
prod = self._f * other._f
v = prod.list()
parent = self.parent()
Q_coeffs = parent._Q_coeffs
n = len(Q_coeffs) - 1
y2 = self.parent()._series_ring_y << 1
for i in range(len(v)-1, n-1, -1):
for j in range(n):
v[i-n+j] -= Q_coeffs[j] * v[i]
v[i-n] += y2 * v[i]
return SpecialHyperellipticQuotientElement(parent, v[0:n])
def _rmul_(self, c):
coeffs = self._f.list(copy=False)
return self.parent()([c*a for a in coeffs], check=False)
def _lmul_(self, c):
coeffs = self._f.list(copy=False)
return self.parent()([a*c for a in coeffs], check=False)
def __lshift__(self, k):
coeffs = self._f.list(copy=False)
return self.parent()([a << k for a in coeffs], check=False)
def __rshift__(self, k):
coeffs = self._f.list(copy=False)
return self.parent()([a >> k for a in coeffs], check=False)
def truncate_neg(self, n):
coeffs = self._f.list(copy=False)
return self.parent()([a.truncate_neg(n) for a in coeffs], check=False)
def _repr_(self):
x = PolynomialRing(QQ, 'x').gen(0)
coeffs = self._f.list()
return repr_lincomb([(x**i, coeffs[i]) for i in range(len(coeffs))])
def _latex_(self):
x = PolynomialRing(QQ, 'x').gen(0)
coeffs = self._f.list()
return repr_lincomb([(x**i, coeffs[i]) for i in range(len(coeffs))], is_latex=True)
def diff(self):
# try:
# return self._diff_x
# except AttributeError:
# pass
# d(self) = A dx + B dy
# = (2y A + BQ') dx/2y
parent = self.parent()
R = parent.base_ring()
x, y = parent.gens()
v = self._f.list()
n = len(v)
A = parent([R(i) * v[i] for i in range(1,n)])
B = parent([a.derivative() for a in v])
dQ = parent._dQ
return parent._monsky_washnitzer( (R(2) * A << 1) + dQ * B )
# self._diff = self.parent()._monsky_washnitzer( two_y * A + dQ * B )
# return self._diff
def extract_pow_y(self, k):
v = [a[k] for a in self._f.list()]
while len(v) < self.parent()._n:
v.append(0)
return v
def min_pow_y(self):
if self._f.degree() == -1:
return 0
return min([a.valuation() for a in self._f.list()])
def max_pow_y(self):
if self._f.degree() == -1:
return 0
return max([a.degree() for a in self._f.list()])
def coeffs(self, R=None):
"""
Returns the raw coefficients of this element.
INPUT:
- ``R`` - an (optional) base-ring in which to cast the coefficients
OUTPUT:
- ``coeffs`` - a list of coefficients of powers of `x` for each power
of `y`
- ``n`` - an offset indicating the power of `y` of the first list
element
EXAMPLES::
sage: R.<x> = QQ['x']
sage: E = HyperellipticCurve(x^5-3*x+1)
sage: x,y = E.monsky_washnitzer_gens()
sage: x.coeffs()
([(0, 1, 0, 0, 0)], 0)
sage: y.coeffs()
([(0, 0, 0, 0, 0), (1, 0, 0, 0, 0)], 0)
sage: a = sum(n*x^n for n in range(5)); a
x + 2*x^2 + 3*x^3 + 4*x^4
sage: a.coeffs()
([(0, 1, 2, 3, 4)], 0)
sage: a.coeffs(Qp(7))
([(0, 1 + O(7^20), 2 + O(7^20), 3 + O(7^20), 4 + O(7^20))], 0)
sage: (a*y).coeffs()
([(0, 0, 0, 0, 0), (0, 1, 2, 3, 4)], 0)
sage: (a*y^-2).coeffs()
([(0, 1, 2, 3, 4), (0, 0, 0, 0, 0), (0, 0, 0, 0, 0)], -2)
Note that the coefficient list is transposed compared to how they
are stored and printed::
sage: a*y^-2
(y^-2)*x + (2*y^-2)*x^2 + (3*y^-2)*x^3 + (4*y^-2)*x^4
A more complicated example::
sage: a = x^20*y^-3 - x^11*y^2; a
(y^-3-4*y^-1+6*y-4*y^3+y^5)*1 + (-12*y^-3+36*y^-1-36*y-y^2+12*y^3+2*y^4-y^6)*x + (54*y^-3-108*y^-1+54*y+6*y^2-6*y^4)*x^2 + (-108*y^-3+108*y^-1-9*y^2)*x^3 + (81*y^-3)*x^4
sage: raw, offset = a.coeffs()
sage: a.min_pow_y()
-3
sage: offset
-3
sage: raw
[(1, -12, 54, -108, 81),
(0, 0, 0, 0, 0),
(-4, 36, -108, 108, 0),
(0, 0, 0, 0, 0),
(6, -36, 54, 0, 0),
(0, -1, 6, -9, 0),
(-4, 12, 0, 0, 0),
(0, 2, -6, 0, 0),
(1, 0, 0, 0, 0),
(0, -1, 0, 0, 0)]
sage: sum(c * x^i * y^(j+offset) for j, L in enumerate(raw) for i, c in enumerate(L)) == a
True
Can also be used to construct elements::
sage: a.parent()(raw, offset) == a
True
"""
zero = self.base_ring()(0) if R is None else R(0)
y_offset = min(self.min_pow_y(), 0)
y_degree = max(self.max_pow_y(), 0)
coeffs = []
n = y_degree - y_offset + 1
for a in self._f.list():
k = a.valuation() - y_offset
z = a.list()
coeffs.append( [zero] * k + z + [zero]*(n - len(z) - k))
while len(coeffs) < self.parent().degree():
coeffs.append( [zero] * n )
V = FreeModule(self.base_ring() if R is None else R, self.parent().degree())
coeffs = transpose_list(coeffs)
return [V(a) for a in coeffs], y_offset
class MonskyWashnitzerDifferentialRing_class(Module):
def __init__(self, base_ring):
r"""
Class for the ring of Monsky--Washnitzer differentials over a given
base ring.
"""
Module.__init__(self, base_ring)
self._cache = {}
def invariant_differential(self):
"""
Returns `dx/2y` as an element of self.
EXAMPLES::
sage: R.<x> = QQ['x']
sage: C = HyperellipticCurve(x^5-4*x+4)
sage: MW = C.invariant_differential().parent()
sage: MW.invariant_differential()
1 dx/2y
"""
return self(1)
def __call__(self, val, offset=0):
return MonskyWashnitzerDifferential(self, val, offset)
def base_extend(self, R):
"""
Return a new differential ring which is self base-extended to `R`
INPUT:
- ``R`` - ring
OUTPUT:
Self, base-extended to `R`.
EXAMPLES::
sage: R.<x> = QQ['x']
sage: C = HyperellipticCurve(x^5-4*x+4)
sage: MW = C.invariant_differential().parent()
sage: MW.base_ring()
SpecialHyperellipticQuotientRing K[x,y,y^-1] / (y^2 = x^5 - 4*x + 4) over Rational Field
sage: MW.base_extend(Qp(5,5)).base_ring()
SpecialHyperellipticQuotientRing K[x,y,y^-1] / (y^2 = (1 + O(5^5))*x^5 + (1 + 4*5 + 4*5^2 + 4*5^3 + 4*5^4 + O(5^5))*x + (4 + O(5^5))) over 5-adic Field with capped relative precision 5
"""
return MonskyWashnitzerDifferentialRing(self.base_ring().base_extend(R))
def change_ring(self, R):
"""
Returns a new differential ring which is self with the coefficient
ring changed to `R`.
INPUT:
- ``R`` - ring of coefficients
OUTPUT:
Self, with the coefficient ring changed to `R`.
EXAMPLES::
sage: R.<x> = QQ['x']
sage: C = HyperellipticCurve(x^5-4*x+4)
sage: MW = C.invariant_differential().parent()
sage: MW.base_ring()
SpecialHyperellipticQuotientRing K[x,y,y^-1] / (y^2 = x^5 - 4*x + 4) over Rational Field
sage: MW.change_ring(Qp(5,5)).base_ring()
SpecialHyperellipticQuotientRing K[x,y,y^-1] / (y^2 = (1 + O(5^5))*x^5 + (1 + 4*5 + 4*5^2 + 4*5^3 + 4*5^4 + O(5^5))*x + (4 + O(5^5))) over 5-adic Field with capped relative precision 5
"""
return MonskyWashnitzerDifferentialRing(self.base_ring().change_ring(R))
def degree(self):
"""
Returns the degree of `Q(x)`, where the model of the underlying
hyperelliptic curve of self is given by `y^2 = Q(x)`.
EXAMPLES::
sage: R.<x> = QQ['x']
sage: C = HyperellipticCurve(x^5-4*x+4)
sage: MW = C.invariant_differential().parent()
sage: MW.Q()
x^5 - 4*x + 4
sage: MW.degree()
5
"""
return self.base_ring().degree()
def dimension(self):
"""
Returns the dimension of self.
EXAMPLES::
sage: R.<x> = QQ['x']
sage: C = HyperellipticCurve(x^5-4*x+4)
sage: K = Qp(7,5)
sage: CK = C.change_ring(K)
sage: MW = CK.invariant_differential().parent()
sage: MW.dimension()
4
"""
return self.base_ring().degree()-1
def Q(self):
"""
Returns `Q(x)` where the model of the underlying hyperelliptic curve
of self is given by `y^2 = Q(x)`.
EXAMPLES::
sage: R.<x> = QQ['x']
sage: C = HyperellipticCurve(x^5-4*x+4)
sage: MW = C.invariant_differential().parent()
sage: MW.Q()
x^5 - 4*x + 4
"""
return self.base_ring().Q()
def x_to_p(self, p):
"""
Returns and caches `x^p`, reduced via the relations coming from the
defining polynomial of the hyperelliptic curve.
EXAMPLES::
sage: R.<x> = QQ['x']
sage: C = HyperellipticCurve(x^5-4*x+4)
sage: MW = C.invariant_differential().parent()
sage: MW.x_to_p(3)
x^3
sage: MW.x_to_p(5)
(-4+y^2)*1 + 4*x
sage: MW.x_to_p(101) is MW.x_to_p(101)
True
"""
try:
return self._cache["x_to_p", p]
except KeyError:
x_to_p = self.base_ring().x() ** p
self._cache["x_to_p", p] = x_to_p
return x_to_p
def frob_Q(self, p):
"""
Returns and caches `Q(x^p)`, which is used in computing the image of
`y` under a `p`-power lift of Frobenius to `A^{\dagger}`.
EXAMPLES::
sage: R.<x> = QQ['x']
sage: C = HyperellipticCurve(x^5-4*x+4)
sage: MW = C.invariant_differential().parent()
sage: MW.frob_Q(3)
(-60+48*y^2-12*y^4+y^6)*1 + (192-96*y^2+12*y^4)*x + (-192+48*y^2)*x^2 + 60*x^3
sage: MW.Q()(MW.x_to_p(3))
(-60+48*y^2-12*y^4+y^6)*1 + (192-96*y^2+12*y^4)*x + (-192+48*y^2)*x^2 + 60*x^3
sage: MW.frob_Q(11) is MW.frob_Q(11)
True
"""
try:
return self._cache["frobQ", p]
except KeyError:
x_to_p = self.x_to_p(p)
frobQ = self.base_ring()._Q.change_ring(self.base_ring())(x_to_p)
self._cache["frobQ", p] = frobQ
return frobQ
def frob_invariant_differential(self, prec, p):
r"""
Kedlaya's algorithm allows us to calculate the action of Frobenius on
the Monsky-Washnitzer cohomology. First we lift `\phi` to `A^{\dagger}`
by setting
.. math::
\phi(x) = x^p
\phi(y) = y^p \sqrt{1 + \frac{Q(x^p) - Q(x)^p}{Q(x)^p}}.
Pulling back the differential `dx/2y`, we get
.. math::
\phi^*(dx/2y) = px^{p-1} y(\phi(y))^{-1} dx/2y
= px^{p-1} y^{1-p} \sqrt{1+ \frac{Q(x^p) - Q(x)^p}{Q(x)^p}} dx/2y
Use Newton's method to calculate the square root.
EXAMPLES::
sage: R.<x> = QQ['x']
sage: C = HyperellipticCurve(x^5-4*x+4)
sage: prec = 2
sage: p = 7
sage: MW = C.invariant_differential().parent()
sage: MW.frob_invariant_differential(prec,p)
((67894400*y^-20-81198880*y^-18+40140800*y^-16-10035200*y^-14+1254400*y^-12-62720*y^-10)*1 + (-119503944*y^-20+116064242*y^-18-43753472*y^-16+7426048*y^-14-514304*y^-12+12544*y^-10-1568*y^-8+70*y^-6+7*y^-4)*x + (78905288*y^-20-61014016*y^-18+16859136*y^-16-2207744*y^-14+250880*y^-12-37632*y^-10+3136*y^-8-70*y^-6)*x^2 + (-39452448*y^-20+26148752*y^-18-8085490*y^-16+2007040*y^-14-376320*y^-12+37632*y^-10-1568*y^-8)*x^3 + (21102144*y^-20-18120592*y^-18+8028160*y^-16-2007040*y^-14+250880*y^-12-12544*y^-10)*x^4) dx/2y
"""
prof = Profiler()
prof("setup")
# TODO, would it be useful to be able to take Frobenius of any element? Less efficient?
x, y = self.base_ring().gens()
prof("x_to_p")
x_to_p_less_1 = x**(p-1)
x_to_p = x*x_to_p_less_1
# cache for future use
self._cache["x_to_p", p] = x_to_p
prof("frob_Q")
a = self.frob_Q(p) >> 2*p # frobQ * y^{-2p}
prof("sqrt")
Q = self.base_ring()._Q
# three_halves = Q.parent().base_ring()(Rational((3,2)))
# one_half = Q.parent().base_ring()(Rational((1,2)))
three_halves = self.base_ring()._series_ring.base_ring()(Rational((3,2)))
one_half = self.base_ring()._series_ring.base_ring()(Rational((1,2)))
half_a = a._rmul_(one_half)
# We are solving for t = a^{-1/2} = (F_pQ y^{-p})^{-1/2}
# Newton's method converges because we know the root is in the same residue class as 1.
# t = self.base_ring()(1)
t = self.base_ring()(three_halves) - half_a # first iteration trivial, start with prec 2
for cur_prec in newton_method_sizes(prec)[2:]: # newton_method_sizes = [1, 2, ...]
y_prec = -(2*cur_prec-1)*p+1 # binomial expansion is $\sum p^{k+1} y^{-(2k+1)p+1} f(x)$
# so if we are only correct mod p^prec, can ignore y powers less than y_prec
t_cube = (t*t*t).truncate_neg(y_prec)
t = t._rmul_(three_halves) - (half_a * t_cube).truncate_neg(y_prec)# t = (3/2) t - (1/2) a t^3
# print "a =", a
# print "t =", t
# prof("verify")
# print "a*t^2 =", a * t**2
prof("compose")
F_dx_y = (p * x_to_p_less_1 * t) >> (p-1) # px^{p-1} sqrt(a) * y^{-p+1}
# print "-----", F_dx_y
# print "-----", x_to_p * F_dx_y
prof("done")
# print prof
return MonskyWashnitzerDifferential(self, F_dx_y)
def frob_basis_elements(self, prec, p):
"""
Returns the action of a `p`-power lift of Frobenius on the basis
.. math::
\{ dx/2y, x dx/2y, ..., x^{d-2} dx/2y \}
where `d` is the degree of the underlying hyperelliptic curve.
EXAMPLES::
sage: R.<x> = QQ['x']
sage: C = HyperellipticCurve(x^5-4*x+4)
sage: prec = 1
sage: p = 5
sage: MW = C.invariant_differential().parent()
sage: MW.frob_basis_elements(prec,p)
[((92000*y^-14-74200*y^-12+32000*y^-10-8000*y^-8+1000*y^-6-50*y^-4)*1 + (-194400*y^-14+153600*y^-12-57600*y^-10+9600*y^-8-600*y^-6)*x + (204800*y^-14-153600*y^-12+38400*y^-10-3200*y^-8)*x^2 + (-153600*y^-14+76800*y^-12-9600*y^-10)*x^3 + (63950*y^-14-18550*y^-12+1600*y^-10-400*y^-8+50*y^-6+5*y^-4)*x^4) dx/2y, ((-1391200*y^-14+941400*y^-12-302000*y^-10+76800*y^-8-14400*y^-6+1320*y^-4-30*y^-2)*1 + (2168800*y^-14-1402400*y^-12+537600*y^-10-134400*y^-8+16800*y^-6-720*y^-4)*x + (-1596800*y^-14+1433600*y^-12-537600*y^-10+89600*y^-8-5600*y^-6)*x^2 + (1433600*y^-14-1075200*y^-12+268800*y^-10-22400*y^-8)*x^3 + (-870200*y^-14+445350*y^-12-63350*y^-10+3200*y^-8-600*y^-6+30*y^-4+5*y^-2)*x^4) dx/2y, ((19488000*y^-14-15763200*y^-12+4944400*y^-10-913800*y^-8+156800*y^-6-22560*y^-4+1480*y^-2-10)*1 + (-28163200*y^-14+18669600*y^-12-5774400*y^-10+1433600*y^-8-268800*y^-6+25440*y^-4-760*y^-2)*x + (15062400*y^-14-12940800*y^-12+5734400*y^-10-1433600*y^-8+179200*y^-6-8480*y^-4)*x^2 + (-12121600*y^-14+11468800*y^-12-4300800*y^-10+716800*y^-8-44800*y^-6)*x^3 + (9215200*y^-14-6952400*y^-12+1773950*y^-10-165750*y^-8+5600*y^-6-720*y^-4+10*y^-2+5)*x^4) dx/2y, ((-225395200*y^-14+230640000*y^-12-91733600*y^-10+18347400*y^-8-2293600*y^-6+280960*y^-4-31520*y^-2+1480+10*y^2)*1 + (338048000*y^-14-277132800*y^-12+89928000*y^-10-17816000*y^-8+3225600*y^-6-472320*y^-4+34560*y^-2-720)*x + (-172902400*y^-14+141504000*y^-12-58976000*y^-10+17203200*y^-8-3225600*y^-6+314880*y^-4-11520*y^-2)*x^2 + (108736000*y^-14-109760000*y^-12+51609600*y^-10-12902400*y^-8+1612800*y^-6-78720*y^-4)*x^3 + (-85347200*y^-14+82900000*y^-12-31251400*y^-10+5304150*y^-8-367350*y^-6+8480*y^-4-760*y^-2-10+5*y^2)*x^4) dx/2y]
"""
F_i = self.frob_invariant_differential(prec, p)
x_to_p = self.x_to_p(p)
F = [F_i]
for i in range(1, self.degree()-1):
F_i *= x_to_p
F.append(F_i)
return F
def helper_matrix(self):
"""
We use this to solve for the linear combination of
`x^i y^j` needed to clear all terms with
`y^{j-1}`.
"""
try:
return self._helper_matrix
except:
AttributeError
# The smallest y term of (1/j) d(x^i y^j) is constant for all j.
L = []
x, y = self.base_ring().gens()
n = self.degree()
for i in range(n):
L.append( (y*x**i).diff().extract_pow_y(0) )
A = matrix(L).transpose()
if not is_IntegralDomain(A.base_ring()):
# must be using integer_mod or something to approximate
self._helper_matrix = (~A.change_ring(QQ)).change_ring(A.base_ring())
else:
self._helper_matrix = ~A
return self._helper_matrix
class MonskyWashnitzerDifferential(ModuleElement):
def __init__(self, parent, val=0, offset=0):
r"""
Create an element of the Monsky-Washnitzer ring of differentials, of
the form `F dx/2y`.
INPUT:
- ``parent`` - Monsky-Washnitzer differential ring (instance of class
:class:`~MonskyWashnitzerDifferentialRing_class`
- ``val`` - element of the base ring, or list of coefficients
- ``offset`` - if non-zero, shift val by `y^\text{offset}` (default 0)
EXAMPLES::
sage: R.<x> = QQ['x']
sage: C = HyperellipticCurve(x^5 - 4*x + 4)
sage: x,y = C.monsky_washnitzer_gens()
sage: MW = C.invariant_differential().parent()
sage: sage.schemes.elliptic_curves.monsky_washnitzer.MonskyWashnitzerDifferential(MW, x)
x dx/2y
sage: sage.schemes.elliptic_curves.monsky_washnitzer.MonskyWashnitzerDifferential(MW, y)
y*1 dx/2y
sage: sage.schemes.elliptic_curves.monsky_washnitzer.MonskyWashnitzerDifferential(MW, x, 10)
y^10*x dx/2y
"""
ModuleElement.__init__(self, parent)
if isinstance(val, MonskyWashnitzerDifferential):
val = val._coeff
self._coeff = self.parent().base_ring()(val, offset)
def _add_(left, right):
"""
Returns the sum of left and right, both elements of the
Monsky-Washnitzer ring of differentials.
EXAMPLES::
sage: R.<x> = QQ['x']
sage: C = HyperellipticCurve(x^5-4*x + 4)
sage: x,y = C.monsky_washnitzer_gens()
sage: w = C.invariant_differential()
sage: w + w
2*1 dx/2y
sage: x*w + w
(1 + x) dx/2y
sage: x*w + y*w
(y*1 + x) dx/2y
"""
return MonskyWashnitzerDifferential(left.parent(),
left._coeff + right._coeff)
def _sub_(left, right):
"""
Returns the difference of left and right, both elements of the
Monsky-Washnitzer ring of differentials.
EXAMPLES::
sage: R.<x> = QQ['x']
sage: C = HyperellipticCurve(x^5-4*x+4)
sage: x,y = C.monsky_washnitzer_gens()
sage: w = C.invariant_differential()
sage: w-w
0 dx/2y
sage: x*w-w
((-1)*1 + x) dx/2y
sage: w - x*w - y*w
((1-y)*1 + (-1)*x) dx/2y
"""
return MonskyWashnitzerDifferential(left.parent(),
left._coeff - right._coeff)
def __neg__(self):
"""
Returns the additive inverse of self.
EXAMPLES::
sage: R.<x> = QQ['x']
sage: C = HyperellipticCurve(x^5-4*x+4)
sage: x,y = C.monsky_washnitzer_gens()
sage: w = C.invariant_differential()
sage: -w
((-1)*1) dx/2y
sage: -((y-x)*w)
((-y)*1 + x) dx/2y
"""
return MonskyWashnitzerDifferential(self.parent(), -self._coeff)
def _lmul_(self, a):
"""
Returns `self * a`.
EXAMPLES::
sage: R.<x> = QQ['x']
sage: C = HyperellipticCurve(x^5-4*x+4)
sage: x,y = C.monsky_washnitzer_gens()
sage: w = C.invariant_differential()
sage: w*x
x dx/2y
sage: (w*x)*2
2*x dx/2y
sage: w*y
y*1 dx/2y
sage: w*(x+y)
(y*1 + x) dx/2y
"""
return MonskyWashnitzerDifferential(self.parent(), self._coeff * a)
def _rmul_(self, a):
"""
Returns `a * self`.
EXAMPLES::
sage: R.<x> = QQ['x']
sage: C = HyperellipticCurve(x^5-4*x+4)
sage: x,y = C.monsky_washnitzer_gens()
sage: w = C.invariant_differential()
sage: x*w
x dx/2y
sage: 2*(x*w)
2*x dx/2y
sage: y*w
y*1 dx/2y
sage: (x+y)*w
(y*1 + x) dx/2y
"""
return MonskyWashnitzerDifferential(self.parent(), a * self._coeff)
def coeff(self):
r"""
Returns `A`, where this element is `A dx/2y`.
EXAMPLES::
sage: R.<x> = QQ['x']
sage: C = HyperellipticCurve(x^5-4*x+4)
sage: x,y = C.monsky_washnitzer_gens()
sage: w = C.invariant_differential()
sage: w
1 dx/2y
sage: w.coeff()
1
sage: (x*y*w).coeff()
y*x
"""
return self._coeff
def __nonzero__(self):
"""
EXAMPLES::
sage: R.<x> = QQ['x']
sage: C = HyperellipticCurve(x^5-4*x+4)
sage: x,y = C.monsky_washnitzer_gens()
sage: w = C.invariant_differential()
sage: not w
False
sage: not 0*w
True
sage: not x*y*w
False
"""
return not not self._coeff
def _repr_(self):
"""
EXAMPLES::
sage: R.<x> = QQ['x']
sage: C = HyperellipticCurve(x^5-4*x+4)
sage: x,y = C.monsky_washnitzer_gens()
sage: w = C.invariant_differential()
sage: w
1 dx/2y
sage: (2*x+y)*w
(y*1 + 2*x) dx/2y
"""
s = self._coeff._repr_()
if s.find("+") != -1 or s.find("-") > 0:
s = "(%s)"%s
return s + " dx/2y"
def _latex_(self):
"""
Returns the latex representation of self.
EXAMPLES::
sage: R.<x> = QQ['x']
sage: C = HyperellipticCurve(x^5-4*x+4)
sage: x,y = C.monsky_washnitzer_gens()
sage: w = C.invariant_differential()
sage: latex(w)
1 \frac{dx}{2y}
sage: latex(x*w)
x \frac{dx}{2y}
"""
s = self._coeff._latex_()
if s.find("+") != -1 or s.find("-") > 0:
s = "\\left(%s\\right)"%s
return s + " \\frac{dx}{2y}"
def extract_pow_y(self, k):
"""
Returns the power of `y` in `A` where self is `A dx/2y`.
EXAMPLES::
sage: R.<x> = QQ['x']
sage: C = HyperellipticCurve(x^5-3*x+1)
sage: x,y = C.monsky_washnitzer_gens()
sage: A = y^5 - x*y^3
sage: A.extract_pow_y(5)
[1, 0, 0, 0, 0]
sage: (A * C.invariant_differential()).extract_pow_y(5)
[1, 0, 0, 0, 0]
"""
return self._coeff.extract_pow_y(k)
def min_pow_y(self):
"""
Returns the minimum power of `y` in `A` where self is `A dx/2y`.
EXAMPLES::
sage: R.<x> = QQ['x']
sage: C = HyperellipticCurve(x^5-3*x+1)
sage: x,y = C.monsky_washnitzer_gens()
sage: w = y^5 * C.invariant_differential()
sage: w.min_pow_y()
5
sage: w = (x^2*y^4 + y^5) * C.invariant_differential()
sage: w.min_pow_y()
4
"""
return self._coeff.min_pow_y()
def max_pow_y(self):
"""
Returns the maximum power of `y` in `A` where self is `A dx/2y`.
EXAMPLES::
sage: R.<x> = QQ['x']
sage: C = HyperellipticCurve(x^5-3*x+1)
sage: x,y = C.monsky_washnitzer_gens()
sage: w = y^5 * C.invariant_differential()
sage: w.max_pow_y()
5
sage: w = (x^2*y^4 + y^5) * C.invariant_differential()
sage: w.max_pow_y()
5
"""
return self._coeff.max_pow_y()
def reduce_neg_y(self):
"""
Use homology relations to eliminate negative powers of `y`.
EXAMPLES::
sage: R.<x> = QQ['x']
sage: C = HyperellipticCurve(x^5-3*x+1)
sage: x,y = C.monsky_washnitzer_gens()
sage: (y^-1).diff().reduce_neg_y()
((y^-1)*1, 0 dx/2y)
sage: (y^-5*x^2+y^-1*x).diff().reduce_neg_y()
((y^-1)*x + (y^-5)*x^2, 0 dx/2y)
"""
S = self.parent().base_ring()
R = S.base_ring()
M = self.parent().helper_matrix()
p = S._p
n = S.degree()
x, y = S.gens()
f = S(0)
reduced = self
for j in range(self.min_pow_y()+1, 0):
if p is not None and p.divides(j):
cs = [a/j for a in reduced.extract_pow_y(j-1)]
else:
j_inverse = ~R(j)
cs = [a*j_inverse for a in reduced.extract_pow_y(j-1)]
lin_comb = M * vector(M.base_ring(), cs)
# print "j =", j, "b =", cs, "lin_comb =", lin_comb
g = self.parent().base_ring()(0)
if not lin_comb.is_zero():
for i in range(n):
if lin_comb[i] != 0:
g += S.monomial(i, j, lin_comb[i])
if not g.is_zero():
f += g
reduced -= g.diff()
# print g, g.diff()
# print "reduced", reduced
return f, reduced
def reduce_neg_y_fast(self, even_degree_only=False):
"""
Use homology relations to eliminate negative powers of `y`.
EXAMPLES::
sage: R.<x> = QQ['x']
sage: E = HyperellipticCurve(x^5-3*x+1)
sage: x, y = E.monsky_washnitzer_gens()
sage: (y^-1).diff().reduce_neg_y_fast()
((y^-1)*1, 0 dx/2y)
sage: (y^-5*x^2+y^-1*x).diff().reduce_neg_y_fast()
((y^-1)*x + (y^-5)*x^2, 0 dx/2y)
It leaves non-negative powers of `y` alone::
sage: y.diff()
((-3)*1 + 5*x^4) dx/2y
sage: y.diff().reduce_neg_y_fast()
(0, ((-3)*1 + 5*x^4) dx/2y)
"""
# prof = Profiler()
# prof("reduce setup")
S = self.parent().base_ring()
R = S.base_ring()
M = self.parent().helper_matrix()
# prof("extract coeffs")
coeffs, offset = self.coeffs(R)
V = coeffs[0].parent()
if offset == 0:
return S(0), self
# prof("loop %s"%self.min_pow_y())
forms = []
p = S._p
for j in range(self.min_pow_y()+1, 0):
if (even_degree_only and j % 2 == 0) or coeffs[j-offset-1].is_zero():
forms.append(V(0))
else:
# this is a total hack to deal with the fact that we're using
# rational numbers to approximate fixed precision p-adics
if p is not None and j % 3 == 1:
try:
v = coeffs[j-offset-1]
for kk in range(len(v)):
a = v[kk]
ppow = p**max(-a.valuation(S._p), 0)
v[kk] = ((a * ppow) % S._prec_cap) / ppow
except AttributeError:
pass
lin_comb = ~R(j) * (M * coeffs[j-offset-1])
forms.append(lin_comb)
for i in lin_comb.nonzero_positions():
# g = lin_comb[i] x^i y^j
# self -= dg
coeffs[j-offset+1] -= lin_comb[i] * S.monomial_diff_coeffs(i, j)[1]
# prof("recreate forms")
f = S(forms, offset+1)
reduced = S._monsky_washnitzer(coeffs[-1-offset:], -1)
# print self - f.diff() - reduced
# prof("done")
# print prof
return f, reduced
def reduce_neg_y_faster(self, even_degree_only=False):
"""
Use homology relations to eliminate negative powers of `y`.
"""
# Timings indicate this isn't any faster after all...
S = self.parent().base_ring()
R = S.base_ring()
M = self.parent().helper_matrix()
coeffs, offset = self.coeffs(R)
V = coeffs[0].parent()
zeroV = V(0)
if offset == 0:
return S(0), self
# See monomial_diff_coeffs
# this is the B_i and x_to_i contributions respectively for all i
d_mat_1, d_mat_2 = S.monomial_diff_coeffs_matrices()
forms = []
for j in range(self.min_pow_y()+1, 0):
if coeffs[j-offset-1].is_zero():
forms.append(zeroV)
else:
# this is a total hack to deal with the fact that we're using
# rational numbers to approximate fixed precision p-adics
if j % 3 == 0:
try:
v = coeffs[j-offset-1]
for kk in range(len(v)):
a = v[kk]
ppow = S._p**max(-a.valuation(S._p), 0)
v[kk] = ((a * ppow) % S._prec_cap) / ppow
except AttributeError:
pass
j_inverse = ~R(j)
lin_comb = (M * coeffs[j-offset-1])
forms.append(j_inverse * lin_comb)
coeffs[j-offset+1] -= (d_mat_1 + j_inverse * d_mat_2) * lin_comb
f = S(forms, offset+1)
reduced = S._monsky_washnitzer(coeffs[-1-offset:], -1)
# reduced = self - f.diff()
return f, reduced
def reduce_pos_y(self):
"""
Use homology relations to eliminate positive powers of `y`.
EXAMPLES::
sage: R.<x> = QQ['x']
sage: C = HyperellipticCurve(x^3-4*x+4)
sage: x,y = C.monsky_washnitzer_gens()
sage: (y^2).diff().reduce_pos_y()
(y^2*1, 0 dx/2y)
sage: (y^2*x).diff().reduce_pos_y()
(y^2*x, 0 dx/2y)
sage: (y^92*x).diff().reduce_pos_y()
(y^92*x, 0 dx/2y)
sage: w = (y^3 + x).diff()
sage: w += w.parent()(x)
sage: w.reduce_pos_y_fast()
(y^3*1 + x, x dx/2y)
"""
S = self.parent().base_ring()
series = S.base_ring()
n = S.Q().degree()
p = S._p
x, y = S.gens()
f = S(0)
reduced = self
for j in range(self.max_pow_y(), 0, -1):
for i in range(n-1, -1, -1):
c = reduced.extract_pow_y(j)[i]
# print "x^%s y^%s"%(i,j), c
if c != 0:
g = S.monomial(0, j+1) if i == n-1 else S.monomial(i+1, j-1)
dg = g.diff()
# print reduced, " - ", dg
denom = dg.extract_pow_y(j)[i]
c /= denom
c = g.parent()(c)
f += c * g
reduced -= c * dg
return f, reduced
def reduce_pos_y_fast(self, even_degree_only=False):
"""
Use homology relations to eliminate positive powers of `y`.
EXAMPLES::
sage: R.<x> = QQ['x']
sage: E = HyperellipticCurve(x^3-4*x+4)
sage: x, y = E.monsky_washnitzer_gens()
sage: y.diff().reduce_pos_y_fast()
(y*1, 0 dx/2y)
sage: (y^2).diff().reduce_pos_y_fast()
(y^2*1, 0 dx/2y)
sage: (y^2*x).diff().reduce_pos_y_fast()
(y^2*x, 0 dx/2y)
sage: (y^92*x).diff().reduce_pos_y_fast()
(y^92*x, 0 dx/2y)
sage: w = (y^3 + x).diff()
sage: w += w.parent()(x)
sage: w.reduce_pos_y_fast()
(y^3*1 + x, x dx/2y)
"""
S = self.parent().base_ring()
R = S.base_ring()
n = S.Q().degree()
coeffs, offset = self.coeffs(R)
V = coeffs[0].parent()
zeroV = V(0)
forms = [V(0), V(0)]
for j in range(self.max_pow_y(), -1, -1):
if (even_degree_only and j % 2 == 1) or (j > 0 and coeffs[j-offset].is_zero()):
forms.append(zeroV)
continue
form = V(0)
i = n-1
c = coeffs[j-offset][i]
if c != 0:
dg_coeffs = S.monomial_diff_coeffs(0, j+1)[0]
c /= dg_coeffs[i]
forms[len(forms)-2][0] = c
# self -= c d(y^{j+1})
coeffs[j-offset] -= c*dg_coeffs
if j == 0:
# the others are basis elements
break
for i in range(n-2, -1, -1):
c = coeffs[j-offset][i]
if c != 0:
dg_coeffs = S.monomial_diff_coeffs(i+1, j-1)
denom = dg_coeffs[1][i]
c /= denom
form[i+1] = c
# self -= c d(x^{i+1} y^{j-1})
coeffs[j-offset] -= c*dg_coeffs[1]
coeffs[j-offset-2] -= c*dg_coeffs[0]
forms.append(form)
forms.reverse()
f = S(forms)
reduced = self.parent()(coeffs[:1-offset], offset)
return f, reduced
def reduce(self):
"""
Use homology relations to find `a` and `f` such that this element is
equal to `a + df`, where `a` is given in terms of the `x^i dx/2y`.
EXAMPLES::
sage: R.<x> = QQ['x']
sage: C = HyperellipticCurve(x^5-4*x+4)
sage: x,y = C.monsky_washnitzer_gens()
sage: w = (y*x).diff()
sage: w.reduce()
(y*x, 0 dx/2y)
sage: w = x^4 * C.invariant_differential()
sage: w.reduce()
(1/5*y*1, 4/5*1 dx/2y)
sage: w = sum(QQ.random_element() * x^i * y^j for i in [0..4] for j in [-3..3]) * C.invariant_differential()
sage: f, a = w.reduce()
sage: f.diff() + a - w
0 dx/2y
"""
# print "max_pow_y = ", self.max_pow_y(), "min_pow_y = ", self.min_pow_y()
n = self.parent().base_ring().Q().degree()
f1, a = self.reduce_neg_y()
f2, a = a.reduce_pos_y()
f = f1 + f2
c = a.extract_pow_y(0)[n-1]
if c != 0:
x, y = self.parent().base_ring().gens()
g = y
dg = g.diff()
c = g.parent()(c/dg.extract_pow_y(0)[n-1])
f += c * g
a -= c * dg
# print g, dg
return f, a
def reduce_fast(self, even_degree_only=False):
"""
Use homology relations to find `a` and `f` such that this element is
equal to `a + df`, where `a` is given in terms of the `x^i dx/2y`.
EXAMPLES::
sage: R.<x> = QQ['x']
sage: E = HyperellipticCurve(x^3-4*x+4)
sage: x, y = E.monsky_washnitzer_gens()
sage: x.diff().reduce_fast()
(x, (0, 0))
sage: y.diff().reduce_fast()
(y*1, (0, 0))
sage: (y^-1).diff().reduce_fast()
((y^-1)*1, (0, 0))
sage: (y^-11).diff().reduce_fast()
((y^-11)*1, (0, 0))
sage: (x*y^2).diff().reduce_fast()
(y^2*x, (0, 0))
"""
# print "max_pow_y = ", self.max_pow_y(), "min_pow_y = ", self.min_pow_y()
f1, reduced = self.reduce_neg_y_fast(even_degree_only)
f2, reduced = reduced.reduce_pos_y_fast(even_degree_only)
# f1, reduced = self.reduce_neg_y()
# f2, reduced = reduced.reduce_pos_y()
v = reduced.extract_pow_y(0)
v.pop()
V = FreeModule(self.base_ring().base_ring(), len(v))
return f1+f2, V(v)
def coeffs(self, R=None):
"""
Used to obtain the raw coefficients of a differential, see
:meth:`SpecialHyperellipticQuotientElement.coeffs`
INPUT:
- R - An (optional) base ring in which to cast the coefficients
OUTPUT:
The raw coefficients of $A$ where self is $A dx/2y$.
EXAMPLES::
sage: R.<x> = QQ['x']
sage: C = HyperellipticCurve(x^5-4*x+4)
sage: x,y = C.monsky_washnitzer_gens()
sage: w = C.invariant_differential()
sage: w.coeffs()
([(1, 0, 0, 0, 0)], 0)
sage: (x*w).coeffs()
([(0, 1, 0, 0, 0)], 0)
sage: (y*w).coeffs()
([(0, 0, 0, 0, 0), (1, 0, 0, 0, 0)], 0)
sage: (y^-2*w).coeffs()
([(1, 0, 0, 0, 0), (0, 0, 0, 0, 0), (0, 0, 0, 0, 0)], -2)
"""
return self._coeff.coeffs(R)
def coleman_integral(self, P, Q):
r"""
Computes the definite integral of self from $P$ to $Q$.
INPUT:
- P, Q - Two points on the underlying curve
OUTPUT:
`\int_P^Q \text{self}`
EXAMPLES::
sage: K = pAdicField(5,7)
sage: E = EllipticCurve(K,[-31/3,-2501/108]) #11a
sage: P = E(K(14/3), K(11/2))
sage: w = E.invariant_differential()
sage: w.coleman_integral(P,2*P)
O(5^6)
sage: Q = E.lift_x(3)
sage: w.coleman_integral(P,Q)
2*5 + 4*5^2 + 3*5^3 + 4*5^4 + 3*5^5 + O(5^6)
sage: w.coleman_integral(2*P,Q)
2*5 + 4*5^2 + 3*5^3 + 4*5^4 + 3*5^5 + O(5^6)
sage: (2*w).coleman_integral(P, Q) == 2*(w.coleman_integral(P, Q))
True
"""
return self.parent().base_ring().curve().coleman_integral(self, P, Q)
integrate = coleman_integral
### end of file