implement Bosma-Lenstra formulas (corrected version from Best)

src/doc/en/reference/references/index.rst

@@ -453,6 +453,9 @@ REFERENCES:
 .. [BeBo2009] Olivier Bernardi and Nicolas Bonichon, *Intervals in Catalan
               lattices and realizers of triangulations*, JCTA 116 (2009)
 
+.. [Best2021] Alex J. Best: Tools and Techniques for Rational Points on Curves.
+              PhD Thesis, Boston University, 2021.
+
 .. [BBGL2008] \A. Blondin Massé, S. Brlek, A. Garon, and S. Labbé,
               Combinatorial properties of f -palindromes in the
               Thue-Morse sequence. Pure Math. Appl.,
@@ -892,6 +895,10 @@ REFERENCES:
             Anal. Appl. 15 (1994) 804-823.
             :doi:`10.1137/S0895479892230031`
 
+.. [BL1995] W. Bosma, H.W. Lenstra: Complete Systems of Addition Laws for
+            Elliptic Curves. Journal of Number Theory, volume 53, issue 2,
+            pages 229-240. 1995.
+
 .. [BHMPW20a] Tom Braden, June Huh, Jacob P. Matherne, Nicholas Proudfoot,
              and Botong Wang, *A semi-small decomposition of the Chow
              ring of a matroid*, :arxiv:`2002.03341` (2020).

src/sage/schemes/elliptic_curves/addition_formulas_ring.py

@@ -0,0 +1,48 @@
+
+def add(E, P, Q):
+    r"""
+    Addition formulas for elliptic curves over general rings
+    with trivial Picard group.
+
+    REFERENCES:
+
+    These formulas were derived by Bosma and Lenstra [BL1995]_,
+    with corrections by Best [Best2021]_ (Appendix A).
+    """
+    a1, a2, a3, a4, a6 = E.a_invariants()
+    b2, b4, b6, b8 = E.b_invariants()
+
+    assert P in E
+    assert Q in E
+    X1, Y1, Z1 = P
+    X2, Y2, Z2 = Q
+
+    #TODO: I've made a half-hearted attempt at simplifying the formulas
+    # by caching common subexpressions. This could almost certainly be
+    # sped up significantly with some more serious optimization effort.
+
+    XYdif = X1*Y2 - X2*Y1
+    XYsum = X1*Y2 + X2*Y1
+    XZdif = X1*Z2 - X2*Z1
+    XZsum = X1*Z2 + X2*Z1
+    YZdif = Y1*Z2 - Y2*Z1
+    YZsum = Y1*Z2 + Y2*Z1
+
+    a1sq, a2sq, a3sq, a4sq = (a**2 for a in (a1, a2, a3, a4))
+
+    X31 = XYdif*YZsum+XZdif*Y1*Y2+a1*X1*X2*YZdif+a1*XYdif*XZsum-a2*X1*X2*XZdif+a3*XYdif*Z1*Z2+a3*XZdif*YZsum-a4*XZsum*XZdif-3*a6*XZdif*Z1*Z2
+
+    Y31 = -3*X1*X2*XYdif-Y1*Y2*YZdif-2*a1*XZdif*Y1*Y2+(a1sq+3*a2)*X1*X2*YZdif-(a1sq+a2)*XYsum*XZdif+(a1*a2-3*a3)*X1*X2*XZdif-(2*a1*a3+a4)*XYdif*Z1*Z2+a4*XZsum*YZdif+(a1*a4-a2*a3)*XZsum*XZdif+(a3sq+3*a6)*YZdif*Z1*Z2+(3*a1*a6-a3*a4)*XZdif*Z1*Z2
+
+    Z31 = 3*X1*X2*XZdif-YZsum*YZdif+a1*XYdif*Z1*Z2-a1*XZdif*YZsum+a2*XZsum*XZdif-a3*YZdif*Z1*Z2+a4*XZdif*Z1*Z2
+
+    yield (X31, Y31, Z31)
+
+    X32 = Y1*Y2*XYsum+a1*(2*X1*Y2+X2*Y1)*X2*Y1+a1sq*X1*X2**2*Y1-a2*X1*X2*XYsum-a1*a2*X1**2*X2**2+a3*X2*Y1*(YZsum+Y2*Z1)+a1*a3*X1*X2*YZdif-a1*a3*XYsum*XZdif-a4*X1*X2*YZsum-a4*XYsum*XZsum-a1sq*a3*X1**2*X2*Z2-a1*a4*X1*X2*(X1*Z2+XZsum)-a2*a3*X1*X2**2*Z1-a3sq*X1*Z2*(Y2*Z1+YZsum)-3*a6*XYsum*Z1*Z2-3*a6*XZsum*YZsum-a1*a3sq*X1*Z2*(XZsum+X2*Z1)-3*a1*a6*X1*Z2*(XZsum+X2*Z1)-a3*a4*(X1*Z2+XZsum)*X2*Z1-b8*YZsum*Z1*Z2-a1*b8*X1*Z1*Z2**2-a3**3*XZsum*Z1*Z2-3*a3*a6*(XZsum+X2*Z1)*Z1*Z2-a3*b8*Z1**2*Z2**2
+
+    Y32 = Y1**2*Y2**2+a1*X2*Y1**2*Y2+(a1*a2-3*a3)*X1*X2**2*Y1+a3*Y1**2*Y2*Z2-(a2sq-3*a4)*X1**2*X2**2+(a1*a4-a2*a3)*(2*X1*Z2+X2*Z1)*X2*Y1+(a1sq*a4-2*a1*a2*a3+3*a3sq)*X1**2*X2*Z2-(a2*a4-9*a6)*X1*X2*XZsum+(3*a1*a6-a3*a4)*(XZsum+X2*Z1)*Y1*Z2+(3*a1sq*a6-2*a1*a3*a4+a2*a3sq+3*a2*a6-a4sq)*X1*Z2*(XZsum+X2*Z1)+(3*a2*a6-a4sq)*X2*Z1*(2*X1*Z2+Z1*X2)+(a1**3*a6-a1sq*a3*a4+a1*a2*a3sq-a1*a4sq+4*a1*a2*a6-a3**3-3*a3*a6)*Y1*Z1*Z2**2+(a1**4*a6-a1**3*a3*a4+5*a1sq*a2*a6+a1sq*a2*a3sq-a1*a2*a3*a4-a1*a3**3-3*a1*a3*a6-a1sq*a4sq+a2sq*a3sq-a2*a4sq+4*a2sq*a6-a3 **2*a4-3*a4*a6)*X1*Z1*Z2**2+(a1sq*a2*a6-a1*a2*a3*a4+3*a1*a3*a6+a2sq*a3sq-a2*a4sq+4*a2sq*a6-2*a3sq*a4-3*a4*a6)*X2*Z1**2*Z2+(a1**3*a3*a6-a1sq*a3sq*a4+a1sq*a4*a6+a1*a2*a3**3+4*a1*a2*a3*a6-2*a1*a3*a4sq+a2*a3sq*a4+4*a2*a4*a6-a3**4-6*a3 **2*a6-a4**3-9*a6**2)*Z1**2*Z2**2
+
+    Z32 = 3*X1*X2*XYsum+Y1*Y2*YZsum+3*a1*X1**2*X2**2+a1*(2*X1*Y2+Y1*X2)*Y1*Z2+a1sq*X1*Z2*(2*X2*Y1+X1*Y2)+a2*X1*X2*YZsum+a2*XYsum*XZsum+a1**3*X1**2*X2*Z2+a1*a2*X1*X2*(2*X1*Z2+X2*Z1)+3*a3*X1*X2**2*Z1+a3*Y1*Z2*(YZsum+Y2*Z1)+2*a1*a3*X1*Z2*YZsum+2*a1*a3*X2*Y1*Z1*Z2+a4*XYsum*Z1*Z2+a4*XZsum*YZsum+(a1sq*a3+a1*a4)*X1*Z2*(XZsum+X2*Z1)+a2*a3*X2*Z1*(2*X1*Z2+X2*Z1)+a3sq*Y1*Z1*Z2**2+(a3sq+3*a6)*YZsum*Z1*Z2+a1*a3sq*(2*X1*Z2+X2*Z1)*Z1*Z2+3*a1*a6*X1*Z1*Z2**2+a3*a4*(XZsum+X2*Z1)*Z1*Z2+(a3**3+3*a3*a6)*Z1**2*Z2**2
+
+    yield (X32, Y32, Z32)
+

src/sage/schemes/elliptic_curves/ell_point.py

@@ -4,7 +4,7 @@
 
 The base class :class:`EllipticCurvePoint` provides support for
 points on elliptic curves defined over general rings, including
-somewhat generic addition formulas. (Not implemented yet.)
+generic addition formulas.
 
 The derived classes :class:`EllipticCurvePoint_field` and its
 child classes :class:`EllipticCurvePoint_number_field` and
@@ -166,8 +166,75 @@ def _add_(self, other):
         Add this point to another point on the same elliptic curve.
 
         This method computes point additions for fairly general rings.
+
+        ALGORITHM:
+
+        Formulas due to Bosma and Lenstra [BL1995]_ with corrections
+        by Best [Best2021]_ (Appendix A).
+        See :mod:`sage.schemes.elliptic_curves.addition_formulas_ring`.
+
+        EXAMPLES::
+
+            sage: N = 1113121
+            sage: E = EllipticCurve(Zmod(N), [1,0])
+            sage: R1 = E(301098, 673883, 644675)
+            sage: R2 = E(411415, 758555, 255837)
+            sage: R3 = E(983009, 342673, 207687)
+            sage: R1 + R2 == R3
+            True
+
+        Checks that :trac:`15964` is fixed::
+
+            sage: N = 1715761513
+            sage: E = EllipticCurve(Integers(N),[3,-13])
+            sage: P = E(2,1)
+            sage: LCM([2..60])*P
+            Traceback (most recent call last):
+            ...
+            ZeroDivisionError: Inverse of 1520944668 does not exist
+            (characteristic = 1715761513 = 26927*63719)
+
+            sage: N = 35
+            sage: E = EllipticCurve(Integers(N),[5,1])
+            sage: P = E(0,1)
+            sage: LCM([2..6])*P
+            Traceback (most recent call last):
+            ...
+            ZeroDivisionError: Inverse of 28 does not exist
+            (characteristic = 35 = 7*5)
         """
-        raise NotImplementedError
+        if self.is_zero():
+            return other
+        if other.is_zero():
+            return self
+
+        E = self.curve()
+        R = E.base_ring()
+
+        from sage.schemes.elliptic_curves.addition_formulas_ring import add
+        from sage.modules.free_module_element import vector
+
+        pts = []
+        for pt in filter(any, add(E, self, other)):
+            if R.one() in R.ideal(pt):
+                return E.point(pt)
+            pts.append(pt)
+        assert len(pts) == 2, 'bug in elliptic-curve point addition'
+
+        #TODO: If the base ring has trivial Picard group, it is known
+        # that some linear combination of the two vectors is a valid
+        # projective point (whose coordinates generate the unit ideal).
+        # Below, we simply try random linear combinations until we
+        # find a good choice. Is there a general method that doesn't
+        # involve guessing?
+
+        pts = [vector(R, pt) for pt in pts]
+        for _ in range(1000):
+            result = tuple(sum(R.random_element() * pt for pt in pts))
+            if R.one() in R.ideal(result):
+                return E.point(result)
+
+        assert False, 'bug: failed to compute elliptic-curve point addition'
 
     def _sub_(self, other):
         """