Changeset 3190:b1427bf737c5

Timestamp:

02/25/07 21:00:36 (6 years ago)

Author:

dmharvey@…

Branch:

default

Message:

algorithm for computing padic E2 in time sqrt(p)

Location:

sage/schemes/elliptic_curves

Files:

• ell_rational_field.py (modified) (3 diffs)
• monsky_washnitzer.py (modified) (5 diffs)

sage/schemes/elliptic_curves/ell_rational_field.py

@@ -40 +40 @@
 from sage.misc.functional import log
 
-# Use some interval arithmetic to guaranty correctness.  We assume
+# Use some interval arithmetic to guarantee correctness.  We assume
 # that alpha is computed to the precision of a float.
 IR = rings.RIF
@@ -4388 +4388 @@
             2 + 4*5 + 2*5^3 + 5^4 + 3*5^5 + 2*5^6 + 5^8 + 3*5^9 + 4*5^10 + 2*5^11 + 2*5^12 + 2*5^14 + 3*5^15 + 3*5^16 + 3*5^17 + 4*5^18 + 2*5^19 + 4*5^20 + 5^21 + 4*5^22 + 2*5^23 + 3*5^24 + 3*5^26 + 2*5^27 + 3*5^28 + O(5^30)
 
+        Here's one using the $p^{1/2}$ algorithm:
+            sage: EllipticCurve([-1, 1/4]).padic_E2(3001, 1, check=True)  # long time (< 10s)
+            1907 + 2819*3001 + 1124*3001^2 + O(3001^3)
+
         """
         if check_hypotheses:
@@ -4428 +4432 @@
                " at p."
 
-        # Need to increase precision a little to compensate for precision
-        # losses during the computation. (See monsky_washnitzer.py
-        # for more details.)
-        adjusted_prec = monsky_washnitzer.adjusted_prec(p, prec)
-        
-        if check:
-            trace = None
+        # If p is large, we use the matrix_of_frobenius_alternate
+        # function, which has better asymptotic performance for large p.
+        # The right cutoff should depend on N, but this has not been
+        # investigated properly yet, so for now we have an arbitrary
+        # cutoff at p = 3000, which works decently on sage.math when
+        # the precision is very low.
+
+        if p < 3000:
+            # Need to increase precision a little to compensate for precision
+            # losses during the computation. (See monsky_washnitzer.py
+            # for more details.)
+            adjusted_prec = monsky_washnitzer.adjusted_prec(p, prec)
+        
+            if check:
+                trace = None
+            else:
+                trace = self.ap(p)
+
+            base_ring = rings.Integers(p**adjusted_prec)
+            output_ring = rings.pAdicField(p, prec)
+        
+            R, x = rings.PolynomialRing(base_ring, 'x').objgen()
+            Q = x**3 + base_ring(X.a4()) * x + base_ring(X.a6())
+            frob_p = monsky_washnitzer.matrix_of_frobenius(
+                             Q, p, adjusted_prec, trace)
+
         else:
-            trace = self.ap(p)
-
-        base_ring = rings.Integers(p**adjusted_prec)
-        output_ring = rings.pAdicField(p, prec)
-        
-        R, x = rings.PolynomialRing(base_ring, 'x').objgen()
-        Q = x**3 + base_ring(X.a4()) * x + base_ring(X.a6())
-        frob_p = monsky_washnitzer.matrix_of_frobenius(
-                        Q, p, adjusted_prec, trace)
+            base_ring = rings.Integers(p**prec)
+            output_ring = rings.pAdicField(p, prec)
+            frob_p = monsky_washnitzer.matrix_of_frobenius_alternate(
+                          X.a4(), X.a6(), p, prec)
+            
 
         if check:

sage/schemes/elliptic_curves/monsky_washnitzer.py

@@ -3 +3 @@
 
 The most interesting functions to be exported here are matrix_of_frobenius()
-and adjusted_prec().
+and adjusted_prec(), and matrix_of_frobenius_alternate().
 
 Currently this code is limited to the case $p \geq 5$ (no $GF(p^n)$
@@ -15 +15 @@
   -- Edixhoven, B., ``Point counting after Kedlaya'', EIDMA-Stieltjes graduate
      course, Lieden (lecture notes?).
+  -- Harvey, D. ``Kedlaya's algorithm in larger characteristic'',
+     http://arxiv.org/abs/math.NT/0610973
 
 AUTHORS:
@@ -23 +25 @@
        more documentation, added Newton iteration method, added more complete
        "trace trick", integrated better into SAGE.
+    -- David Harvey (Feb 2007): added algorithm with sqrt(p) complexity
 
 """
@@ -40 +43 @@
 from sage.rings.ring import CommutativeAlgebra
 from sage.structure.element import CommutativeAlgebraElement
-
-from sage.rings.arith import binomial
-from sage.misc.functional import log, ceil
+from sage.matrix.matrix_space import MatrixSpace
+
+from sage.rings.arith import binomial, floor
+from sage.misc.functional import log, ceil, sqrt
 
 
@@ -1256 +1260 @@
 
 
+
+#*****************************************************************************
+# From here on is an implementation of D. Harvey's modification of Kedlaya's
+# algorithm, ``Kedlaya's algorithm in larger characteristic'' (see arXiv).
+# Actually the algorithm implemented here is better than the one in the paper;
+# the $N^3$ contribution is reduced to $N^{5/2}$, and incorporates some ideas
+# from ``Linear recurrences with polynomial coefficients and application to
+# integer factorization and Cartier-Manin operator'' (Bostan, Gaudry, Schost)
+# to reduce the running time by another factor of $\log p$.
+#
+# todo: I haven't yet done the precision analysis (i.e. the analogue of
+# lemmas 2 and 3 from Kedlaya's original paper), so for now I'm playing it
+# safe, using more precision than is presumably necessary.
+#
+# AUTHOR:
+#    -- David Harvey (2007-02)
+#
+#*****************************************************************************
+
+
+def shift_evaluation_points(R, values, b, r, p, cache):
+    r"""
+    Given the values of a polynomial $P(x)$ on an arithmetic progression,
+    this function computes the values of $P(x)$ on a shift of that
+    arithmetic progression.
+    
+    INPUT:
+      p -- a prime
+      R -- coefficient ring, of the form $\ZZ/p^m \ZZ$
+      values -- a list of $d+1$ values of a polynomial $P(x)$ of degree $d$,
+          at the $d+1$ points $x = 0, r, 2r, \ldots, dr$.
+      b -- an integer
+      r -- a nonzero integer
+      cache -- a dictionary
+
+    OUTPUT:
+      A list of $d+1$ values $P(b)$, $P(b+r)$, ..., $P(b+dr)$.
+      (Also the cache might get updated.)
+
+    NOTE: The polynomial $P$ is *not* part of the input.
+
+    PRECONDITIONS: $d$ must be less than $p$
+      
+    EXAMPLES:
+    Run it on some random input, and compare to a naive implementation:
+        sage: from sage.schemes.elliptic_curves.monsky_washnitzer import shift_evaluation_points
+        sage: R = Integers(7**3)
+        sage: S.<x> = PolynomialRing(R)
+        sage: r = 3
+        sage: for d in range(7):
+        ...      for b in range(50):
+        ...         p = S([R.random_element() for i in range(d+1)])
+        ...         L = [p(i*r) for i in range(d+1)]
+        ...         shifted = shift_evaluation_points(R, L, b, r, 7, {})
+        ...         if shifted != [p(b + i*r) for i in range(d+1)]:
+        ...            print "failed"
+
+    ALGORITHM:
+    The basic idea (as described in the Bostan/Gaudry/Schost paper) is as
+    follows. The Lagrange interpolation formula says that
+      $$ P(x) = \sum_{i=0}^d P(ir) \prod_{j \neq i} (x-jr)/(ir-jr). $$
+    We substitute $x = b+kr$, and rearrange to eventually obtain
+      $$ P(b+kr) = r^{-d} \prod_{j=0}^d (b+(k-j)r)
+        \sum_{i=0}^d [P(ir) \prod_{j \neq i} (i-j)^{-1}] \frac{1}{b+(k-i)r}. $$
+    The point is that the latter sum is a convolution and so can be
+    evaluated for all $k$ simultaneously by a single polynomial multiplication.
+
+    There are some complications related to non-invertible elements of $R$,
+    but we just deal with it :-)
+        
+    """
+
+    d = len(values) - 1
+    assert p > d
+
+    try:
+        # try to retrieve some precomputed stuff from the cache
+        bad, input_scale, S, kernel, output_scale = cache[d, b, r]
+
+    except KeyError:
+        # that didn't work; need to compute all of it
+        
+        # figure out for which i in [-d, d] is b+ir not invertible;
+        # put these i into a "bad" list
+        i = (-R(b)/r).lift() % p
+        bad = []
+        if i <= d:
+            bad.append(i)
+        if i >= p-d:
+            bad.append(i-p)
+
+        # compute inverse factorials [1/0!, 1/1!, 1/2!, ..., 1/d!] in R
+        total = R(1)
+        for k in range(2, d+1):
+            total = total * R(k)
+        inv = total**(-1)
+        inv_fac = [inv]
+        for k in range(d, 1, -1):
+            inv_fac.append(inv_fac[-1] * R(k))
+        inv_fac.append(R(1))
+        inv_fac.reverse()
+
+        # compute some coefficients to be used in the preprocessing step
+        scale = R(r)**(-d)
+        input_scale = [scale * inv_fac[i] * inv_fac[d-i] * \
+                       (-1 if (d-i) & 1 else 1) for i in range(d+1)]
+
+        # Let B_i = b+ir for each i in [-d, d], except we declare B_i = 1
+        # if i is "bad".
+        B = [(R(1) if i in bad else R(b + i*r)) for i in range(-d, d+1)]
+
+        # compute partial products of B_i
+        products = [R(1)]
+        for z in B:
+            products.append(products[-1] * z)
+
+        # compute partial products of inverses of B_i
+        inverse = ~(products[-1])
+        inv_products = [inverse]
+        for z in reversed(B):
+            inv_products.append(inv_products[-1] * z)
+        inv_products.reverse()
+        
+        # compute inverses of each B_i
+        B_inv = [inv_products[i] * products[i-1] \
+                 for i in range(1, len(products))]
+
+        # replace the "bad" B_i with 0, so they don't contribute to the
+        # convolution at all
+        for i in bad:
+            B_inv[i+d] = R(0)
+
+        # kernel is the polynomial we need to multiply by in the convolution
+        S = PolynomialRing(R, "x")
+        kernel = S(B_inv, check=False)
+
+        # coefficients for postprocessing step
+        output_scale = [products[k+d+1] * inv_products[k] for k in range(d+1)]
+
+        # save everything useful to the cache
+        cache[d, b, r] = (bad, input_scale, S, kernel, output_scale)
+    
+
+    # preprocessing step:
+    # multiply each P(ir) by \prod_{j \neq i} r (i-j)^{-1}
+    values = [values[i] * input_scale[i] for i in range(d+1)]
+
+    # do the convolution
+    H = S(values, check=False) * kernel
+    output = H.list()[d:(2*d+1)]
+    # (need to zero-pad output up to length d+1)
+    output.extend([R(0)] * (d+1 - len(output)))
+
+    # postprocessing:
+    # multiply each output by \prod_{j=0}^d B_{k-j}
+    output = [output[k] * output_scale[k] for k in range(d+1)]
+
+
+    # Now need to deal with "bad" indices. This shouldn't happen too often,
+    # especially when $d$ is much smaller than $p$.
+
+    # The first problem is that in the postprocessing step we left out
+    # multiplying by the bad factors. Now we put them back in.
+    for i in bad:
+        factor = b + i*r
+        for k in range(max(i, 0), min(d+i, d) + 1):
+            output[k] = output[k] * factor
+
+    # The second problem is that in the polynomial S(B_inv) the monomials
+    # corresponding to the bad i were completely wrong.
+    for i in bad:
+        if i >= 0:
+            correction = ([R(0)] * i) + [values[k-i] * output_scale[k] \
+                          for k in range(i, d+1)]
+        else:
+            correction = [values[k-i] * output_scale[k] \
+                          for k in range(0, i+d+1)] + ([R(0)] * (-i))
+
+        # again, we left out all the bad factors; put them back in
+        for ii in bad:
+            if ii != i:
+                factor = b + ii*r
+                for k in range(max(ii, 0), min(d+ii, d) + 1):
+                    correction[k] = correction[k] * factor
+
+        # add correction term to output
+        for k in range(d+1):
+            output[k] = output[k] + correction[k]
+
+    return output
+
+
+
+def naive_shift_evaluation_points(R, values, b, r, p):
+    r"""
+    Same as shift_evaluation_points, but uses naive algorithm (naive
+    Lagrange interpolation, followed by naive evaluation).
+
+    For testing/debugging purposes only.
+    """
+    d = len(values) - 1
+    S, x = PolynomialRing(R, "x").objgen()
+    total = S(0)
+    for j in range(d+1):
+        prod = S(1)
+        for i in range(d+1):
+            if i != j:
+                prod = prod * (x - i*r) / (r*(j-i))
+        total = total + values[j] * prod
+
+    return [total(b + i*r) for i in range(d+1)]
+
+
+
+def matrix_polys(A, R, n, b, r, q, p, cache):
+    r"""
+    INPUT:
+      p -- a prime
+      R -- coefficient ring, of the form $\ZZ/p^m \ZZ$
+      n -- a positive integer
+      A -- an $n \times n$ matrix of linear polys over $R$ (stored as a list
+           of lists)
+      b -- an integer
+      r -- a positive integer
+      q -- a positive integer
+      cache -- dictionary (passed through to shift_evaluation_points)
+
+    OUTPUT:
+      Let $B(x) = A(x+1) A(x+2) ... A(x+r-1) A(x+r)$, a matrix of degree
+      $r$ polynomials.
+
+      Output an $n \times n$ matrix, the $(i, j)$-entry is a list of the
+      $(i, j)$-entries of the evaluations $B(b), B(b+q), \ldots, B(b+rq)$;
+      i.e. it evaluates $B(x)$ at $r+1$ points spaced out by $q$.
+
+    ALGORITHM:
+      Umm, it's probably the one described in the Bostan/Gaudry/Schost paper,
+      although I haven't looked at that part of the paper too closely yet.
+
+    EXAMPLES:
+    Run it on some random input, and compare to a naive implementation:
+        sage: from sage.schemes.elliptic_curves.monsky_washnitzer import matrix_polys, naive_matrix_polys
+        sage: R = Integers(7**3)
+        sage: S.<x> = PolynomialRing(R, "x")
+        sage: b = 5
+        sage: for r in range(1, 14):
+        ...       A = [[R.random_element() * x + R.random_element() \
+        ...             for i in range(3)] for j in range(3)]
+        ...       X = matrix_polys(A, R, 3, b, r, 3, 7, {})
+        ...       Y = naive_matrix_polys(A, R, 3, b, r, 3, 7)
+        ...       if X != Y: print "failed"
+
+    """
+    assert r >= 1
+    assert r < 2*p
+
+    if r == 1:
+        # base case; i.e. evaluate original polys at x = b+1, b+q+1
+        return [[[A[j][i](b+1), A[j][i](b+q+1)] for i in range(n)]
+                for j in range(n)]
+
+    # let r' = floor(r/2)
+    half_r = r // 2
+    
+    # evaluate A(x+1) A(x+2) ... A(x+r')
+    # at x = b, b + q, ..., b + r'q.
+    X = matrix_polys(A, R, n, b, half_r, q, p, cache)
+
+    # shift to obtain evaluations of A(x+r'+1) A(x+r'+2) ... A(x+2r')
+    # at x = b, b + q, ..., b + r'q.
+    shift1 = [[shift_evaluation_points(R, values, half_r, q, p, cache) \
+               for values in row] for row in X]
+
+    # shift to obtain evaluations of A(x+1) A(x+2) ... A(x+r')
+    # at x = b + (r'+1)q, b + (r'+2)q, ..., b + (2r'+1)q
+    shift2 = [[shift_evaluation_points(R, values, (half_r+1)*q, q, p, cache) \
+               for values in row] for row in X]
+
+    # shift to obtain evaluations of A(x+r'+1) A(x+r'+2) ... A(x+2r')
+    # at x = b + (r'+1)q, b + (r'+2)q, ..., b + (2r'+1)q
+    shift3 = [[shift_evaluation_points(R, values, half_r, q, p, cache) \
+               for values in row] for row in shift2]
+
+    # multiply matrices to obtain evaluations of A(x+1) A(x+2) ... A(x+2r')
+    # at x = b, b + q, ..., b + r'q.
+    prod1 = [[[sum([X[i][m][k] * shift1[m][j][k] for m in range(n)]) \
+               for k in range(half_r+1)] for j in range(n)] for i in range(n)]
+               
+    # multiply matrices to obtain evaluations of A(x+1) A(x+2) ... A(x+2r')
+    # at x = b + (r'+1)q, ..., b + (2r'+1)q.
+    prod2 = [[[sum([shift2[i][m][k] * shift3[m][j][k] for m in range(n)]) \
+               for k in range(half_r+1)] for j in range(n)] for i in range(n)]
+               
+    if r & 1 == 0:
+        # if r is even, then the polys are the right degree, but we have one
+        # extra evaluation point (namely x = b + (r+1)q). Throw it away.
+        for i in range(n):
+            for j in range(n):
+                prod2[i][j].pop()
+                prod1[i][j].extend(prod2[i][j])
+        output = prod1
+        
+    else:
+        # if r is odd, then we have precisely the right number of evaluation
+        # points, but our polynomials are short by one term, namely A(x+r).
+        # Add in the effect of that term.
+        for i in range(n):
+            for j in range(n):
+                prod1[i][j].extend(prod2[i][j])
+
+        output = \
+               [[[sum([prod1[i][m][k] * A[m][j](b+k*q+r) for m in range(n)]) \
+                  for k in range(r+1)] for j in range(n)] for i in range(n)]
+
+    return output
+
+
+
+def naive_matrix_polys(A, R, n, b, r, q, p):
+    r"""
+    Same as matrix_polys, but uses a naive algorithm.
+
+    For testing/debugging purposes only.
+    """
+    result = [[[] for _ in range(n)] for _ in range(n)]
+    for i in range(r+1):
+        prod = MatrixSpace(R, n, n).identity_matrix()
+        for j in range(1, r+1):
+            A_eval = matrix(R, [[f(b+q*i+j) for f in row] for row in A])
+            prod = prod * A_eval
+        for k1 in range(n):
+            for k2 in range(n):
+                result[k1][k2].append(prod[k1][k2])
+    return result
+            
+
+
+
+def reduce_one(T, s, n, M, D, R):
+    r"""
+    A single step reduction.
+    
+    INPUT:
+      R -- coefficient ring, of the form $\ZZ/p^m \ZZ$
+      n -- a positive integer
+      T -- an $n$-tuple of elements of $R$.
+      M -- an $n \times n$ matrix of linear polys over $R$ (stored as a list
+           of lists)
+      D -- a linear polynomial over $R$.
+
+    OUTPUT:
+      Let $A(x) = M(x)/D(x)$.
+      Return value is the matrix-vector product $A(s) T$, as a list.
+      
+    """
+    d = D(s).lift()
+    return [R(sum([M[i][j](s) * T[j] for j in range(n)]).lift() / d) \
+            for i in range(n)]
+
+
+
+def reduce_sequence(T, d, s, n, M, D, p, R):
+    r"""
+    Reduces a sequence of terms.
+    
+    INPUT:
+      p -- a prime
+      R -- coefficient ring, of the form $\ZZ/p^m \ZZ$
+      T -- a nonempty list of $n$-tuples of elements of $R$. Denote by $L$
+           the length of $T$.
+      d -- a list of positive integers, of length $L$; they must be in
+           strictly ascending order
+      s -- an integer >= 0
+      n -- an integer >= 1
+      M -- an $n \times n$ matrix of linear polynomials over $R$,
+           stored as a list of lists
+      D -- a linear polynomial over $R$
+
+    PRECONDITIONS:
+       The largest $d_i$ must be less than $\sqrt p$.
+       (or something like that... todo: check this...)
+
+    OUTPUT:
+       Let $A(x) = M(x)/D(x)$, and let $B(t)$ denote the matrix product
+         $$ A(s+1) A(s+2) \cdots A(s+t). $$
+
+       This function returns the sum
+         $$ \sum_{i=1}^L  B(d_i) T_i. $$
+
+       In other words, treating each term $T_i$ as appearing in position
+       $s + d_i$, it applies the reduction matrices to reduce them all down to
+       a single $n$-tuple in position $s$.
+
+       Running time is something like
+          $O(M(e^{1/2}) + L M(e^{1/4}))$
+       where $e = \max(d_i)$, and where $M(k)$ is the time for degree $k$
+       polynomial multiplication over $R$.
+
+    """
+    # L = number of terms to reduce
+    L = len(T)
+    # e = total distance to reduce
+    e = d[-1]
+
+    if e <= 30:
+        # naive algorithm when the reduction length is small enough
+        output = [R(0)] * n
+        for i in range(L):
+            term = T[i]
+            for t in range(d[i], 0, -1):
+                term = reduce_one(term, s + t, n, M, D, R)
+            for j in range(n):
+                output[j] = output[j] + term[j]
+        return output
+
+    # We do the reduction by splitting into segments. There are $w+1$
+    # segments of length $w$, where $w$ is the smallest integer such
+    # that $w(w+1) \geq e$ (so in particular $w \approx \sqrt e$).
+    w = ceil(sqrt(e))
+    while w*(w+1) >= e:
+        w = w - 1
+    if w*(w+1) < e:
+        w = w + 1
+
+    # The BGS method only works if the terms are on segment boundaries.
+    # So first we (recursively) reduce each term separately to make it
+    # lie on a segment boundary.
+    # This is the step that takes time $O(L M(e^{1/4}))$.
+    new_T = []
+    new_d = []
+    for i in range(L):
+        this_d = d[i]
+        reduced_d = w * (this_d // w)
+        if reduced_d == this_d:
+            reduced_T = T[i]
+        else:
+            reduced_T = reduce_sequence([T[i]], [this_d - reduced_d],
+                                        s + reduced_d, n, M, D, p, R)
+        new_T.append(reduced_T)
+        new_d.append(reduced_d)
+    T = new_T
+    d = new_d
+
+    # Now all terms are on segment boundaries. Next we compute a reduction
+    # matrix for each segment, using BGS method.
+    cache = {}
+    M_segments = matrix_polys([[M[i][j] for j in range(n)] for i in range(n)],
+                              R, n, s, w, w, p, cache)
+    D_segments = matrix_polys([[D]], R, 1, s, w, w, p, cache)[0][0]
+
+    # Now we use the reduction matrices to sweep up all the terms.
+    total = [R(0), R(0), R(0)]
+    for segment in range(w+1, -1, -1):
+        # add in terms at the current segment boundary
+        while len(d) and (d[-1] == segment*w):
+            for j in range(n):
+                total[j] = total[j] + T[-1][j]
+            T.pop()
+            d.pop()
+
+        # apply reduction matrix
+        if segment > 0:
+            M_here = [[M_segments[i][j][segment-1] \
+                       for j in range(n)] for i in range(n)]
+            D_here = D_segments[segment-1].lift()
+            total = [R(sum([M_here[i][j] * total[j]
+                            for j in range(n)]).lift() / D_here) \
+                     for i in range(n)]
+    
+    return total
+
+
+
+def starting_differentials(R, a, b, p, M):
+    r"""
+    Computes the initial terms in the expansion of $\sigma(x^i dx/y)$
+    (where $\sigma$ is frobenius) to which the cohomology reduction
+    algorithm will be applied.
+    
+    INPUT:
+      p -- a prime $\geq 5$
+      R -- coefficient ring, of the form $\ZZ/p^m \ZZ$
+      a, b -- elements of $R$, coefficients of an elliptic curve
+              $y^2 = x^3 + ax + b$
+      M -- number of terms being computed (see ``Kedlaya's algorithm in
+           larger characteristic'' for the exact meaning)
+
+    OUTPUT:
+      A list of length $M$; the $j$th entry is a list of length $3j+1$;
+      its $m$th entry is the coefficient of
+         $ x^{p(m+i+1)-1} y^{-p(2j+1)+1} dx/y $
+      in the expansion of $\sigma(x^i dx/y)$.
+
+      Complexity is $O(M^2)$ operations in $R$.
+
+    """
+    # compute B[k][j] = binomial(k, j) in R, for 0 <= j <= k < 2M
+    B = [[R(1)]]
+    for k in range(1, 2*M):
+        B.append([R(1)] + [B[-1][j] + B[-1][j-1] for j in range(1,k)] + [R(1)])
+
+    # compute temp[k] = binomial(2k, k) / 4^k in R, for 0 <= k < M
+    four_to_minus_1 = R(4)**(-1)
+    four_to_minus_k = R(1)
+    temp = [R(1)]
+    for k in range(1, M):
+        four_to_minus_k = four_to_minus_k * four_to_minus_1
+        temp.append(four_to_minus_k * B[2*k][k])
+
+    # compute constants
+    # C[j] = (-1)^j p \sum_{k=j}^{M-1} 4^{-k} (2k \choose k) (k \choose j)
+    # in R, for 0 <= j < M
+    C = [p * sum([temp[k] * B[k][j] for k in range(j, M)]) for j in range(M)]
+    for j in range(1, M, 2):
+        C[j] = -C[j]
+
+    # compute A[j][m] = u^m coefficient of (u^3 + a*u + b)^j
+    # for 0 <= j < M
+    Ru, u = PolynomialRing(R, "u").objgen()
+    Q = u**3 + a*u + b
+    Q_to_j = Ru(1)
+    A = [[R(1)]]
+    for j in range(1, M):
+        Q_to_j = Q_to_j * Q
+        A.append(Q_to_j.list())
+
+    # compute products C[j] * A[j][m]
+    return [[C[j] * x for x in A[j]] for j in range(M)]
+
+
+
+
+def matrix_of_frobenius_alternate(a, b, p, N):
+    r"""
+    Computes the matrix of Frobenius on Monsky-Washnitzer cohomology,
+    with respect to the basis $(dx/y, x dx/y)$.
+
+    INPUT:
+        a, b -- rational numbers, coefficients of elliptic curve
+                $y^2 = x^3 + ax + b$
+        p -- a prime for which the curve has good reduction
+        N -- desired output precision for matrix (i.e. modulo $p^N$).
+
+    OUTPUT:
+        2x2 matrix of frobenius, with coefficients in Integers(p^N)
+
+    TODO:
+    Haven't yet checked yet exactly which combinations of $p$ and $N$ are
+    valid. If $p$ is much larger than $N$, there shouldn't be a problem.
+
+    EXAMPLES:
+        sage: from monsky_washnitzer import matrix_of_frobenius_alternate
+        sage: M = matrix_of_frobenius_alternate(-1, 1/4, 101, 3); M
+        [824362 606695]
+        [641451 205942]
+        sage: M.det()
+        101
+        sage: M.trace()
+        3
+    
+    """
+    assert p > N    # todo: check this is the right condition....
+
+    # choose minimal M such that M >= floor(log_p(2M+1)) + N
+    M = N
+    while M < floor(log(2*M+1, p)) + N:
+      M = M + 1
+
+    # I haven't done the precision analysis properly yet, so for now
+    # we work to precision 3M, with M digits at the bottom and M spare
+    # digits at the top, which should be pretty safe.
+    R = Integers(p**(3*M))
+    R_output = Integers(p**N)
+    S, x = PolynomialRing(R, "x").objgen()
+    a = R(a)
+    b = R(b)
+    
+    differentials = starting_differentials(R, a, b, p, M)
+
+    # multiply through by p^M for safety, need to come back and fix this...
+    differentials = [[item * p**M for item in row] for row in differentials]
+
+    output = []
+
+    for i in range(2):
+        column = []
+        
+        for j in range(M):
+            t = p*(2*j+1)-1
+
+            terms = [R(0)] if i == 1 else []
+            terms.extend(differentials[j])
+            terms = [[y, R(0), R(0)] for y in terms]
+
+            # reduce this row to the left
+            triple = reduce_sequence(
+              terms, [p*(m+1)-1 for m in range(len(terms))],
+              0, 3,
+              [[S(0),            S(0),            -2*b*x          ],
+               [2*x - 3*t + 3,   S(0),            -a*(2*x - t + 1)],
+               [S(0),            2*x - 3*t + 3,   S(0)            ]],
+              2*x - 3*t + 3,
+              p, R)
+
+            # do one reduction to make it only two terms instead of three
+            triple = reduce_one(triple, t//2, 3,
+              [[ 9*b*(6*x-5), 2*a*a*(6*x-5), -3*a*b*(6*x-5)],
+               [-6*a*(6*x-7),   9*b*(6*x-7),  2*a*a*(6*x-7)],
+               [        S(0),          S(0),           S(0)]],
+              (27*b**2 + 4*a**3)*(2*x-1), R)
+
+            column.append([triple[0], triple[1]])
+
+        # reduce the column to the bottom
+        result = reduce_sequence(
+              column, [(p*(2*m+1)-3)//2 for m in range(len(column))],
+              0, 2,
+              [[ 9*b*(6*x-5), 2*a*a*(6*x-5)],
+               [-6*a*(6*x-7),   9*b*(6*x-7)]],
+              (27*b**2 + 4*a**3)*(2*x-1),
+              p, R)
+
+        # undo stupid p^M adjustment, and reduce down to Z/p^N
+        for k in range(2):
+            result[k] = R_output(result[k].lift() / p**M)
+
+        output.append(result)
+
+    return matrix(R_output, 2, 2, [output[0][0], output[1][0],
+                                   output[0][1], output[1][1]])
+            
+
 ### end of file