gh-36637: sums of elliptic-curve morphisms
    
Here we add a new `EllipticCurveHom` child class representing formal
sums of elliptic-curve morphisms with the same domain and codomain. One
of the main features is a conversion function from formal sums to a more
explicit representation of the same morphism as an isogeny chain, which
will come in handy for endomorphism-ring-related computations.

In the process, we also generalize `point_of_order()` from primes to
prime powers.

There is certainly a lot of potential for optimizations; some ideas are
mentioned in the code. Since I think it's fair to say that this version
is already much better than what we currently have for sums of isogenies
(i.e., nothing), I'd suggest merging this as is and dealing with
bottlenecks later as they arise.

Cc: @JohnCremona @defeo @GiacomoPope @remyoudompheng
    
URL: #36637
Reported by: Lorenz Panny
Reviewer(s): John Cremona, Lorenz Panny

Author: Release Manager

src/doc/en/reference/arithmetic_curves/index.rst

@@ -18,10 +18,11 @@ Maps between them
    :maxdepth: 1
 
    sage/schemes/elliptic_curves/hom
+   sage/schemes/elliptic_curves/hom_composite
+   sage/schemes/elliptic_curves/hom_sum
    sage/schemes/elliptic_curves/weierstrass_morphism
    sage/schemes/elliptic_curves/ell_curve_isogeny
    sage/schemes/elliptic_curves/hom_velusqrt
-   sage/schemes/elliptic_curves/hom_composite
    sage/schemes/elliptic_curves/hom_scalar
    sage/schemes/elliptic_curves/hom_frobenius
    sage/schemes/elliptic_curves/isogeny_small_degree

src/sage/schemes/elliptic_curves/ell_curve_isogeny.py

@@ -2879,28 +2879,19 @@ def kernel_polynomial(self):
             self.__init_kernel_polynomial()
         return self.__kernel_polynomial
 
-    def is_separable(self):
+    def inseparable_degree(self):
         r"""
-        Determine whether or not this isogeny is separable.
+        Return the inseparable degree of this isogeny.
 
-        Since :class:`EllipticCurveIsogeny` only implements
-        separable isogenies, this method always returns ``True``.
+        Since this class only implements separable isogenies,
+        this method always returns one.
 
-        EXAMPLES::
-
-            sage: E = EllipticCurve(GF(17), [0,0,0,3,0])
-            sage: phi = EllipticCurveIsogeny(E,  E((0,0)))
-            sage: phi.is_separable()
-            True
-
-        ::
+        TESTS::
 
-            sage: E = EllipticCurve('11a1')
-            sage: phi = EllipticCurveIsogeny(E, E.torsion_points())
-            sage: phi.is_separable()
-            True
+            sage: EllipticCurveIsogeny.inseparable_degree(None)
+            1
         """
-        return True
+        return Integer(1)
 
     def _set_pre_isomorphism(self, preWI):
         """

src/sage/schemes/elliptic_curves/ell_field.py

@@ -2082,13 +2082,13 @@ def compute_model(E, name):
 
     raise NotImplementedError(f'cannot compute {name} model')
 
-def point_of_order(E, l):
+def point_of_order(E, n):
     r"""
     Given an elliptic curve `E` over a finite field or a number field
-    and an integer `\ell \geq 1`, construct a point of order `\ell` on `E`,
+    and an integer `n \geq 1`, construct a point of order `n` on `E`,
     possibly defined over an extension of the base field of `E`.
 
-    Currently only prime values of `\ell` are supported.
+    Currently only prime powers `n` are supported.
 
     EXAMPLES::
 
@@ -2103,6 +2103,13 @@ def point_of_order(E, l):
         sage: P.curve().a_invariants()
         (1, 2, 3, 4, 5)
 
+    ::
+
+        sage: Q = point_of_order(E, 8); Q
+        (69*x^5 + 24*x^4 + 100*x^3 + 65*x^2 + 88*x + 97 : 65*x^5 + 28*x^4 + 5*x^3 + 45*x^2 + 42*x + 18 : 1)
+        sage: 8*Q == 0 and 4*Q != 0
+        True
+
     ::
 
         sage: from sage.schemes.elliptic_curves.ell_field import point_of_order
@@ -2111,10 +2118,23 @@ def point_of_order(E, l):
         (x : -Y : 1)
         sage: P.base_ring()
         Number Field in Y with defining polynomial Y^2 - x^3 - 7*x - 7 over its base field
+        sage: P.base_ring().base_field()
+        Number Field in x with defining polynomial x^4 + 14*x^2 + 28*x - 49/3
         sage: P.order()
         3
         sage: P.curve().a_invariants()
         (0, 0, 0, 7, 7)
+
+    ::
+
+        sage: Q = point_of_order(E, 4); Q  # random
+        (x : Y : 1)
+        sage: Q.base_ring()
+        Number Field in Y with defining polynomial Y^2 - x^3 - 7*x - 7 over its base field
+        sage: Q.base_ring().base_field()
+        Number Field in x with defining polynomial x^6 + 35*x^4 + 140*x^3 - 245*x^2 - 196*x - 735
+        sage: Q.order()
+        4
     """
     # Construct the field extension defined by the given polynomial,
     # in such a way that the result is recognized by Sage as a field.
@@ -2127,14 +2147,16 @@ def ffext(poly):
             return poly.splitting_field(rng.variable_name())
         return fld.extension(poly, rng.variable_name())
 
-    l = ZZ(l)
-    if l == 1:
+    n = ZZ(n)
+    if n == 1:
         return E(0)
 
-    if not l.is_prime():
-        raise NotImplementedError('composite orders are currently unsupported')
+    l,m = n.is_prime_power(get_data=True)
+    if not m:
+        raise NotImplementedError('only prime-power orders are currently supported')
 
-    xpoly = E.division_polynomial(l)
+    xpoly = E.division_polynomial(n).radical()
+    xpoly //= E.division_polynomial(n//l).radical()
     if xpoly.degree() < 1:  # supersingular and l == p
         raise ValueError('curve does not have any points of the specified order')
 
@@ -2149,4 +2171,6 @@ def ffext(poly):
         xx = FF(xx)
 
     EE = E.change_ring(FF)
-    return EE.lift_x(xx)
+    pt = EE.lift_x(xx)
+    pt.set_order(n, check=False)
+    return pt

src/sage/schemes/elliptic_curves/ell_point.py

@@ -1359,9 +1359,10 @@ def set_order(self, value=None, *, multiple=None, check=True):
                 raise ValueError('Value %s illegal for point order' % value)
             E = self.curve()
             q = E.base_ring().cardinality()
-            low, hi = Hasse_bounds(q)
-            if value > hi:
-                raise ValueError('Value %s illegal: outside max Hasse bound' % value)
+            if q < oo:
+                _, hi = Hasse_bounds(q)
+                if value > hi:
+                    raise ValueError('Value %s illegal: outside max Hasse bound' % value)
             if value * self != E(0):
                 raise ValueError('Value %s illegal: %s * %s is not the identity' % (value, value, self))
             if hasattr(self, '_order') and self._order != value:  # already known