Add BMSS algorithm for isogenies of elliptic curves

[defeo]

Add the algorithm by Bostan, Morain, Salvy and Schost for normalized isogenies in characteristic 0 or large enough.

References:

[BMSS08] Alin Bostan, François Morain, Bruno Salvy, and Éric Schost, Fast algorithms for computing isogenies between elliptic curves, Mathematics of Computation 77 (2008), 1755–1778.

Apply: attachment: trac_11095_isogenies_char_0_v3.patch​

CC: @categorie @jpflori @pjbruin

Component: elliptic curves

Keywords: isogenies ecc2011

Work Issues: merge conflicts

Author: Luca De Feo

Branch/Commit: public/ticket/11095 @ 59417df

Reviewer: John Cremona, Marco Streng

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

[defeo]

Description changed:

--- 
+++ 
@@ -1,4 +1,4 @@
-Add the algorithm by Bostain, Morain, Salvy and Schost for normalized cyclic isogenies in characteristic 0 or large enough.
+Add the algorithm by Bostain, Morain, Salvy and Schost for normalized isogenies in characteristic 0 or large enough.

 References:

[defeo]

comment:2

Besides implementing BMSS, I cleaned a bit ell_curve_isogeny, moving the algorithms that do not work in any characteristic to a new module ell_isogeny_char_zero. In my view, only the EllipticCurveIsogeny class and Vélu/Kohel formulae should remain in ell_curve_isogeny` in the end.

BMSS runs significantly faster than Stark's, and works for slightly larger characteristics, I thus set it as the default choice. I haven't produced any serious benchmarks, but I will if the reviewers ask me to.

[JohnCremona]

comment:3

I like this a lot, and it is a good idea to start moving stuff out of the huge isogenies file. I am planning to do the same for the isogenies of small degree code --currently the genus 0 prime degrees 2,3,5,7,13 and the others which can occur over Q. Kimi Tsukazaki and I are working on extending that to all the ell for which X_0(ell)^+ has genus 0.

I tested this against 4.7.rc1 and had some doctest failures. I expect these will be easy to fix.

sage -t  sage/schemes/elliptic_curves/ell_isogeny_char_zero.py
// ** nInitExp failed: using Z/2^2
**********************************************************************
File "/home/jec/sage-4.7.rc1/devel/sage-main/sage/schemes/elliptic_curves/ell_isogeny_char_zero.py", line 165:
    sage: isogeny_Stark(E1pr, E2pr, 11)
Exception raised:
    Traceback (most recent call last):
      File "/home/jec/sage-current/local/bin/ncadoctest.py", line 1231, in run_one_test
        self.run_one_example(test, example, filename, compileflags)
      File "/home/jec/sage-current/local/bin/sagedoctest.py", line 38, in run_one_example
        OrigDocTestRunner.run_one_example(self, test, example, filename, compileflags)
      File "/home/jec/sage-current/local/bin/ncadoctest.py", line 1172, in run_one_example
        compileflags, 1) in test.globs
      File "<doctest __main__.example_2[7]>", line 1, in <module>
        isogeny_Stark(E1pr, E2pr, Integer(11))###line 165:
    sage: isogeny_Stark(E1pr, E2pr, 11)
    NameError: name 'E1pr' is not defined
**********************************************************************
File "/home/jec/sage-4.7.rc1/devel/sage-main/sage/schemes/elliptic_curves/ell_isogeny_char_zero.py", line 334:
    sage: kernel = isogeny_BMSS(E1pr, E2pr, 9)
Exception raised:
    Traceback (most recent call last):
      File "/home/jec/sage-current/local/bin/ncadoctest.py", line 1231, in run_one_test
        self.run_one_example(test, example, filename, compileflags)
      File "/home/jec/sage-current/local/bin/sagedoctest.py", line 38, in run_one_example
        OrigDocTestRunner.run_one_example(self, test, example, filename, compileflags)
      File "/home/jec/sage-current/local/bin/ncadoctest.py", line 1172, in run_one_example
        compileflags, 1) in test.globs
      File "<doctest __main__.example_3[23]>", line 1, in <module>
        kernel = isogeny_BMSS(E1pr, E2pr, Integer(9))###line 334:
    sage: kernel = isogeny_BMSS(E1pr, E2pr, 9)
    NameError: name 'E1pr' is not defined
**********************************************************************
File "/home/jec/sage-4.7.rc1/devel/sage-main/sage/schemes/elliptic_curves/ell_isogeny_char_zero.py", line 335:
    sage: kernel
Expected:
    x^8 + 36*x^7 + 95*x^6 + 84*x^4 + 16*x^3 + 70*x^2 + 12*x + 78
Got:
    x^4 + 14*x^3 + x^2 + 34*x + 21
**********************************************************************
File "/home/jec/sage-4.7.rc1/devel/sage-main/sage/schemes/elliptic_curves/ell_isogeny_char_zero.py", line 337:
    sage: kernel.is_square()
Expected:
    False
Got:
    True
**********************************************************************
File "/home/jec/sage-4.7.rc1/devel/sage-main/sage/schemes/elliptic_curves/ell_isogeny_char_zero.py", line 339:
    sage: isogeny_kernel(E1pr, E2pr, 9, algorithm="bmss")
Expected:
    Traceback (most recent call last):
    ...
    ValueError: The two curves are not linked by a normalized isogeny of degree 9
Got:
    Traceback (most recent call last):
      File "/home/jec/sage-current/local/bin/ncadoctest.py", line 1231, in run_one_test
        self.run_one_example(test, example, filename, compileflags)
      File "/home/jec/sage-current/local/bin/sagedoctest.py", line 38, in run_one_example
        OrigDocTestRunner.run_one_example(self, test, example, filename, compileflags)
      File "/home/jec/sage-current/local/bin/ncadoctest.py", line 1172, in run_one_example
        compileflags, 1) in test.globs
      File "<doctest __main__.example_3[26]>", line 1, in <module>
        isogeny_kernel(E1pr, E2pr, Integer(9), algorithm="bmss")###line 339:
    sage: isogeny_kernel(E1pr, E2pr, 9, algorithm="bmss")
    NameError: name 'E1pr' is not defined
**********************************************************************
2 items had failures:
   1 of  22 in __main__.example_2
   4 of  28 in __main__.example_3
***Test Failed*** 5 failures.
For whitespace errors, see the file /home/jec/.sage//tmp/.doctest_ell_isogeny_char_zero.py
     [3.8 s]

[JohnCremona]

Reviewer: John Cremona

[defeo]

comment:4

Ooops! I had made some last minute change to the doctest and must have tested the wrong file afterwards. The new patch passes all tests, including -testall -long

[defeo]

comment:6

I realized that my code for BMSS fails when the isogeny has strictly positive valuation at x, which happens for example when the origin is a point of order 2 and the isogeny has odd degree as in the code below

sage: from sage.schemes.elliptic_curves.ell_isogeny_char_zero import *
sage: E = EllipticCurve([-1,0])
sage: E2 = EllipticCurve([-3^4,0])
sage: E.isogeny(kernel=None, codomain=E2, degree=9)
Traceback (most recent call last)
...
ValueError: The two curves are not linked by a normalized isogeny of degree 9
sage: isogeny_Stark(E,E2,9)
x^8 - 4*x^6 + 10/3*x^4 + 4/3*x^2 + 1/9
sage: isogeny_BMSS(E,E2,9)
x^7 - 4*x^5 + 10/3*x^3 + 4/3*x

The fix is easy, I hope I'll have time to send a revised patch next week.

[defeo]

comment:7

Attachment: trac_11095_isogenies_char_0.patch.gz

Finally, here's the easy fix! I hope it will be ok now.

[mstreng]

comment:8

On sage 4.7.1.alpha3 with this patch:

sage -t -long devel/sage/sage/interfaces/ecm.py
*** *** Error: TIMED OUT! PROCESS KILLED! *** ***

         [1800.5 s]

[mstreng]

Changed reviewer from John Cremona to John Cremona, Marco Streng

[defeo]

comment:9

Thanks for testing, Marco. This is wierd: I can't seem to reproduce the error; besides, I don't see how the ecm.py module interacts with the elliptic_curve package.

I tested on sage 4.7.1.alpha4 with and without the patch and I have no failure in any case. Does anyone understand what's going on? Can anyone confirm the bug?

[mstreng]

comment:10

Replying to @defeo:

Thanks for testing, Marco. This is wierd: I can't seem to reproduce the error; besides, I don't see how the ecm.py module interacts with the elliptic_curve package.

I tested on sage 4.7.1.alpha4 with and without the patch and I have no failure in any case. Does anyone understand what's going on? Can anyone confirm the bug?

I'll see if I can reproduce it, and will let you know tomorrow.

[mstreng]

Changed reviewer from John Cremona, Marco Streng to John Cremona

[mstreng]

comment:11

Replying to @mstreng:

I'll see if I can reproduce it, and will let you know tomorrow.

The problem vanished. Sorry for the trouble.

[defeo]

comment:12

Added depracation warnings, as suggested by Marco. Apply second patch only (is this going to confuse the patchbot?).

[mstreng]

comment:13

1. I'm confused by the ValueErrors raised by isogeny_kernel.

ValueError, "The two curves are not linked by a normalized isogeny of degree %s"%ell

First of all, the comment "or a non-sense polynomial if no such isogeny exists" (line 110) is inconsistent with the documentation of isogeny_BMSS. Which of the two is correct?

Second, the test "ker_poly.degree() != ell-1" (line 113) implicitly assumes that there is only one possible degree for a normalized isogeny. It would help future developers to have a comment explaining that this is always the case (it is, isn't it?)

Third, why is the test "not check" there in the Stark case (line 100)? Could isogeny_Stark return non-sense polynomials just like you say isogeny_BMSS does? Does a non-sense polynomial always imply that no normalized isogeny exist?

2. The word "separable" seems to be missing from the ValueErrors

3. I noticed that (with the default options) the function isogeny_kernel will first test if algorithm equals "Stark", "stark", "Starks", or "starks", before detecting that it actually equals "BMSS". It may speed things up by epsilon to put the most-used case first. (this is not the reason for "needs work")

ps. The patch bot seems to have understood that you only want the second patch.

[JohnCremona]

comment:14

Over a finite field theer are always lots of isogeny degrees between any two isogenous curves. For ordinary curves the set of degrees is the set of all positive integers represented by a ceratin positive definite quadratic form; for supersingular curves every degree except a finite number are possible (see Kohel's thesis).

[mstreng]

Changed reviewer from John Cremona to John Cremona, Marco Streng

[mstreng]

comment:15

Replying to @JohnCremona:

Over a finite field theer are always lots of isogeny degrees between any two isogenous curves.

The word "normalized" is important here. That should limit the possibilities.

[JohnCremona]

comment:16

Replying to @mstreng:

Replying to @JohnCremona:

Over a finite field theer are always lots of isogeny degrees between any two isogenous curves.

The word "normalized" is important here. That should limit the possibilities.

Not if you do not mind replacing the codomain with an isomorphic curve, since composing with a scaling map ([u,r,s,t]=[u,0,0,0]) multiplies the constant by u (or maybe u^-1) so u can always be chosen to make the isogeny normalised!

[mstreng]

comment:17

Replying to @JohnCremona:

Not if you do not mind replacing the codomain with an isomorphic curve

I don't mind, but the code in this patch does mind.

[defeo]

comment:18

Replying to @mstreng:

First of all, the comment "or a non-sense polynomial if no such isogeny exists" (line 110) is inconsistent with the documentation of isogeny_BMSS. Which of the two is correct?

I fixed the docs to be more explicit on this, hope it is clear now.

Second, the test "ker_poly.degree() != ell-1" (line 113) implicitly assumes that there is only one possible degree for a normalized isogeny. It would help future developers to have a comment explaining that this is always the case (it is, isn't it?)

Marco is right: the "normalized" requirement is very strong. Indeed if two isogenies are normalized for the same models, then they are solution to the same differential equation with the same initial conditions. In characteristic 0, this means that their power series expansion at infinity is the same, which I think should be enough to prove that there is only one such isogeny.

In positive characteristic, this means that the power series are equal at least up to the (p-1)/2-th term (which is the limit up to which the BMSS algorithm can compute the solution). I am less sure about unicity in this case, but I am inclined to think that only one isogeny of degree less than p/4 will satisfy the condition (this is true for multiplication-by-m maps, at least).

As an example, take the multiplication maps by 1 and by p-1 on a curve E: they are both normalized for the models (E,E), which implies that the expansion at infinity of the multiplication-by-(p-1) map is x + O(x^((p-1)/2)). Of course the BMSS algorithm will refuse to compute the second map, because it has degree (p-1)^2 > p/4.

As a final remark, if we know the right (u,r,s,t) we can normalize the models, as John suggests. But, except for the case of scalar multiplication maps, computing such an isomorphisms is very expensive. I'm working on this in another patch.

Third, why is the test "not check" there in the Stark case (line 100)? Could isogeny_Stark return non-sense polynomials just like you say isogeny_BMSS does? Does a non-sense polynomial always imply that no normalized isogeny exist?

Maybe not, but these algorithms all tend to have some special cases which were not considered in the literature and which yield false answers. So why not test as long as the test comes for free (we have to compute the square root anyway).

2. The word "separable" seems to be missing from the ValueErrors

The constraints on the degree imply that if such an error is returned, then p does not divide ell (otherwise a ZeroDivisionError is returned instead), thus no unseparable isogeny of degree ell exists either.

3. I noticed that (with the default options) the function isogeny_kernel will first test if algorithm equals "Stark", "stark", "Starks", or "starks", before detecting that it actually equals "BMSS". It may speed things up by epsilon to put the most-used case first. (this is not the reason for "needs work")

done

ps. The patch bot seems to have understood that you only want the second patch.

Yes, but it also used to mention a failure that I was not able to reproduce. It seems to be down now, anyway.

[defeo]

Attachment: trac_11095_isogenies_char_0.2.patch.gz

[zimmermann6]

comment:19

what are the work issues for that ticket?

Paul

[zimmermann6]

Changed keywords from isogenies to isogenies ecc2011

[zimmermann6]

Description changed:

--- 
+++ 
@@ -1,4 +1,4 @@
-Add the algorithm by Bostain, Morain, Salvy and Schost for normalized isogenies in characteristic 0 or large enough.
+Add the algorithm by Bostan, Morain, Salvy and Schost for normalized isogenies in characteristic 0 or large enough.

 References:

[fchapoton]

comment:21

here is a patch rebased on 5.12.beta4

for the bot :

apply only trac_11095_isogenies_char_0_v3.patch​

[fchapoton]

Description changed:

--- 
+++ 
@@ -4,3 +4,4 @@

 [BMSS08] Alin Bostan, François Morain, Bruno Salvy, and Éric Schost, Fast algorithms for computing isogenies between elliptic curves, Mathematics of Computation 77 (2008), 1755–1778.

+Apply: [attachment: trac_11095_isogenies_char_0_v3.patch​](https://github.com/sagemath/sage/files/ticket11095/0e2630fbe9449c551f50bca100e732dd.gz)

[fchapoton]

comment:22

apply only trac_11095_isogenies_char_0_v3.patch

[fchapoton]

Attachment: trac_11095_isogenies_char_0_v3.patch.gz

[fchapoton]

comment:23

apply only trac_11095_isogenies_char_0_v3.patch

[JohnCremona]

comment:24

This is going to conflict with #13615 -- can we agree on an order to do them in?

[defeo]

comment:25

Oh, definitely #13615 first.

This ticket has been bitrotting for years, although I've never abandoned it and I still use the code from time to time. The new git workflow might be the occasion to go back to it... maybe. In the meantime, I'll put the milestone to sage-pending.

[pjbruin]

comment:27

I converted the latest patch into a Git branch; there were some conflicts that I hopefully resolved correctly, and all doctests pass.

Would it make sense to call the new module isogeny_char_zero instead of ell_isogeny_char_zero? There is also isogeny_small_degree; it would be nice if we could avoid having modules with three different prefixes (isogeny*, ell_isogeny* and ell_curve_isogeny*)...

[pjbruin]

Branch: u/pbruin/11095-isogeny_char_zero

[pjbruin]

Commit: 12d8a9c

[JohnCremona]

comment:28

Good point about naming of modules. there are 14 with names starting ell_ in the directory called ellipticcurves and the ell prefix could be removed from all of them! (There is also lseries_ell). But doing so would cause a lot of problems for any other branches affecting any of those files...

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

8cae4a7

Trac 11095: rename ell_isogeny_char_zero to isogeny_char_zero

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 12d8a9c to 8cae4a7

[pjbruin]

comment:30

Since John and I are for it and nobody seems to be against it, I removed the ell_ prefix from isogeny_char_zero. (For existing files that would of course be too drastic to do here.)

[defeo]

comment:31

Thanks Peter for resurrecting this old ticket of mine. I have no objections to the file renaming.

Back in 2011 Marco Streng and I bypassed Trac and discussed this ticket via mail. I feel I must make an update to explain why the ticket was set to needs work back then.

Essentially, Marco was unhappy about the documentation for two reasons.

1. It wasn't clear whether Stark's algorithm could fail silently. I just had a look at it. It can. I'm commiting the necessary modifications to the docs.

2. The docstring talks about THE normalized isogeny of degree at most degree, but it wasn't immediately clear that this isogeny should be unique. Marco himself gave the proof of unicity. This is explained in the note in the docstring of isogeny_BMSS.

I can't find any other issue in my mail logs, so it is good to go for me.

[defeo]

Changed branch from u/pbruin/11095-isogeny_char_zero to u/defeo/ticket/11095

[pjbruin]

Changed commit from 8cae4a7 to 3c978f8

[pjbruin]

comment:33

OK. I cannot review it right now, so I just changed the milestone to "undiscourage" people from reviewing it.

---

New commits:

3c978f8

Modifed docstrings and error messages

[fchapoton]

Changed commit from 3c978f8 to 8257108

[fchapoton]

Changed branch from u/defeo/ticket/11095 to public/ticket/11095

[fchapoton]

comment:35

I have made a few minor changes, including:

• put new raise statements into python3 syntax
• alignement changes in doc

---

New commits:

9f71709	Merge branch 'u/defeo/ticket/11095' of ssh://trac.sagemath.org:22/sage into 11095
8257108	trac #11095 minor details (doc formatting and python3 raise syntax)

[pjbruin]

comment:36

I've been experimenting a bit and it seems to work well. Would it be possible to make EllipticCurve.isogeny() and EllipticCurveIsogeny accept algorithm='bmss'? In the following example I wanted to compute a 20-isogeny between two elliptic curves over a quadratic field (knowing a priori that such an isogeny exists). Trying E.isogeny() gives a NotImplementedError and it seems that I have to import isogeny_kernel:

sage: K.<a> = NumberField(x^2 - 17*x - 120)
sage: E = EllipticCurve([0, 1, 0, -89423170640*a + 1999987182576, 87486176701141216*a - 1956665490639249900])
sage: E_isog = EllipticCurve([0, 0, 0, 1/3*(6255584200920576*a - 139908796923911440), 1/27*(-5210907833850999276681216*a + 116544166379886781944583040)])
sage: E.isogeny(kernel=None, codomain=E_isog, degree=20)
...
NotImplementedError: For basic Kohel's algorithm, if the kernel degree is even then the kernel must be contained in the two torsion.
sage: from sage.schemes.elliptic_curves.isogeny_char_zero import isogeny_kernel 
sage: isogeny_kernel(E, E_isog, 20)
x^10 + (-2152280*a + 48136659)*x^9 + (-148051244547864/5*a + 3311228950404238/5)*x^8 + (-129467694651804753040*a + 2895599965960126731024)*x^7 + (-131438756167338709694710304*a + 2939683593714796730731625088)*x^6 + (196151384137228169247793952562368*a - 4387008996787386307912516668742368)*x^5 + (211127593344837320952211060140401067136*a - 4721958274971357551356075375634080673472)*x^4 + (-676089273725371340027412364537916253836377856/5*a + 15121023690506428142245429629695675550575584512/5)*x^3 + (-58033791954670621184885097803503676816467494405632*a + 1297950399599077493253508317741911096974395766973184)*x^2 + (26614381917191620756309912883562400685500952121788824576*a - 595241952679626281929246822605533328842111805768309632256)*x + 665314371570672918992258944764049661423822960574864980480000*a - 14880038428536031844363469294864157021058877011246176591280640

By the way, in this example Stark's algorithm fails; I haven't checked exactly why.

sage: isogeny_kernel(E, E_isog, 20, algorithm='stark')
...
RuntimeError: Stark's algorithm returned an unexpected result. Please report this bug.

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

Changed commit from 8257108 to c7f35aa

[7ed8c4ca-6d56-4ae9-953a-41e42b4ed313]

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

c7f35aa

Merge branch 'develop' into ticket/11095-isogeny_char_zero

[pjbruin]

comment:39

The merge conflicts appear to result from #11327.