require 6 to be a unit in elliptic-curve point addition over non-fields

src/sage/schemes/elliptic_curves/ell_point.py

@@ -252,7 +252,7 @@ def _add_(self, other):
         prime divisors, for which the result is computed using the "old",
         much simpler formulas for fields.) ::
 
-            sage: N = ZZ(randrange(2, 10**5))
+            sage: N = randrange(10**4) * 6 + choice([5, 7])  # coprime to 6
             sage: E = None
             sage: while True:
             ....:     try:
@@ -272,19 +272,19 @@ def _add_(self, other):
             ....:     if xs:
             ....:         pts.append(E(choice(xs), y, z))
             sage: P, Q = pts
-            sage: R = P + Q  # not tested (:issue:`39191`)
-            sage: for d in N.divisors():  # not tested (:issue:`39191`)
+            sage: R = P + Q
+            sage: for d in N.divisors():
             ....:     if d > 1:
             ....:         assert R.change_ring(Zmod(d)) == P.change_ring(Zmod(d)) + Q.change_ring(Zmod(d))
         """
-        if self.is_zero():
-            return other
-        if other.is_zero():
-            return self
-
         E = self.curve()
         R = E.base_ring()
 
+        # According to https://cr.yp.to/bib/1987/lenstra-ecnta.pdf, §3,
+        # the formulas require 6 to be a unit. See #39191 for details.
+        if not R(6).is_unit():
+            raise NotImplementedError('addition of elliptic-curve points over non-fields is only supported when 6 is a unit')
+
         # We handle Euclidean domains modulo principal ideals separately.
         # Important special cases of this include quotient rings of the
         # integers as well as of univariate polynomial rings over fields.
@@ -360,6 +360,9 @@ def _add_(self, other):
         # Below, we simply try random linear combinations until we
         # find a good choice. Is there a general method that doesn't
         # involve guessing?
+        # Answer: Yes.
+        # See pages 7-8 of Lenstra's "Elliptic Curves and Number-Theoretic Algorithms".
+        # https://cr.yp.to/bib/1987/lenstra-ecnta.pdf
 
         pts = [vector(R, pt) for pt in pts]
         for _ in range(1000):