isogeny efficiency improvement

[JohnCremona]

Computation of isogenies of prime degree p is expensive when the degree is neither a "genus zero" prime [2,3,5,7,13] or a "hyperelliptic prime" [11, 17, 19, 23, 29, 31, 41, 47, 59, 71] (for these there is special code written). In one situation we can save time, after factoring the degree (p^2-1)/2 division polynomial, if there is exactly one factor of degree (p-1)/2, or one subset of factors whose product has that degree, then the factor of degree (p-1)/2 must be a kernel polynomial. Then we do not need to check consistency, which is very expensive.

The example which led me to this was with p=89 over a quadratic number field, where E.isogeny_class() was taking days. After the change here that goes down to 3 hours. (There are 4 curves in the isogeny class and the code requires factoring the 89-division polynomial of each!) I used a less extreme example for a doctest: a 37-isogeny.

CC: @defeo

Component: elliptic curves

Keywords: isogeny

Author: John Cremona

Branch/Commit: 3d687e5

Reviewer: Jeroen Demeyer

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

[JohnCremona]

Commit: 47ccfd5

[JohnCremona]

New commits:

47ccfd5

#18589 isogeny improvement

[JohnCremona]

Branch: u/cremona/18589

[jdemeyer]

comment:4

Interesting, the hardest part of the computation is checking the consistency, not computing the actual result.

I had a look at what PARI can do for isogenies: ellisogeny computes the isogeny, so it could be used to implement E.isogeny(), but I assume it does no checking.

Minor comment: in the doctest, you don't need

from sage.schemes.elliptic_curves.isogeny_small_degree import isogenies_prime_degree_general

[JohnCremona]

comment:5

Thanks for the comments, I noticed you over at pari-dev and wondered what you were up to.

I don't know how general the pari code is -- here we can handle arbitrary fields (but only separable isogenies). Over number fields the really hard part is determining which prime degrees to test to get the whole class, and I put a lot of work into this (including CM and "potential CM" cases).

Some checking is necessary: I have an example where the 13-division poly factors as 14 degree 6 factors, but only 2 of them are kernel polys! For the rest you need to make a quadratic extension which splits each into two cubics, and then you have to match the factors. The Sage code does this!

[jdemeyer]

comment:6

Replying to @JohnCremona:

I have an example where the 13-division poly factors as 14 degree 6 factors, but only 2 of them are kernel polys!

Is this example in Sage? It would be good to have (if it takes too long, just add # not tested (10 minutes) or whatever)

[jdemeyer]

comment:7

Replying to @JohnCremona:

For the rest you need to make a quadratic extension which splits each into two cubics, and then you have to match the factors.

I'm wondering if you could use resultants instead...

[JohnCremona]

comment:8

Replying to @jdemeyer:

Replying to @JohnCremona:

I have an example where the 13-division poly factors as 14 degree 6 factors, but only 2 of them are kernel polys!

Is this example in Sage? It would be good to have (if it takes too long, just add # not tested (10 minutes) or whatever)

No, but it is Example 3.1.1 on page 28 of [KT2013] which is here: http://wrap.warwick.ac.uk/57568/
The example is over F_3 so is not slow.

[JohnCremona]

comment:9

Replying to @jdemeyer:

Replying to @JohnCremona:

For the rest you need to make a quadratic extension which splits each into two cubics, and then you have to match the factors.

I'm wondering if you could use resultants instead...

I thought that we did use resultants... but looking again at Kimi's code, the relevant function is mult() on line 1960: the effect of this function is to take the polynomial f and the rational function m, under the assumption that f is coprime to the denominator of m, and return the poly whose roots are m(alpha) as alpha runs over the roots of f. If m were a polynomial that would be res_Y(X-m(Y),f(Y)). I can see now that we could stil use that formula with m(Y) replaced by the polynomial p(Y) which satisfies p*d=n (mod f) where m=n/d which would amount to one extended gcd and one resultant. This may be better than the current code since it does not use rational functions.

Can we keep this idea for another ticket?

[jdemeyer]

comment:10

Replying to @JohnCremona:

Can we keep this idea for another ticket?

Sure...

[jdemeyer]

comment:11

In your code, "special case 1" is really a special case of "special case 2". For simplicity, I would therefore keep only special case 2.

[jdemeyer]

comment:12

There is a pre-existing typo: otain instead of obtain.

[jdemeyer]

comment:13

There is a redundant

from sage.misc.all import prod

[JohnCremona]

comment:14

OK, I will fix the redundant imports; I agree that Case 1 is a special case of Case 2. And I am working on the resultant improvement -- so setting to "needs work" is fine. (Without the patch though I would not have been able to compute this: http://beta.lmfdb.org/EllipticCurve/2.2.89.1/81.1/a/ !!

[jdemeyer]

comment:15

Since the constant (l-1)//2 appears in several places in the algorithm, I would prefer to see a variable assigned to it (let's call it D or whatever).

[jdemeyer]

comment:16

Also, can't you avoid some conversions? In

F(f(m(x)).numerator())

the call m(x) is really a conversion which should be done just once, probably best as F(m). Then, I don't think the F() is needed and neither do you need x.

[jdemeyer]

comment:17

In this loop,

        for i in range((l-1)/(2*d)-1):
            g = mult(S[i])
            S.append(g)
            if g in factors:
                factors.remove(g)

when can it happen that g is not in factors?

[jdemeyer]

comment:18

Replying to @JohnCremona:

I thought that we did use resultants... but looking again at Kimi's code, the relevant function is mult() on line 1960: the effect of this function is to take the polynomial f and the rational function m, under the assumption that f is coprime to the denominator of m, and return the poly whose roots are m(alpha) as alpha runs over the roots of f. If m were a polynomial that would be res_Y(X-m(Y),f(Y)). I can see now that we could stil use that formula with m(Y) replaced by the polynomial p(Y) which satisfies p*d=n (mod f) where m=n/d which would amount to one extended gcd and one resultant. This may be better than the current code since it does not use rational functions.

I don't quite get what you're saying here. I think I said "resultants" too quickly before looking at the actual code.

Mathematically, we need to compute

gcd( numerator( f(a/b) ), psi)

where f, a, b, psi are univariate polynomials and only f changes. I guess we know that the denominator of f(a/b) is b^d with d the degree of f. So we need

gcd( b^d f(a/b), psi)

Now suppose c is such that a/b = c mod psi, then we really need

gcd( b^d f(c), psi)

which equals

gcd(f(c), psi)

Now compute this in R[x]/psi(x) and this should be the most efficient way.

(this was written up quickly, I haven't really checked that it's correct)

[jdemeyer]

comment:19

Actually, my idea works, but it's not so efficient since you're working mod psi which is huge.

[jdemeyer]

comment:20

I'm going to set this ticket to needs_work such that you can address the several minor points I mentioned (and answer [comment:17], which I need to understand the algorithm better). Never mind the computation of mult().

[jdemeyer]

comment:21

Replying to @JohnCremona:

No, but it is Example 3.1.1 on page 28 of [KT2013] which is here: http://wrap.warwick.ac.uk/57568/
The example is over F_3 so is not slow.

Wouldn't it be a good testcase for the current code (both over F_3 and F_9)?

[jdemeyer]

comment:39

I think you can simplify

factors = [h for h,e in psi_l.factor() if l2 % h.degree() == 0]

to

factors = [h for h,e in psi_l.factor()]

(in the line below, you are only selecting the degrees you care about anyway)

[JohnCremona]

comment:40

Replying to @jdemeyer:

Replying to @JohnCremona:

Its computation is short, and its factoring takes < 1 hour.

Then why did the computation of the isogenies without this patch take days? What was the bottleneck? I find it hard to believe that the calls to mult() can take so much time.

Me too, but as I tried to explain the script I was running was computing the isogeny class many times. I am rerunning just E.isogeny_class() under 6.7 now and will report back (but this ticket need not wait).

[JohnCremona]

comment:41

Replying to @jdemeyer:

I think you can simplify

factors = [h for h,e in psi_l.factor() if l2 % h.degree() == 0]

to

factors = [h for h,e in psi_l.factor()]

(in the line below, you are only selecting the degrees you care about anyway)

Fine, I will make these simplifications!

[sagetrac-git]

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

2f83c20

#18589: 2 simplifications following review

[sagetrac-git]

Changed commit from ad56369 to 2f83c20

[jdemeyer]

Reviewer: Jeroen Demeyer

[jdemeyer]

comment:43

Replying to @JohnCremona:

Me too, but as I tried to explain the script I was running was computing the isogeny class many times. I am rerunning just E.isogeny_class() under 6.7 now and will report back (but this ticket need not wait).

I agree that this ticket makes sense by itself. It's just that it only does some easy optimizations, it doesn't fundamentally improve the algorithm.

I have been thinking a bit about the computation of mult() this afternoon and I might have an idea to compute mult^-1(f) (that is, find g such that mult(g) == f) much faster than the current mult(). Of course, theoretically, mult() and mult^-1() play the same role.

I'm wondering if it's worth pursuing further and for that I need to know if mult() is ever the bottleneck of the algorithm. (this is of course for a hypothetical follow-up ticket)

[jdemeyer]

comment:44

Replying to @sagetrac-git:

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

2f83c20

#18589: 2 simplifications following review

The check not factors_by_degree only checks that the dict has no keys, but you really need to check that all values are empty lists:

all(factors == [] for factors in factors_by_degree.values())

[JohnCremona]

comment:45

Replying to @jdemeyer:

Replying to @JohnCremona:

Me too, but as I tried to explain the script I was running was computing the isogeny class many times. I am rerunning just E.isogeny_class() under 6.7 now and will report back (but this ticket need not wait).

I agree that this ticket makes sense by itself. It's just that it only does some easy optimizations, it doesn't fundamentally improve the algorithm.

Agreed.

I have been thinking a bit about the computation of mult() this afternoon and I might have an idea to compute mult^-1(f) (that is, find g such that mult(g) == f) much faster than the current mult(). Of course, theoretically, mult() and mult^-1() play the same role.

That is true -- you can use either. You can see in Prop 3.3.1 of the thesis that both are natural! It would be good to improve this.

I'm wondering if it's worth pursuing further and for that I need to know if mult() is ever the bottleneck of the algorithm. (this is of course for a hypothetical follow-up ticket)

[sagetrac-git]

Changed commit from 2f83c20 to 3d687e5

[sagetrac-git]

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

3d687e5

#18589 further one-liner

[JohnCremona]

comment:47

Replying to @jdemeyer:

Replying to @sagetrac-git:

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

2f83c20	#18589: 2 simplifications following review

The check not factors_by_degree only checks that the dict has no keys, but you really need to check that all values are empty lists:

all(factors == [] for factors in factors_by_degree.values())

You are right, I had forgotten that in creating the dict it might have had some empty lists. I should have known better than to do more simplifaction than you recommended!

---

New commits:

3d687e5

#18589 further one-liner

---

New commits:

3d687e5

#18589 further one-liner

[JohnCremona]

comment:48

Any further changes will have to wait until tomorrow....

[jdemeyer]

comment:49

Further work: #18611

[JohnCremona]

comment:50

Many thanks for a fantastic job reviewing this ( excelletn in both mathematical and coding aspects). I will continue the discussion over at #18611.

I suspect that there have been other changes since 6.7 which have sped up my l=89 example, to do with factoring the 89-division polynomial over a quadratic field. As you say, one computation of mult() could not take such a long time even before this patch, yet a single run of E.isogenies_prime_degree(89) is stil running after 14 hours, and all that has to do is run isogenies_prime_degree_general once with l=89.

[jdemeyer]

comment:51

Replying to @JohnCremona:

yet a single run of E.isogenies_prime_degree(89) is stil running after 14 hours

Can you interrupt (CTRL-C) it and look at the traceback to see what it's doing?

[JohnCremona]

comment:52

Replying to @jdemeyer:

Replying to @JohnCremona:

yet a single run of E.isogenies_prime_degree(89) is stil running after 14 hours

Can you interrupt (CTRL-C) it and look at the traceback to see what it's doing?

As expected: it is in the lines

-> 1939         if mult(S[-1]) == f:

and

-> 1928         return gcd(F(f(m(x)).numerator()),psi_l).monic()

and

-> 1588             return a.gcd(b, **kwargs)

and below that

/usr/local/sage/sage-2/src/sage/rings/polynomial/polynomial_number_field.pyx in sage.rings.poly
nomial.polynomial_number_field.Polynomial_absolute_number_field_dense.gcd (build/cythonized/sag
e/rings/polynomial/polynomial_number_field.c:1761)()
    203         h1 = self._pari_with_name('x')
    204         h2 = other._pari_with_name('x')
--> 205         g = h1.gcd(h2)

and at the bottom, pari/gen.pyx

[jdemeyer]

comment:53

Thanks for the info.

Interesting, it looks like the gcd() computation in mult() takes so much time. I would assume that some kind of coefficient explosion occurs.

[JohnCremona]

comment:54

Replying to @jdemeyer:

Thanks for the info.

Interesting, it looks like the gcd() computation in mult() takes so much time. I would assume that some kind of coefficient explosion occurs.

I think that the real point here is #18461 which re-implemented univariate polynomial gcd! So it is possible that without any of the changes on this ticket or #18611 we would have seen anoticeable improvement. Anyway, between these two tickets we certainly have more efficient code so the exercise was worth carrying out. Thanks!

[jdemeyer]

comment:55

Replying to @JohnCremona:

Replying to @jdemeyer:

Thanks for the info.

Interesting, it looks like the gcd() computation in mult() takes so much time. I would assume that some kind of coefficient explosion occurs.

I think that the real point here is #18461 which re-implemented univariate polynomial gcd!

From the traceback, you see that PARI is used to compute this gcd, so the generic Sage code doesn't matter.

[vbraun]

Changed branch from u/cremona/18589 to 3d687e5