implement sparse strategy for prime-power isogenies

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

@@ -1980,6 +1980,11 @@ REFERENCES:
 .. [DeVi1984] \M.-P. Delest, and G. Viennot, *Algebraic Languages and
               Polyominoes Enumeration.* Theoret. Comput. Sci. 34, 169-206, 1984.
 
+.. [DJP2014] Luca De Feo, David Jao and Jérôme Plût: Towards quantum-resistant
+             cryptosystems from supersingular elliptic curve isogenies. Journal
+             of Mathematical Cryptology, vol. 8, no. 3, 2014, pp. 209-247.
+             https://eprint.iacr.org/2011/506.pdf
+
 .. [DFMS1996] Philipppe Di Francesco, Pierre Mathieu, and David Sénéchal.
               *Conformal Field Theory*. Graduate Texts in Contemporary
               Physics, Springer, 1996.

src/sage/schemes/elliptic_curves/hom_composite.py

@@ -89,8 +89,6 @@
 from sage.schemes.elliptic_curves.ell_curve_isogeny import EllipticCurveIsogeny
 from sage.schemes.elliptic_curves.weierstrass_morphism import WeierstrassIsomorphism, identity_morphism
 
-# TODO: Implement sparse strategies? (cf. the SIKE cryptosystem)
-
 
 def _eval_factored_isogeny(phis, P):
     """
@@ -117,33 +115,95 @@ def _eval_factored_isogeny(phis, P):
     return P
 
 
-def _compute_factored_isogeny_prime_power(P, l, e):
-    """
-    This method takes a point `P` of order `l^e` and returns
-    a sequence of degree-`l` isogenies whose composition has
+def _compute_factored_isogeny_prime_power(P, l, n, split=.8):
+    r"""
+    This method takes a point `P` of order `\ell^n` and returns
+    a sequence of degree-`\ell` isogenies whose composition has
     the subgroup generated by `P` as its kernel.
 
+    The optional argument ``split``, a real number between
+    `0` and `1`, controls the *strategy* used to compute the
+    isogeny: In general, the algorithm performs a number of
+    scalar multiplications `[\ell]` and a number of
+    `\ell`-isogeny evaluations, and there exist tradeoffs
+    between them.
+
+    - Setting ``split`` to `0` skews the algorithm towards
+      isogenies, minimizing multiplications.
+      The asymptotic complexity is `O(n \log(\ell) + n^2 \ell)`.
+
+    - Setting ``split`` to `1` skews the algorithm towards
+      multiplications, minimizing isogenies.
+      The asymptotic complexity is `O(n^2 \log(\ell) + n \ell)`.
+
+    - Values strictly between `0` and `1` define *sparse*
+      strategies, which balance the number of isogenies and
+      multiplications according to the parameter.
+      The asymptotic complexity is `O(n \log(n) \ell)`.
+
+    .. NOTE::
+
+        As of July 2022, good values for ``split`` range somewhere
+        between roughly `0.6` and `0.9`, depending on the size of
+        `\ell` and the cost of base-field arithmetic.
+
+    REFERENCES:
+
+    Sparse strategies were introduced in [DJP2014]_, §4.2.2.
+
+    ALGORITHM:
+
+    The splitting rule implemented here is a simple heuristic which
+    is usually not optimal. The advantage is that it does not rely
+    on prior knowledge of degrees and exact costs of elliptic-curve
+    arithmetic, while still working reasonably well for a fairly
+    wide range of situations.
+
     EXAMPLES::
 
         sage: from sage.schemes.elliptic_curves import hom_composite
         sage: E = EllipticCurve(GF(8191), [1,0])
         sage: P = E.random_point()
-        sage: (l,e), = P.order().factor()
-        sage: phis = hom_composite._compute_factored_isogeny_prime_power(P,l,e)
+        sage: (l,n), = P.order().factor()
+        sage: phis = hom_composite._compute_factored_isogeny_prime_power(P,l,n)
         sage: hom_composite._eval_factored_isogeny(phis, P)
         (0 : 1 : 0)
-        sage: [phi.degree() for phi in phis] == [l]*e
+        sage: [phi.degree() for phi in phis] == [l]*n
+        True
+
+    All choices of ``split`` produce the same result, albeit
+    not equally fast::
+
+        sage: E = EllipticCurve(GF(2^127 - 1), [1,0])
+        sage: P, = E.gens()
+        sage: (l,n), = P.order().factor()
+        sage: phis = hom_composite._compute_factored_isogeny_prime_power(P,l,n)
+        sage: for s in (0, 0.1, 0.5, 0.9, 1):
+        ....:     phis == hom_composite._compute_factored_isogeny_prime_power(P,l,n, split=s)
+        True
+        True
+        True
+        True
         True
     """
-    E = P.curve()
-    phis = []
-    for i in range(e):
-        K = l**(e-1-i) * P
-        phi = EllipticCurveIsogeny(E, K)
-        E = phi.codomain()
-        P = phi(P)
-        phis.append(phi)
-    return phis
+    def rec(Q, k):
+        assert k
+        if k == 1:  # base case
+#            assert Q and not l*Q
+            return [EllipticCurveIsogeny(Q.curve(), Q)]
+
+        k1 = int(k * split + .5)
+        k1 = max(1, min(k-1, k1))  # clamp to [1, k-1]
+
+        Q1 = l**k1 * Q
+        L = rec(Q1, k-k1)
+
+        Q2 = _eval_factored_isogeny(L, Q)
+        R = rec(Q2, k1)
+
+        return L + R
+
+    return rec(P, n)
 
 
 def _compute_factored_isogeny_single_generator(P):
@@ -193,7 +253,7 @@ def _compute_factored_isogeny(kernel):
     phis = []
     ker = list(kernel)
     while ker:
-        K, ker = ker[0], ker[1:]
+        K = ker.pop(0)
         psis = _compute_factored_isogeny_single_generator(K)
         ker = [_eval_factored_isogeny(psis, P) for P in ker]
         phis += psis