This worksheet is by Luca De Feo.

The implementation of isogenies was started in Sage 4.1 by Dan Shumow, then taken up by many authors. This worksheet tries to review the known limits of the current implementation (linking to relevant trac tickets), point out its weaknesses and suggest improvements.

The idea is to have a work plan for the upcoming Sage Days that will take place during the ECC summer school in Nancy on Sept. 12-16.

Current status

So, here's the status of morphisms of elliptic curves. Sage implements isomorphisms of Weierstrass models:

Sage also implements isogenies. For example, if  a generator of the kernel is given, a canonical isogeny is computed via Vélu formulae:

Alternatively we can pass the polynomial vanishing on the abscissae of the kernel, Vélu's formulae are applied in this case too. And here's the first thing I hate: notice the difference between the timings. This is because two completely separate implementations of Vélu's formulae are used in the two cases (see Dan Shumow's thesis).

Unfortunately, the implementation of Vélu's formulae from the kernel polynomial is incomplete and does not handle even-degree isogenies other than degree 2 ones and the multiplication-by-2 map.

We can test isogenies for equality. Isogenies are not unique objects, though; but elliptic curves aren't either, so I guess that this problem can wait.

however, equality test is buggy: only the __hash__ method takes consideration of the domain and codomain (thus == will work as expected), but the __cmp__ function only compares kernel polynomials, really! I should make a ticket for this.

Because of the poor implementation of the __cmp__ method, I suspect there are other ways of making equality tests fail.

Isogenies internally use WeierstrassIsomorphism to reduce to some canonical forms, but here's one of the things I hate the most: isogenies cannot be composed with isomorphisms.

An important subclass of isogenies are the endomorphisms. Unfortunately what Sage really implements are Vélu isogenies (also called normalized isogenies), and endomorphisms are in general not normalized. I am currently working on this and will soon post a patch that allows to construct non-normalized isogenies.

Partly because of this, there is no endomorphism object for elliptic curves and there are two different functions to compute the multiplication-by-m map. Yes, contraddicting what I said before, the second returns a non-normalized isogeny; the point is that this and EllipticCurveIsogeny.dual() are the only two ways of constructing non-normalized isogenies.

Also observe the huge difference between the timings (don't even consider computing the multiplication-by-13!). I don't even understand what exactly takes so long, I believe there might be a serious bug in the way Vélu's formulae are implemented. This issue is already pointed out in ticket #8014 and is a serious obstacle to seriously implementing an endomorphism object. The solution would be to store only the information relative to the kernel of the isogeny, and to lazily compute all the other quantities upon request.

It is also worth mentioning that the kernel of a normalized isogeny can be computed from the knowledge of its domain, codomain and degree. Currently only Stark's algorithm is implemented; ticket 11095 adds the BMSS algorithm, which is faster and slightly more general, and does some cleanup.

Unfortunately, the situation where the isogeny is not normalized is more complex. As said before, I am currently working on this and shall soon post a patch.

Another major feature missing in sage is unseparable isogenies: there's simply no way of constructing them. Besides that, the multiplication_by_m method also fails to compute unseparable endomorphisms. See also ticket 6413.

Finally, it should also be mentioned that Sage does not implement any algorithm to compute the kernel polynomial of an isogeny from its domain and codomain (even when the isogeny is normalized) in case the characteristic of the base field is too small. Ticket 11094 aims at solving this problem at least for finite fields, but these algorithms are quite complex and have some big prerequisites; in particular:

• Lercier-Sirvent requires arithmetic in unramified extensions of \mathbb{Q}_p, in particular Hensel lifting;
• Couveignes '96 requires (at least) finite field lattices (see ticket 8335);
• Couveignes '94 should be feasible (although slow) in the current version.

Implementing these algorithms is a long-term goal and shouldn't be considered for the next Sage Days.

Now I will go through some proposed improvements to the current situation. Each of these should translate to one or more tickets to be implemented at the Sage Days.

The computational part: efficient unique representation of isogenies via Vélu's formulae

There are many ways of representing (aka computing) an isogeny, but here we are mainly interested in what is usually called explicit isogenies. Given an isogeny between two curves in (affine) Weierstrass form, it is always possible to write it as a pair of rational fractions

If the isogeny is separable, Vélu's formulae show that a monic polynomial h uniquely determines g, p_1, p_2 and q up to isomorphism of the codomain curve. Vélu's formulae are not the only possible canonical choice to represent separable isogenies, but they are historically relevant and have some advantages over other choices. Their main disadvantage is that the Vélu isogeny associated to the m-torsion subgroup is not the multiplication-by-m endomorphism (in general, not even after rescaling).

If the isogeny is unseparable, then it is the composition of a Frobenius map followed by a separable isogeny (Silverman AEC, Corollary II.2.12, I believe this is also true for non-perfect fields, right?). Hence any isogeny can be uniquely decomposed as

So, I propose to implement isogenies with a class looking like this:

Notice that this representation encompasses well all sort of isogenies:

• non-normalized isogenies (in this case u\ne1),
• endomorphisms (for example, for the multiplication-by-m map we have h vanishing on the m-torsion, u=m plus some other conditions on r,s,t),
• isomorphisms (in this case h=1),
• purely inseparable isogenies (in this case h=1 and u_degree>0).

However, the constant zero isogeny requires a special treatment: we can use h=0, but this does not uniquely define the codomain, which must be given explicitely.

For efficiency, no big computation should be performed in the __init__ method. The rational maps can be computed as lazy attributes using Vélu's formulae; if implemented correctly, this can be almost as fast as computing the multiplication-by-m maps starting from the m-division polynomial.

For a good implementation of Vélu's formulae, I suggest following these guidelines:

• Use as few case splittings as possible. This is good for maintainability, readability and ease of debugging.
• Frown upon polynomial GCDs (and thus taking square roots of polynomials), because not all fields implement them (or not efficiently, at least). [Bostan, Morain, Salvy, Schost, Fast algorithms for computing isogenies between elliptic curves, Prop 4.1] gives formulas that need no GCDs (they only work in characteristic not 2 or 3, but they are not too hard to generalise).
  
• Avoid field inversions as much as possible: we want to keep the ability to compute multiplication-by-m maps over the integers.

On the inputs to the constructor

Some people may not be happy with my choice of inputs to the __init__ method: while u_degree and urst will probably please anyone, h may not be consensual. The polynomial h, in fact, factors as a product ab^2, where a|(x^3+b_2x^2+b_4x+b_6) vanishes on the abscissae of the 2-torsion part of the kernel, and b vanishes on all the other abscissae with multiplicity exactly one. Following the convention for division polynomials, it is traditional to input the polynomial ab, instead of h to Vélu formulae (this is the input currently expected in Sage, for example). Unfortunately, I don't think there's any formula using this input that can avoid a GCD between ab and x^3+b_2x^2+b_4x+b_6. For maximum compatibility, and to make everyone happy, I suggest the following convention:

• E.isogeny(kernel_polynomial=h) assumes h is given as above
• E.isogeny(kernel_polynomial=(a,b)) assumes h=ab^2 (this could be possibly used for efficiency)
• E.isogeny(kernel_polynomial=h, degree=n) compares the degree of h with n-1, if they are different it uses GCDs to see if h is of the form ab with a|(x^3+b_2x_2+b_4x+b_6).

We may also allow to input an inseparable polynomial h. For example, if h=H(x^q) for some separable polynomial H, we can consider that the user is asking to compose the q-th power Frobenius with the separable isogeny of E^{(q)} having kernel polynomial H.

For compatibility with division polynomials, we may also accept a non-monic h. If c is the leading coefficient of h, this would be equivalent to create an isogeny with parameters h'=h/c and (u,r,s,t)=(\sqrt{c},0,0,0). If c has no square root, we could extend the base field or move to a quadratic twist, although it is probably better to raise an exception.

Constructing isogenies from kernel points

The current implementation allows to input a list of kernel points instead of a kernel polynomial to the constructor. Although seldom used in real-world computations, and limited to separable isogenies, this feature is often handy and I suggest keeping it:

• The current implementation uses different formulae, depending whether the input is a kernel polynomial or kernel points. I think it would be better to always reduce to the implementation based on the kernel polynomial, as it is more general and more efficient. Supposing only one point is given, verify that it has finite order, enumerate its multiples, and construct the polynomial h vanishing on their abscissae.
• If a list of points is given, the current implementation constructs the subgroup generated by all the points, then applies Vélu's formulas; I think it would be more efficient to construct a chain of isogenies: construct a Vélu isogeny for the first point, push the other points through the isogeny, discard points sent to zero, then start again until no point is left. The result would be a composition of all the isogenies, or maybe a Factorization object (which would use much less resources).
• Maybe, we could also add a "factor" option, so that the order of each single point is factored and a complete factorization into isogenies of prime degree is returned. (Nota: this could be done for kernel polynomials too, but it would be much slower and the factorization wouldn't necessary be into prime factors)

The engineering part: class hierarchy

The current implementation of EllipticCurveIsogeny inherits from sage.categories.morphism.Morphism. This is not very tight and looses most of the benefits of the coercion model. I'd like some advice from coercion model experts, here. I think that isogenies should inherit at least from

• sage.structure.element.AdditiveGroupElement (this, of course, means that we have to implement addition)
• sage.categories.morphism.Morphism

but maybe sage.schemes contains some possibly better superclasses:

• sage.schemes.generic.morphism.SchemeMorphism,
• sage.schemes.generic.morphism.SchemeMorphism_abelian_variety_coordinates_field

By looking at sage.schemes.elliptic_curves.ell_point, it seems that at some point it was considered making EllipticCurvePoint_field inherit from it, but then this has been abandoned in favour AdditiveGroupElement.

Isogenies with given domain and codomain should have a parent. Currently this is

which seems a reasonable choice, but this should be extended to inherit from AbelianGroup at least. Properly speaking, this parent is \mathrm{Hom}_K(E_1,E_2), where K is the base field of the elliptic curves.

I think it would be useful to implement some subclasses of isogenies and their respective parents:

• Endomorphisms
• Isomorphisms
• Automorphisms

with methods and optimizations specific to each class. For example, the endomorphisms should form a ring, containing the automorphisms as its unit group.

Some more interesting methods could be added to these classes. For example, currently the Frobenius endomorphism is given as the generator of an order in an imaginary quadratic field

in this implementation, the Frobenius could be given as an element of the endomorphism ring, instead, and its parent could implement a method 'minimal_order' (I cannot think of a better name right now) returning the order above, and registering a coercion from this order to the endomorphism ring (and maybe a conversion in the other direction, but this is computationally harder). In the future, a method "order" could also be implemented, computing the order isomorphic to End(E).

Concluding remarks

I hope to have given a sound and practical specification for isogenies. This document probably contains errors and inconstincencies. I'd appreciate feedback, especially on the class hierarchy.

The degree to which this project will be completed during the summer school will largely depend on the programming skills of the students, but I'd be surprised if at the end of the week we had a patch mature enough to be merged. If this project succeeds, I'll take care of coordinating the work after the summer school and uploading the relevant patches.