working Coleman integrals

src/sage/rings/multi_polynomial_element.py

@@ -237,6 +237,9 @@ def __rpow__(self, n):
     def element(self):
         return self.__element
 
+    def change_ring(self, R):
+        return self.parent().change_ring(R)(self)
+
 
 class MPolynomial_macaulay2_repr:
     """

src/sage/rings/multi_polynomial_ring.py

@@ -329,6 +329,10 @@ def latex_variable_names(self):
         self.__latex_variable_names = names
         return names
 
+    def change_ring(self, R):
+        from polynomial_ring_constructor import PolynomialRing
+        return PolynomialRing(R, self.variable_names(), order=self.term_order())
+
 class MPolynomialRing_polydict( MPolynomialRing_macaulay2_repr, MPolynomialRing_generic):
     """
     Multivariable polynomial ring.

src/sage/rings/power_series_ring_element.pyx

@@ -837,6 +837,8 @@ cdef class PowerSeries(AlgebraElement):
             t^5 + O(t^25)
         """
         denom = <PowerSeries>denom_r
+        if denom.is_zero():
+            raise ZeroDivisionError, "Can't divide by something indistinguishable from 0"
         u = denom.valuation_zero_part()
         inv = ~u  # inverse
 
@@ -980,6 +982,8 @@ cdef class PowerSeries(AlgebraElement):
             1 + t + t^2 + 7*t^3 + O(t^50)
             sage: sqrt(K(0))
             0
+            sage: sqrt(t^2)
+            t
 
         AUTHOR:
             -- Robert Bradshaw
@@ -991,8 +995,12 @@ cdef class PowerSeries(AlgebraElement):
         if val is not infinity and val % 2 == 1:
             raise ValueError, "Square root not defined for power series of odd valuation."
 
+        if self.degree() == 0:
+            x = self.valuation_zero_part()[0].sqrt()
+            return self.parent().base_extend(x.parent())([x], val/2)
+
         prec = self.prec()
-        if prec is infinity and self.degree() > 0:
+        if prec == infinity:
             prec = self._parent.default_prec()
 
         a = self.valuation_zero_part()
@@ -1045,6 +1053,8 @@ cdef class PowerSeries(AlgebraElement):
             raise ValueError, "Square root not defined for power series of odd valuation."
         elif not self[val].is_square():
             raise ValueError, "Square root does not live in this ring."
+        elif is_Field(self.base_ring()):
+            return self.sqrt()
         else:
             try:
                 return self.parent()(self.sqrt())

src/sage/rings/rational.pyx

@@ -163,6 +163,7 @@ cdef class Rational(sage.structure.element.FieldElement):
 
     def __set_value(self, x, unsigned int base):
         cdef int n
+        cdef Rational temp_rational
 
         if isinstance(x, Rational):
             set_from_Rational(self, x)
@@ -236,6 +237,10 @@ cdef class Rational(sage.structure.element.FieldElement):
             if n or mpz_cmp_si(mpq_denref(self.value), 0) == 0:
                 raise TypeError, "Unable to coerce %s (%s) to Rational"%(x,type(x))
 
+        elif hasattr(x, 'rational_reconstruction'):
+            temp_rational = x.rational_reconstruction()
+            mpq_set(self.value, temp_rational.value)
+
         else:
 
             raise TypeError, "Unable to coerce %s (%s) to Rational"%(x,type(x))

src/sage/schemes/elliptic_curves/ell_padic_field.py

@@ -25,10 +25,15 @@
 from sage.rings.all import PowerSeriesRing, PolynomialRing, IntegerModRing, ZZ, QQ
 from sage.misc.functional import ceil, log
 
+# Elliptic curves are very different than genus > 1 hyperelliptic curves,
+# there is an "is a" relationship here, and common implementation with regard
+# Coleman integration.
+from sage.schemes.hyperelliptic_curves.hyperelliptic_padic_field import HyperellipticCurve_padic_field
+
 import sage.databases.cremona
 
 
-class EllipticCurve_padic_field(EllipticCurve_field):
+class EllipticCurve_padic_field(EllipticCurve_field, HyperellipticCurve_padic_field):
     """
     Elliptic curve over a padic field.
     """
@@ -40,7 +45,7 @@ def __init__(self, x, y=None):
         else:
             if isinstance(y, str):
                 field = x
-                X = sage.databases.cremona.CremonaDatabase()[label]
+                X = sage.databases.cremona.CremonaDatabase()[y]
                 ainvs = [field(a) for a in X.a_invariants()]
             else:
                 field = x
@@ -51,6 +56,7 @@ def __init__(self, x, y=None):
         EllipticCurve_field.__init__(self, [field(x) for x in ainvs])
 
         self._point_class = ell_point.EllipticCurvePoint_field
+        self._genus = 1
 
     def _pari_(self):
         try:
@@ -62,63 +68,6 @@ def _pari_(self):
         return self.__pari
 
 
-# The functions below were prototyped at the 2007 Arizona Winter School by
-# Robert Bradshaw and Ralf Gerkmann, working with Miljan Brakovevic and
-# Kiran Kedlaya
-# All of the below is with respect to the Monsky Washnitzer cohomology.
-
-    def local_analytic_interpolation(self, P, Q):
-        """
-        Construct a linear interpolation from P to Q in a power series
-        in the local parameter t, with precision equal to the p-adic
-        precision of the underlying ring.
-
-        Returns a point $X(t) = ( x(t) : y(t) : z(t) )$ such that
-        $X(0) = P and X(1) = Q$
-
-        This is implemented by linearly interpolating x, solving formally for y,
-        and letting z(t) = 1.
-
-        For this to make sense, P and Q must be in the same residue series, neither equal to infinity.
-
-        TODO: remove the last condition?
-        """
-        prec = self.base_ring().precision_cap()
-        t = PowerSeriesRing(self.base_ring(), 't', prec).gen(0)
-        x = P[0]+t*(Q[0]-P[0])
-        pts = self.lift_x(x)
-        if (pts[0][1] - P[1]).valuation() > 0:
-            return pts[0]
-        else:
-            return pts[1]
-
-    def tiny_integrals(self, F, P, Q):
-        """
-        Evaluate the integrals of $f_i dx/y$ from P to Q for each f_i in F
-        by formally integrating a power series in a local parameter $t$
-
-        P and Q MUST be in the same residue disk for this result to make sense.
-        """
-        x, y, z = self.local_analytic_interpolation(P, Q)
-        dt = x.derivative() / y
-        integrals = []
-        for f in F:
-            f = f(x,y)
-            I = (f*dt).integral()
-            integrals.append(I(1))
-        return integrals
-
-    def tiny_integrals_on_basis(self, P, Q):
-        """
-        Evaluate the integrals of $dx/y$ and $x dx/y$
-        by formally integrating a power series in a local parameter $t$
-
-        P and Q MUST be in the same residue disk for this result to make sense.
-        """
-        R = PolynomialRing(self.base_ring(), ['x', 'y'])
-        return self.tiny_integrals([R(1), R.gen(0)], P, Q)
-
-
     def frobenius(self, P=None):
 
         try:
@@ -152,35 +101,6 @@ def frob(P):
         else:
             return frob(P)
 
-    def teichmuller(self, P):
-        """
-        Find a Teichm\:uller point in the same residue class of P.
-
-        Because this lift of frobenius acts as $x \mapsto x^p$,
-        take the Teichmuler lift of $x$ and then find a matching y
-        from that.
-
-        EXAMPLES:
-            sage: K = pAdicField(7, 5)
-            sage: E = EllipticCurve(K, [-31/3, -2501/108]) # 11a
-            sage: P = E(K(14/3), K(11/2))
-            sage: E.frobenius(P) == P
-            False
-            sage: TP = E.teichmuller(P); TP
-            (0 : 2 + 3*7 + 3*7^2 + 3*7^4 + O(7^5) : 1 + O(7^5))
-            sage: E.frobenius(TP) == TP
-            True
-            sage: (TP[0] - P[0]).valuation() > 0, (TP[1] - P[1]).valuation() > 0
-            (True, True)
-        """
-        x = padic_teichmuller(P[0])
-        pts = self.lift_x(x)
-        if (pts[0][1] - P[1]).valuation() > 0:
-            return pts[0]
-        else:
-            return pts[1]
-
-
     def coleman_integrals_on_basis(self, P, Q):
         """
         Return the coleman integral of dx/y and x dx/y from P to Q.
@@ -254,16 +174,3 @@ def coleman_integrals_on_basis(self, P, Q):
         return P_to_TP + TP_to_TQ + TQ_to_Q
 
 
-
-# TODO: add this to padics (if it isn't there in the new version already).
-def padic_teichmuller(a):
-    K = a.parent()
-    p = K.prime()
-    p_less_1_inverse = ~K(p-1)
-    one = K(1)
-    x = K(ZZ(a) % p)
-    if x == 0:
-        return x
-    for _ in range(ceil(log(K.precision_cap(),2))):
-        x = ((p-2)*x + x**(2-p)) * p_less_1_inverse
-    return x

src/sage/schemes/elliptic_curves/ell_rational_field.py

@@ -4404,8 +4404,7 @@ def padic_sigma(self, p, N=20, E2=None, check=False, check_hypotheses=True):
 
 
 
-    def padic_E2(self, p, prec=20, check=False, check_hypotheses=True,
-                 algorithm="auto"):
+    def padic_E2(self, p, prec=20, check=False, check_hypotheses=True, algorithm="auto"):
         r"""
         Returns the value of the $p$-adic modular form $E2$ for $(E, \omega)$
         where $\omega$ is the usual invariant differential
@@ -4531,6 +4530,38 @@ def padic_E2(self, p, prec=20, check=False, check_hypotheses=True,
             1907 + 2819*3001 + 1124*3001^2 + O(3001^3)
 
         """
+        frob_p = self.matrix_of_frobenius(p, prec, check, check_hypotheses, algorithm).change_ring(Integers(p**prec))
+
+        frob_p_n = frob_p**prec
+
+        # todo: write a coercion operator from integers mod p^n to the
+        # p-adic field (it doesn't seem to currently exist)
+        # see trac #4
+
+        # todo: think about the sign of this. Is it correct?
+
+        E2_of_X = output_ring( (-12 * frob_p_n[0,1] / frob_p_n[1,1]).lift() ) \
+                  + O(p**prec)
+
+        # Take into account the coordinate change.
+        fudge_factor = (X.discriminant() / self.discriminant()).nth_root(6)
+
+        # todo: here I should be able to write:
+        #  return E2_of_X / fudge_factor
+        # However, there is a bug in SAGE (#51 on trac) which makes this
+        # crash sometimes when prec == 1. For example,
+        #    EllipticCurve([1, 1, 1, 1, 1]).padic_E2(5, 1)
+        # makes it crash. I haven't figured out exactly what the bug
+        # is yet, but for now I use the following workaround:
+        fudge_factor_inverse = Qp(p, prec=(E2_of_X.precision_absolute() + 1))(1 / fudge_factor)
+        return output_ring(E2_of_X * fudge_factor_inverse)
+
+
+    def matrix_of_frobenius(self, p, prec=20, check=False, check_hypotheses=True, algorithm="auto"):
+        """
+        See the parameters and documentation for padic_E2.
+        """
+        # TODO, add lots of comments like the above
         if check_hypotheses:
             p = self.__check_padic_hypotheses(p)
 
@@ -4617,29 +4648,7 @@ def padic_E2(self, p, prec=20, check=False, check_hypotheses=True,
                     "Consistency check failed! (correct = %s, actual = %s)" % \
                     (correct_trace, trace_of_frobenius)
 
-        frob_p_n = frob_p**prec
-
-        # todo: write a coercion operator from integers mod p^n to the
-        # p-adic field (it doesn't seem to currently exist)
-        # see trac #4
-
-        # todo: think about the sign of this. Is it correct?
-
-        E2_of_X = output_ring( (-12 * frob_p_n[0,1] / frob_p_n[1,1]).lift() ) \
-                  + O(p**prec)
-
-        # Take into account the coordinate change.
-        fudge_factor = (X.discriminant() / self.discriminant()).nth_root(6)
-
-        # todo: here I should be able to write:
-        #  return E2_of_X / fudge_factor
-        # However, there is a bug in SAGE (#51 on trac) which makes this
-        # crash sometimes when prec == 1. For example,
-        #    EllipticCurve([1, 1, 1, 1, 1]).padic_E2(5, 1)
-        # makes it crash. I haven't figured out exactly what the bug
-        # is yet, but for now I use the following workaround:
-        fudge_factor_inverse = Qp(p, prec=(E2_of_X.precision_absolute() + 1))(1 / fudge_factor)
-        return output_ring(E2_of_X * fudge_factor_inverse)
+        return frob_p.change_ring(Zp(p, prec))
 
 
     # This is the old version of padic_E2 that requires MAGMA:

src/sage/schemes/elliptic_curves/monsky_washnitzer.py

@@ -2021,9 +2021,17 @@ def matrix_of_frobenius_alternate(a, b, p, N):
 from sage.misc.profiler import Profiler
 from sage.misc.misc import repr_lincomb
 
-def matrix_of_frobenius_hyperelliptic(Q, p, prec, M=None):
+def matrix_of_frobenius_hyperelliptic(Q, p=None, prec=None, M=None):
     prof = Profiler()
     prof("setup")
+    if p is None:
+        try:
+            K = Q.base_ring()
+            p = K.prime()
+            prec = K.precision_cap()
+            #prec -= (adjusted_prec(p, prec) - prec)
+        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) # pAdicField(p, M) # SLOW!
@@ -2037,6 +2045,7 @@ def matrix_of_frobenius_hyperelliptic(Q, p, prec, M=None):
     # 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 coefficent size)
     rational_S._prec_cap = p**M
     rational_S._p = p
@@ -2060,7 +2069,7 @@ def matrix_of_frobenius_hyperelliptic(Q, p, prec, M=None):
     M += real_prec_ring(0).add_bigoh(prec)
     print prof
 #    print len(S._monomials)
-    return M.transpose()
+    return M.transpose(), [f for f, a in reduced]
 
 
 
@@ -2091,21 +2100,24 @@ def MonskyWashnitzerDifferentialRing(base_ring):
 
 
 class SpecialHyperellipticQuotientRing_class(CommutativeAlgebra):
-    def __init__(self, Q, R, invert_y=False):
+    def __init__(self, Q, R=None, invert_y=True):
+        if R is None:
+            R = Q.base_ring()
+
         CommutativeAlgebra.__init__(self, R)
 
         x = PolynomialRing(R, 'x').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.defining_polynomial()(x,0,1)
+            Q = (-E.defining_polynomial()).change_ring(R)(x,0,1)
 
         elif is_HyperellipticCurve(Q):
             C = Q
             if C.hyperelliptic_polynomials()[1] != 0:
                 raise NotImplementedError, "Curve must be of form y^2 = Q(x)."
-            Q = E.hyperelliptic_polynomials()[0]()(x)
+            Q = C.hyperelliptic_polynomials()[0].change_ring(R)
 
         if is_Polynomial(Q):
             self._Q = Q.change_ring(R)
@@ -2262,7 +2274,7 @@ def __init__(self, parent, val=0, offset=0):
         if isinstance(val, tuple):
             val, offset = val
         if isinstance(val, SpecialHyperellipticQuotientElement):
-            val = val.coeffs()
+            val, offset = val.coeffs()
         if isinstance(val, list) and len(val) > 0 and is_FreeModuleElement(val[0]):
             val = transpose_list(val)
         self._f = parent._poly_ring(val)
@@ -2272,6 +2284,9 @@ def __init__(self, parent, val=0, offset=0):
     def change_ring(self, R):
         return self.parent().change_ring(R)(self.coeffs())
 
+    def __call__(self, *x):
+        return self._f(*x)
+
     def __invert__(self):
         """
         The general element in our ring is not invertible, but y may be.