Opened 7 years ago

Closed 7 years ago

#18589 closed enhancement (fixed)

isogeny efficiency improvement

Reported by: cremona Owned by:
Priority: major Milestone: sage-6.8
Component: elliptic curves Keywords: isogeny
Cc: defeo Merged in:
Authors: John Cremona Reviewers: Jeroen Demeyer
Report Upstream: N/A Work issues:
Branch: 3d687e5 (Commits, GitHub, GitLab) Commit: 3d687e5225f67808eb6c5af5fbf4cb93f2000c62
Dependencies: Stopgaps:

Status badges

Description (last modified by cremona)

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.

Change History (56)

comment:1 Changed 7 years ago by cremona

  • Description modified (diff)

comment:2 Changed 7 years ago by cremona

  • Branch set to u/cremona/18589
  • Commit set to 47ccfd587402b953c612fcd3cddaa541a6847bd3
  • Description modified (diff)
  • Status changed from new to needs_review

New commits:

47ccfd5#18589 isogeny improvement

comment:3 Changed 7 years ago by defeo

  • Cc defeo added

comment:4 Changed 7 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: Changed 7 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: Changed 7 years ago by jdemeyer

Replying to 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)

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

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

comment:8 in reply to: ↑ 6 ; follow-ups: Changed 7 years ago by cremona

Replying to jdemeyer:

Replying to 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: Changed 7 years ago by cremona

Replying to jdemeyer:

Replying to 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 7 years ago by jdemeyer

Replying to cremona:

Can we keep this idea for another ticket?

Sure...

comment:11 follow-up: Changed 7 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: Changed 7 years ago by jdemeyer

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

comment:13 follow-up: Changed 7 years ago by jdemeyer

There is a redundant

from sage.misc.all import prod

comment:14 follow-up: Changed 7 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: Changed 7 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 7 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: Changed 7 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 7 years ago by jdemeyer

Replying to cremona:

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 7 years ago by jdemeyer (previous) (diff)

comment:19 Changed 7 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 7 years ago by jdemeyer

  • Status changed from needs_review to 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 7 years ago by jdemeyer (previous) (diff)

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

Replying to cremona:

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: Changed 7 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 7 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 7 years ago by cremona

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

Agreed -- done.

comment:25 in reply to: ↑ 12 Changed 7 years ago by cremona

Replying to jdemeyer:

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

Done.

comment:26 in reply to: ↑ 13 Changed 7 years ago by cremona

Replying to jdemeyer:

There is a redundant

from sage.misc.all import prod

Done.

comment:27 in reply to: ↑ 15 Changed 7 years ago by cremona

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

Done -- I called it l2.

comment:28 in reply to: ↑ 17 ; follow-up: Changed 7 years ago by cremona

Replying to 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?

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 7 years ago by cremona (previous) (diff)

comment:29 in reply to: ↑ 22 Changed 7 years ago by cremona

Replying to 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?

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: Changed 7 years ago by jdemeyer

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

comment:31 in reply to: ↑ 30 Changed 7 years ago by cremona

Replying to jdemeyer:

Replying to 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 7 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: Changed 7 years ago by jdemeyer

Replying to 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?

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

comment:34 in reply to: ↑ 33 ; follow-up: Changed 7 years ago by cremona

Replying to jdemeyer:

Replying to 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 7 years ago by git

  • Commit changed from 47ccfd587402b953c612fcd3cddaa541a6847bd3 to ad56369f5def8060575b6567342dedb28a447293

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

ad56369#18589 simplified code + further improvements

comment:36 Changed 7 years ago by cremona

  • Status changed from needs_work to 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: Changed 7 years ago by jdemeyer

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

comment:38 Changed 7 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: Changed 7 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: Changed 7 years ago by cremona

Replying to jdemeyer:

Replying to 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 7 years ago by cremona

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!

comment:42 follow-up: Changed 7 years ago by git

  • Commit changed from ad56369f5def8060575b6567342dedb28a447293 to 2f83c20a3753b9141c635ab26ed88015a9fed499

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: Changed 7 years ago by jdemeyer

  • Reviewers set to Jeroen Demeyer

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

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: Changed 7 years ago by jdemeyer

  • Status changed from needs_review to needs_work

Replying to 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())

comment:45 in reply to: ↑ 43 Changed 7 years ago by cremona

Replying to jdemeyer:

Replying to 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 7 years ago by git

  • Commit changed from 2f83c20a3753b9141c635ab26ed88015a9fed499 to 3d687e5225f67808eb6c5af5fbf4cb93f2000c62

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

3d687e5#18589 further one-liner

comment:47 in reply to: ↑ 44 Changed 7 years ago by cremona

Replying to jdemeyer:

Replying to 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

comment:48 Changed 7 years ago by cremona

  • Status changed from needs_work to needs_review

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

comment:49 Changed 7 years ago by jdemeyer

  • Status changed from needs_review to positive_review

Further work: #18611

comment:50 follow-up: Changed 7 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: Changed 7 years ago by jdemeyer

Replying to 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?

comment:52 in reply to: ↑ 51 Changed 7 years ago by cremona

Replying to jdemeyer:

Replying to 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: Changed 7 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: Changed 7 years ago by cremona

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 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 7 years ago by jdemeyer

Replying to cremona:

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

comment:56 Changed 7 years ago by vbraun

  • Branch changed from u/cremona/18589 to 3d687e5225f67808eb6c5af5fbf4cb93f2000c62
  • Resolution set to fixed
  • Status changed from positive_review to closed
Note: See TracTickets for help on using tickets.