Opened 8 years ago

Closed 8 years ago

# isogeny efficiency improvement

Reported by: Owned by: cremona major sage-6.8 elliptic curves isogeny defeo John Cremona Jeroen Demeyer N/A 3d687e5 3d687e5225f67808eb6c5af5fbf4cb93f2000c62

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.

### comment:1 Changed 8 years ago by cremona

Description: modified (diff)

### comment:2 Changed 8 years ago by cremona

Branch: → u/cremona/18589 → 47ccfd587402b953c612fcd3cddaa541a6847bd3 modified (diff) new → needs_review

New commits:

 ​47ccfd5 `#18589 isogeny improvement`

### comment:4 Changed 8 years ago by jdemeyer

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
```

### comment:5 follow-ups:  6  7 Changed 8 years ago by cremona

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!

### comment:6 in reply to:  5 ; follow-up:  8 Changed 8 years ago by jdemeyer

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)

### comment:7 in reply to:  5 ; follow-up:  9 Changed 8 years ago by jdemeyer

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...

### comment:8 in reply to:  6 ; follow-ups:  15  21 Changed 8 years ago by cremona

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.

### comment:9 in reply to:  7 ; follow-ups:  10  18 Changed 8 years ago by cremona

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?

### comment:10 in reply to:  9 Changed 8 years ago by jdemeyer

Can we keep this idea for another ticket?

Sure...

### comment:11 follow-up:  24 Changed 8 years ago by jdemeyer

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.

### comment:12 follow-up:  25 Changed 8 years ago by jdemeyer

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

### comment:13 follow-up:  26 Changed 8 years ago by jdemeyer

There is a redundant

```from sage.misc.all import prod
```

### comment:14 follow-up:  33 Changed 8 years ago by cremona

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/ !!

### comment:15 in reply to:  8 ; follow-up:  27 Changed 8 years ago by jdemeyer

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).

### comment:16 Changed 8 years ago by jdemeyer

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`.

### comment:17 follow-up:  28 Changed 8 years ago by jdemeyer

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`?

### comment:18 in reply to:  9 Changed 8 years ago by jdemeyer

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)

Last edited 8 years ago by jdemeyer (previous) (diff)

### comment:19 Changed 8 years ago by jdemeyer

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

### comment:20 Changed 8 years ago by jdemeyer

Status: needs_review → needs_work

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

Last edited 8 years ago by jdemeyer (previous) (diff)

### comment:21 in reply to:  8 Changed 8 years ago by jdemeyer

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)?

### comment:22 follow-up:  29 Changed 8 years ago by jdemeyer

Since you write `Every kernel polynomial is a product of irreducible factors of the division polynomial of the same degree`, could special case 2 apply by just considering the polynomials of a given degree?

### comment:23 Changed 8 years ago by cremona

Thanks for all the comments -- I was busy with other things so have not replied to them indiviually yet, but perhaps it would be clearer if I now did so!

### comment:24 in reply to:  11 Changed 8 years ago by cremona

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.

Agreed -- done.

### comment:25 in reply to:  12 Changed 8 years ago by cremona

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

Done.

### comment:26 in reply to:  13 Changed 8 years ago by cremona

There is a redundant

```from sage.misc.all import prod
```

Done.

### comment:27 in reply to:  15 Changed 8 years ago by cremona

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).

Done -- I called it l2.

### comment:28 in reply to:  17 ; follow-up:  30 Changed 8 years ago by cremona

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`?

Very good question! I think it always happens, since the operation f --> mult(f) applies the multiplication-by-a map to the roots, so must take f to another irreducible factor of the same degree; the fact being checked is that the resulting cycle has the correct length, which is (l-1)/2 divided by deg(f), so that the product of the f's in the cycle gives the kernel polynomial. [I will delete the check.] I tried removing the check and got run-time errors, so my mathematics must be wrong...

Last edited 8 years ago by cremona (previous) (diff)

### comment:29 in reply to:  22 Changed 8 years ago by cremona

Since you write `Every kernel polynomial is a product of irreducible factors of the division polynomial of the same degree`, could special case 2 apply by just considering the polynomials of a given degree?

Yes: if for some degree d dividing l2=(l-1)/2 there are exactly (l-1)/2d factors then their product must be a kernel poly.

Some more theory: each subgroup has l2 x-coordinates and the rational function m(x) permutes these in a single cycle (by definition of "semi-primitive root" as a class mod l which generates the multiplcative group modulo <-1>). A kernel poly is a poly of degree l2 whose set of roots is invariant under this map. Since the map is defined over the ground field, each of the x-coordinates in any subgroup has the same degree d over the base, so d is a divisor of l2 and the kernel poly is a product of l2/d polys of degree d. On the other hand, an irreducible factor f of degree d (where d didvides l2) need not be a factor of a kernel poly, as the F_3 example with l=13 shows. If that is the case then the orbit of f under mult() will be longer than l2/d and we will detect that after l2/d steps we have not returned to the start. In this circumstance we will have removed l2/d factors from the list but this is only *part of* an orbit, so later on when we start with one of the remaining factors in that orbit, we may well iterate onto a factor which is already removed.

There, I have proved that g need not be in factors! And more -- if we find a g which is not in factors then we have hit a part-orbit already seen so can break.

### comment:30 in reply to:  28 ; follow-up:  31 Changed 8 years ago by jdemeyer

Very good question! I think it always happens, since the operation f --> mult(f) applies the multiplication-by-a map to the roots, so must take f to another irreducible factor of the same degree; the fact being checked is that the resulting cycle has the correct length, which is (l-1)/2 divided by deg(f), so that the product of the f's in the cycle gives the kernel polynomial. [I will delete the check.] I tried removing the check and got run-time errors, so my mathematics must be wrong...

I also used the same thinking as you, and I also don't understand the problem...

### comment:31 in reply to:  30 Changed 8 years ago by cremona

Very good question! I think it always happens, since the operation f --> mult(f) applies the multiplication-by-a map to the roots, so must take f to another irreducible factor of the same degree; the fact being checked is that the resulting cycle has the correct length, which is (l-1)/2 divided by deg(f), so that the product of the f's in the cycle gives the kernel polynomial. [I will delete the check.] I tried removing the check and got run-time errors, so my mathematics must be wrong...

I also used the same thinking as you, and I also don't understand the problem...

See comment 29!

### comment:32 Changed 8 years ago by cremona

OK, I have done quite a rewrite which makes the logic much clearer (I think) and also treats all the special cases together. New commit soon...

### comment:33 in reply to:  14 ; follow-up:  34 Changed 8 years ago by jdemeyer

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/ !!

I'm actually curious why you get a huge speed-up here. I can see that the patch here improves things, but to go from "not being able to compute" to "being able to compute" surprises me. Doesn't the computation and factoring of the division polynomial dominate the whole computation anyway?

Last edited 8 years ago by jdemeyer (previous) (diff)

### comment:34 in reply to:  33 ; follow-up:  37 Changed 8 years ago by cremona

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/ !!

I'm actually curious why you get a huge speed-up here. I can see that the patch here improves things, but to go from "not being able to compute" to "being able to compute" surprises me. Doesn't the computation and factoring of the division polynomial dominate the whole computation anyway?

Its computation is short, and its factoring takes < 1 hour. In the code I had, computing the complete isogeny class of 4 curves, I had to do that 4 times. (I am working on not having to recompute the kernel polys for other curves in the isogeny class, but that will definitely be another ticket). I also had some lines like E1.is_isogenous(E2,proof=True) which triggered computation of the whole isogeny class of E1. (With proof=False it just checks traces of Frobenius but I also had curves which were not isogenous but whose traces agreed at a lot of primes, so was playing safe).

### comment:35 Changed 8 years ago by git

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

 ​ad56369 `#18589 simplified code + further improvements`

### comment:36 Changed 8 years ago by cremona

Status: needs_work → needs_review

Please take a look at this. the code is (I hope) simpler to understand, deals with factors of each degree separately, and has more explanatory comments. I experimented with a couple of other version of mult(), one of which is there but commented out, but they were no faster. I also added a new example as discussed above.

### comment:37 in reply to:  34 ; follow-up:  40 Changed 8 years ago by jdemeyer

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.

### comment:38 Changed 8 years ago by jdemeyer

In the current implementation, you are computing

```    from sage.rings.arith import gcd
a = _least_semi_primitive(l)
m = E.multiplication_by_m(a, x_only=True)
F = psi_l.parent()
x = F.gen()
d = F(m.denominator())
n = F(m.numerator())
```

even when `factors_by_degree` is empty.

An easy optimization would be

```if all(factors == [] for factors in factors_by_degree):
return ...
```

just before that block.

### comment:39 follow-up:  41 Changed 8 years ago by 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)

### comment:40 in reply to:  37 ; follow-up:  43 Changed 8 years ago by cremona

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).

### comment:41 in reply to:  39 Changed 8 years ago by cremona

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!

### comment:42 follow-up:  44 Changed 8 years ago by git

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

 ​2f83c20 `#18589: 2 simplifications following review`

### comment:43 in reply to:  40 ; follow-up:  45 Changed 8 years ago by jdemeyer

Reviewers: → Jeroen Demeyer

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)

### comment:44 in reply to:  42 ; follow-up:  47 Changed 8 years ago by jdemeyer

Status: needs_review → needs_work

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())
```

### comment:45 in reply to:  43 Changed 8 years ago by cremona

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)

### comment:46 Changed 8 years ago by git

Commit: 2f83c20a3753b9141c635ab26ed88015a9fed499 → 3d687e5225f67808eb6c5af5fbf4cb93f2000c62

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

 ​3d687e5 `#18589 further one-liner`

### comment:47 in reply to:  44 Changed 8 years ago by cremona

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`

### comment:48 Changed 8 years ago by cremona

Status: needs_work → needs_review

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

### comment:49 Changed 8 years ago by jdemeyer

Status: needs_review → positive_review

Further work: #18611

### comment:50 follow-up:  51 Changed 8 years ago by cremona

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.

### comment:51 in reply to:  50 ; follow-up:  52 Changed 8 years ago by jdemeyer

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?

### comment:52 in reply to:  51 Changed 8 years ago by cremona

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

### comment:53 follow-up:  54 Changed 8 years ago by 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.

### comment:54 in reply to:  53 ; follow-up:  55 Changed 8 years ago by cremona

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 http://trac.sagemath.org/ticket/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!

### comment:55 in reply to:  54 Changed 8 years ago by jdemeyer

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