implement elliptic-curve point addition for general rings

yyyyx4 commented 2 years ago

Currently, point addition for elliptic curves is only defined in EllipticCurvePoint_field and its children. This patch implements point arithmetic for EllipticCurvePoint, while also keeping the specialized (faster) implementation for fields. This is useful for symbolic computations, such as the ones that show up in a basic version of Schoof's algorithm (where we compute in a polynomial ring modulo the ideal generated by a division polynomial and the curve equation).

Another implication of this is that we can remove part of the "temporary" hack for elliptic-curve points over ℤ/n that was introduced in #1975.

Depends on #34417 Depends on #34463

CC: @JohnCremona @roed314 @edgarcosta @alexjbest

Component: elliptic curves

Author: Lorenz Panny

Branch/Commit: public/generalize_elliptic_curve_point_addition_to_rings @ 940f8a2

Issue created by migration from https://trac.sagemath.org/ticket/33228

fchapoton commented 2 years ago

comment:1

some failing doctests, see patchbot report

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 2 years ago

Branch pushed to git repo; I updated commit sha1. New commits:

915b725 make (P in E) return False when P is defined over an extension field

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 2 years ago

Changed commit from 4d67190 to 915b725

fchapoton commented 2 years ago

comment:4

Elliptic gurus, please say your word on these changes!

JohnCremona commented 2 years ago

comment:5

This is interesting -- as well as being the author of the "temporary hack", I also came across this issue when trying to implement Schoof's algorithm from basics wih a student. I'll take a look.

edgarcosta commented 2 years ago

comment:7

In case you are interested, Alex Best corrected some typos Bosma–Lenstra paper formulas for addition law for elliptic curves over a general ring, see Appendix A of https://open.bu.edu/bitstream/handle/2144/43150/Best_bu_0017E_16371.pdf

I bet he must have them implemented somewhere...

yyyyx4 commented 2 years ago

comment:8

Replying to @edgarcosta:

In case you are interested, Alex Best corrected some typos Bosma–Lenstra paper formulas for addition law for elliptic curves over a general ring, see Appendix A of https://open.bu.edu/bitstream/handle/2144/43150/Best_bu_0017E_16371.pdf

Interesting, thanks! Generalizing the formulas to not fail when non-units are encountered was a next step on my list (and the original motivation to move + up in the class hierarchy). I'll try this out as soon as I find time.

Replying to @JohnCremona:

This is interesting -- as well as being the author of the "temporary hack", I also came across this issue when trying to implement Schoof's algorithm from basics wih a student. I'll take a look.

Sadly, this change alone does not yet allow implementing a straightforward Schoof (e.g., because the current addition formulas attempt to divide by zero divisors rather frequently in that context). But we'll get there! :-)

JohnCremona commented 2 years ago

comment:9

I reealise that it may be annoying to make comments about chunks of code which were just copied over but I would change the lines

        x1, y1 = self[0], self[1]
        x2, y2 = right[0], right[1]
        if x1 == x2 and y1 == -y2 - a1*x2 - a3:

to

        x1, y1, _ = self
        x2, y2, _ = right
        if x1 == x2 and y1 != y2:

where the first two are just python while the 3rd is simpler and mathematically equivalent.

However, are you sure that the z-coordinate will always be 1 for nonzero points over non-fields? I will keep reading...

JohnCremona commented 2 years ago

comment:10

Another simplification: after the preceding check, if x1==x2 then also y1==y2 (unless you do not want me to assume that quadratics over non-fields have no more than 2 roots, in which case my previous suggestion may not be valid either):

    if x1==x2:
       num = 3*x1*x1 + 2*a2*x1 + a4 - a1*y1
       den = 2*y1 + a1*x1 + a3
    else:
       num = y1-y2
       den = x1-x2
   try:
      m = num/den
   except ZeroDivisionError:
              R = E.base_ring()
              if R.is_finite():
                  N = R.characteristic()
                  N1 = N.gcd(Integer(den))
                  N2 = N//N1
                  raise ZeroDivisionError("Inverse of %s does not exist (characteristic = %s = %s*%s)" % (den, N, N1, N2))
              else:
                   raise ZeroDivisionError("Inverse of %s does not exist" % (den))

Apart from these (which I know don't have much to do with your chenges) I think this looks good. Can you add at least one test to show that what used to work (over Z/NZ say) still does, and also something which used not to work but now does?

yyyyx4 commented 2 years ago

Dependencies: #34417

yyyyx4 commented 2 years ago

Description changed:

--- 
+++ 
@@ -1,3 +1,3 @@
-Currently, point addition for elliptic curves is only defined in `EllipticCurvePoint_field` and its children. This patch moves it up to `EllipticCurvePoint`. This is useful for symbolic computations, such as the ones that show up in a basic version of Schoof's algorithm (where we compute in a polynomial ring modulo the ideal generated by a division polynomial and the curve equation).
+Currently, point addition for elliptic curves is only defined in `EllipticCurvePoint_field` and its children. This patch implements point arithmetic for `EllipticCurvePoint`, while also keeping the specialized (faster) implementation for fields. This is useful for symbolic computations, such as the ones that show up in a basic version of Schoof's algorithm (where we compute in a polynomial ring modulo the ideal generated by a division polynomial and the curve equation).

 Another implication of this is that we can remove part of the "temporary" hack for elliptic-curve points over ℤ/n that was introduced in #1975.

yyyyx4 commented 2 years ago

Changed commit from 915b725 to 2b4509e

yyyyx4 commented 2 years ago

comment:13

I've included Alex' version of the Bosma-Lenstra formulas to compute point additions over fairly general rings. One thing I'm not quite sure about is how to combine the results of the two formulas in case neither gives a valid projective point — in the current code, this is done by brute-forcing random linear combinations, but there must be a more algebraic way. (It's fairly straightforward for Bézout rings.)

New commits:

`020c9e5`	`move elliptic-curve point addition to superclass`
`52c622b`	`make (P in E) return False when P is defined over an extension field`
`c73444e`	`undo ancient hack for elliptic-curve points modulo composites`
`45ba391`	`update documentation to reflect this change`
`d0d6100`	`implement Bosma-Lenstra formulas (corrected version from Best)`
`f890d6c`	`add random test for point addition in ℤ/n`
`8b56f25`	`undo ancient hack for ECM examples`
`2b4509e`	`add changes to history`

yyyyx4 commented 2 years ago

Changed branch from public/move_elliptic_curve_point_addition_to_parent to public/generalize_elliptic_curve_point_addition_to_rings

roed314 commented 2 years ago

comment:14

Looks pretty good to me. Here are some questions/suggestions.

Are there other cases where we can use simpler formulas? For example, over a Euclidean domain you could use the field-formula and then clear denominators. And over Z/N you could try to use the field-formula and only fall back on the more general one when you hit a zero divisor. I haven't done any profiling, but the general formula looks pretty slow, so this change will probably slow down any toy ECM implementations out there.
Relatedly, I'm not sure exactly how to handle deprecation warnings for this change. Making arithmetic work when it didn't before is great, but this is one of the few cases where people may be relying on having a ZeroDivisionError rather than getting an actual point out. I know that my ECM demos would break with this change. Perhaps one way forward would be to print a deprecation warning whenever someone does arithmetic on an elliptic curve over Z/N with a pointer to this ticket so that at least they understand what has changed when their code just runs forever (as it might if they just do things until catching a ZeroDivisionError).

I'll also give John Cremona a chance to take a look at this, since he was commenting on an earlier version.

yyyyx4 commented 2 years ago

comment:15

Thanks for your comments.

Indeed, some quick experiments suggest a ≈10× slowdown. Perhaps we should try to adapt the standard formulas to quotients of Euclidean domains? This would cover both of your examples, and more.
Getting the ZeroDivisionError back should be achievable by explicitly constructing an EllipticCurvePoint_field(E, ...) rather than using E.point(...). The standard ZeroDivisionError error message is slightly less well-adjusted to the ECM use case, but it does contain the relevant non-unit so you can still easily factor the modulus from the exception string. Is this what we should recommend to users in a warning?

JohnCremona commented 2 years ago

comment:16

Since I was asked, I'll comment. Two points:

I myself have no need to ever try to define elliptic curves over general rings. The reason I put in the hack which pretended that elliptic curves over \ZZ/N were elliptic curves over fields, was to enable a simple-minded ECM to work for teaching purposes, making sure that the error message produced, when zero-divisors were encountered, displayed the nontrivial factor of N. And I did this after a complaint by Ken Ribet (eminent mathematician and AMS President 2017-2018) that the Sage code he used to teach exactly this to undergraduates at Berkeley no longer worked. It's a very good thing to have such people promoting Sage.

I think it could work to have some parameter to the elliptic curve constructor like assume_base_ring_is_field, which by default is False, but which could be set to True for this precise purpose -- with some doctests to show exactly how this could be used for such a simple ECM implementation. In (just about) all other situations, the default would do the mathematically correct thing.

Secondly, it is unacceptable to have a 10* slowdown for point operations over fields! The vast majority of point operations would be over fields. Surely the thing to do would be to put in the generic stuff introduced here for EllipticCurvePoint, while keeping the old field-only code (or whatever is fastest over fields) for the class EllipticCurvePoint_field?

yyyyx4 commented 2 years ago

comment:17

Replying to @JohnCremona:

Surely the thing to do would be to put in the generic stuff introduced here for EllipticCurvePoint, while keeping the old field-only code (or whatever is fastest over fields) for the class EllipticCurvePoint_field?

Of course; the patch does it like this already. The 10× difference is between the field-specific formulas and the general formulas, not between before and after. The only case which is now slower (but more correct!) than before is ℤ/n for composite n.

JohnCremona commented 2 years ago

comment:18

Replying to @yyyyx4:

Replying to @JohnCremona:

Surely the thing to do would be to put in the generic stuff introduced here for EllipticCurvePoint, while keeping the old field-only code (or whatever is fastest over fields) for the class EllipticCurvePoint_field?

Of course; the patch does it like this already. The 10× difference is between the field-specific formulas and the general formulas, not between before and after. The only case which is now slower (but more correct!) than before is ℤ/n for composite n.

Thanks for clarifying (and I know that I could also have looked at the code, but it is good to have this stated explicitly on the ticket).

This deals with my second point, at least!

roed314 commented 2 years ago

comment:19

Two plugin issues:

There's a pyflakes error: src/sage/schemes/elliptic_curves/ell_generic.py:56:1 'sage.rings.abc' imported but unused
The add function in addition_formulas_ring.py needs doctests.

There's also a doctest failure:

File "src/sage/schemes/elliptic_curves/ell_point.py", line 214, in sage.schemes.elliptic_curves.ell_point.EllipticCurvePoint._add_
Failed example:
    for d in N.divisors():
        if d > 1:
            assert R.change_ring(Zmod(d)) == P.change_ring(Zmod(d)) + Q.change_ring(Zmod(d))
Exception raised:
    Traceback (most recent call last):
      File "/home/sage-patchbot/sage/src/sage/doctest/forker.py", line 695, in _run
        self.compile_and_execute(example, compiler, test.globs)
      File "/home/sage-patchbot/sage/src/sage/doctest/forker.py", line 1093, in compile_and_execute
        exec(compiled, globs)
      File "<doctest sage.schemes.elliptic_curves.ell_point.EllipticCurvePoint._add_[14]>", line 3, in <module>
        assert R.change_ring(Zmod(d)) == P.change_ring(Zmod(d)) + Q.change_ring(Zmod(d))
      File "/home/sage-patchbot/sage/src/sage/schemes/generic/morphism.py", line 1977, in change_ring
        return S.point(Q, check=check)
      File "/home/sage-patchbot/sage/src/sage/schemes/projective/projective_subscheme.py", line 122, in point
        return self._point(self.point_homset(), v, check=check)
      File "/home/sage-patchbot/sage/src/sage/schemes/projective/projective_point.py", line 203, in __init__
        X.extended_codomain()._check_satisfies_equations(v)
      File "/home/sage-patchbot/sage/src/sage/schemes/generic/algebraic_scheme.py", line 984, in _check_satisfies_equations
        raise TypeError("Coordinates %s do not define a point on %s"%(coords,self))
    TypeError: Coordinates [60, 58, 13] do not define a point on Elliptic Curve defined by y^2 + 40*x*y + 12*y = x^3 + 57*x^2 + 76*x + 55 over Ring of integers modulo 87

roed314 commented 2 years ago

comment:20

Adapting the standard formulas to quotients of Euclidean domains seems reasonable.

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 2 years ago

Changed commit from 2b4509e to 940f8a2

7ed8c4ca-6d56-4ae9-953a-41e42b4ed313 commented 2 years ago

Branch pushed to git repo; I updated commit sha1. This was a forced push. Last 10 new commits:

`ca4ac6b`	`move elliptic-curve point addition to superclass`
`d4a0a40`	`make (P in E) return False when P is defined over an extension field`
`327a902`	`undo ancient hack for elliptic-curve points modulo composites`
`38a8756`	`update documentation to reflect this change`
`0f73f57`	`implement Bosma-Lenstra formulas (corrected version from Best)`
`5fc88be`	`add random test for point addition in ℤ/n`
`5955a1e`	`undo ancient hack for ECM examples`
`9701928`	`add changes to history`
`22a9a7f`	`add doctest for generic addition formulas`
`940f8a2`	`faster isogeny formulas for quotients of Euclidean domains`

yyyyx4 commented 2 years ago

comment:22

I had a shot at adapting the standard formulas to quotients of Euclidean domains.

Performance comparison for a single point addition modulo 2^127-1:

field-specific code:  625 loops, best of 3: 21.8 μs per loop
Euclidean quotients:  625 loops, best of 3: 37.8 μs per loop
generic formulas:     625 loops, best of 3: 123 μs per loop

I think we don't need to add another hack to the EllipticCurve constructor for this because Zmod already has one that does the trick:

sage: R = Zmod(13 * 17, is_field=True)
sage: E = EllipticCurve(R, [1,0])
sage: P = E.lift_x(4)
sage: 123*P
Traceback (most recent call last):
...
ZeroDivisionError: inverse of Mod(34, 221) does not exist

Replying to David Roe:

There's also a doctest failure:

    TypeError: Coordinates [60, 58, 13] do not define a point on Elliptic Curve defined by y^2 + 40*x*y + 12*y = x^3 + 57*x^2 + 76*x + 55 over Ring of integers modulo 87

That's #34417.

yyyyx4 commented 2 years ago

Changed dependencies from #34417 to #34417, #34463

tornaria commented 1 year ago

comment:23

Replying to Lorenz Panny:

I think we don't need to add another hack to the EllipticCurve constructor for this because Zmod already has one that does the trick:
sage: R = Zmod(13 * 17, is_field=True)
sage: E = EllipticCurve(R, [1,0])
sage: P = E.lift_x(4)
sage: 123*P
Traceback (most recent call last):
...
ZeroDivisionError: inverse of Mod(34, 221) does not exist

That's nice, however:

sage: E = EllipticCurve(Zmod(35, is_field=True), [5,1])
sage: E.base_field()
...
AttributeError: 'EllipticCurve_generic_with_category' object has no attribute 'base_field'

in fact

sage: parent(E)
<class 'sage.schemes.elliptic_curves.ell_generic.EllipticCurve_generic_with_category'>

so it would still need "the hack" (broken as of 9.6 and fixed in #34681). Compare to

sage: parent(EllipticCurve(Zmod(7, is_field=True), [5,1]))
<class 'sage.schemes.elliptic_curves.ell_finite_field.EllipticCurve_finite_field_with_category'>

Actually:

sage: Zmod(35, is_field=True).is_field()
False

so, what's the effect of is_field=True?

tornaria commented 1 year ago

comment:24

Also, this is potentially dangerous:

sage: Zmod(35, is_field=True) is Zmod(35, is_field=False)
True

This means the example above works after clearing the cache:

sage: Zmod._cache.clear()
sage: E = EllipticCurve(Zmod(35, is_field=True), [5,1])
sage: E.base_field()
Ring of integers modulo 35

But then one might get a scary error message

sage: Zmod(35).is_field(proof=True)
...
ValueError: THIS SAGE SESSION MIGHT BE SERIOUSLY COMPROMISED!
The order 35 is not prime, but this ring has been put
into the category of fields. This may already have consequences
in other parts of Sage. Either it was a mistake of the user,
or a probabilistic primality test has failed.
In the latter case, please inform the developers.

sagemath / sage

implement elliptic-curve point addition for general rings #33228