Faster Kohel isogenies without bivariate polynomials

This patch accelerates computation of Kohel formulas by replacing
internal bivariate polynomials k[x,y] by a tower of polynomial
rings k[x][y].

Because the y-coordinate of isogenies are always defined by
a polynomial of y-degree 1, this is equivalent to working
with a pair of univariate polynomials, which often have efficient
representations especially over finite fields.

The public API still exposes bivariate rational functions and is
not modified.

The resulting representation is several times faster.

Author: remyoudompheng

src/sage/schemes/elliptic_curves/ell_curve_isogeny.py

@@ -891,7 +891,7 @@ class EllipticCurveIsogeny(EllipticCurveHom):
     #
     __base_field = None
     __poly_ring = None # univariate in x over __base_field
-    __mpoly_ring = None # bivariate in x, y over __base_field
+    __mpoly_ring = None # __base_field[x][y], internal use only
 
     #
     # Rational Maps
@@ -931,7 +931,7 @@ class EllipticCurveIsogeny(EllipticCurveHom):
     #
     __psi = None # psi polynomial
     __phi = None # phi polynomial
-    __omega = None # omega polynomial
+    __omega = None # omega polynomial, an element of k[x][y]
 
 
     #
@@ -1466,7 +1466,7 @@ def __init_algebraic_structs(self, E):
             sage: phi._EllipticCurveIsogeny__poly_ring
             Univariate Polynomial Ring in x over Finite Field of size 17
             sage: phi._EllipticCurveIsogeny__mpoly_ring
-            Multivariate Polynomial Ring in x, y over Finite Field of size 17
+            Univariate Polynomial Ring in y over Univariate Polynomial Ring in x over Finite Field of size 17
 
         Now, calling the initialization function does nothing more::
 
@@ -1478,7 +1478,7 @@ def __init_algebraic_structs(self, E):
             sage: phi._EllipticCurveIsogeny__poly_ring
             Univariate Polynomial Ring in x over Finite Field of size 17
             sage: phi._EllipticCurveIsogeny__mpoly_ring
-            Multivariate Polynomial Ring in x, y over Finite Field of size 17
+            Univariate Polynomial Ring in y over Univariate Polynomial Ring in x over Finite Field of size 17
 
             sage: E = EllipticCurve(QQ, [0,0,0,1,0])
             sage: phi = EllipticCurveIsogeny(E, E((0,0)))
@@ -1490,7 +1490,7 @@ def __init_algebraic_structs(self, E):
             sage: phi._EllipticCurveIsogeny__poly_ring
             Univariate Polynomial Ring in x over Rational Field
             sage: phi._EllipticCurveIsogeny__mpoly_ring
-            Multivariate Polynomial Ring in x, y over Rational Field
+            Univariate Polynomial Ring in y over Univariate Polynomial Ring in x over Rational Field
 
             sage: F = GF(19); R.<x> = F[]
             sage: E = EllipticCurve(j=GF(19)(0))
@@ -1503,14 +1503,16 @@ def __init_algebraic_structs(self, E):
             sage: phi._EllipticCurveIsogeny__poly_ring
             Univariate Polynomial Ring in x over Finite Field of size 19
             sage: phi._EllipticCurveIsogeny__mpoly_ring
-            Multivariate Polynomial Ring in x, y over Finite Field of size 19
+            Univariate Polynomial Ring in y over Univariate Polynomial Ring in x over Finite Field of size 19
         """
         self._domain = E
         self.__base_field = E.base_ring()
         self.__poly_ring = PolynomialRing(self.__base_field, ['x'])
-        self.__mpoly_ring = PolynomialRing(self.__base_field, ['x','y'])
+        self.__mpoly_ring = PolynomialRing(self.__poly_ring, ['y'])
+        # The fraction fields are implicitly part of the public API, being the parents
+        # of the rational maps.
         self.__xfield = FractionField(self.__poly_ring)
-        self.__xyfield = FractionField(self.__mpoly_ring)
+        self.__xyfield = FractionField(PolynomialRing(self.__base_field, ['x', 'y']))
 
     def __compute_codomain(self):
         r"""
@@ -2046,7 +2048,7 @@ def __initialize_rational_maps_via_velu(self):
             ((x^4 + 5*x^3 + x^2 + 4*x)/(x^3 + 5*x^2 + 3*x + 5), (x^5*y - 2*x^3*y - x^2*y - 2*x*y + 2*y)/(x^5 + 3*x^3 + 3*x^2 + x - 1))
         """
         x = self.__poly_ring.gen()
-        y = self.__mpoly_ring.gen(1)
+        y = self.__xyfield.gen(1)
         return self.__compute_via_velu(x,y)
 
     def __init_kernel_polynomial_velu(self):
@@ -2197,7 +2199,7 @@ def __init_even_kernel_polynomial(self, E, psi_G):
             sage: from sage.schemes.elliptic_curves.ell_curve_isogeny import two_torsion_part
             sage: psig = two_torsion_part(E,x)
             sage: phi._EllipticCurveIsogeny__init_even_kernel_polynomial(E,psig)
-            (x^3 + 6*x, x^3*y + x*y, 6, 0, 1, 2)
+            (x^3 + 6*x, (x^3 + x)*y, 6, 0, 1, 2)
 
             sage: F = GF(2^4, 'alpha'); R.<x> = F[]
             sage: E = EllipticCurve(F, [1,1,0,1,0])
@@ -2206,7 +2208,7 @@ def __init_even_kernel_polynomial(self, E, psi_G):
 
             sage: psig = two_torsion_part(E,x)
             sage: phi._EllipticCurveIsogeny__init_even_kernel_polynomial(E,psig)
-            (x^3 + x, x^3*y + x^2 + x*y, 1, 0, 1, 2)
+            (x^3 + x, (x^3 + x)*y + x^2, 1, 0, 1, 2)
 
             sage: E = EllipticCurve(GF(7), [0,-1,0,0,1])
             sage: R.<x> = GF(7)[]
@@ -2216,7 +2218,7 @@ def __init_even_kernel_polynomial(self, E, psi_G):
             sage: psig = two_torsion_part(E,f)
             sage: psig = two_torsion_part(E,f)
             sage: phi._EllipticCurveIsogeny__init_even_kernel_polynomial(E,psig)
-            (x^7 + 5*x^6 + 2*x^5 + 6*x^4 + 3*x^3 + 5*x^2 + 6*x + 3, x^9*y - 3*x^8*y + 2*x^7*y - 3*x^3*y + 2*x^2*y + x*y - y, 1, 6, 3, 4)
+            (x^7 + 5*x^6 + 2*x^5 + 6*x^4 + 3*x^3 + 5*x^2 + 6*x + 3, (x^9 + 4*x^8 + 2*x^7 + 4*x^3 + 2*x^2 + x + 6)*y, 1, 6, 3, 4)
         """
         # check if the polynomial really divides the two_torsion_polynomial
         if self.__check and E.division_polynomial(2, x=self.__poly_ring.gen()) % psi_G  != 0 :
@@ -2228,7 +2230,7 @@ def __init_even_kernel_polynomial(self, E, psi_G):
         a1, a2, a3, a4, a6 = E.a_invariants()
         b2, b4, _, _ = E.b_invariants()
         x = self.__poly_ring.gen()
-        y = self.__mpoly_ring.gen(1)
+        y = self.__mpoly_ring.gen()
 
         if n == 1:
             x0 = -psi_G.constant_coefficient()
@@ -2307,7 +2309,7 @@ def __init_odd_kernel_polynomial(self, E, psi):
 
             sage: R.<x> = GF(7)[]
             sage: phi._EllipticCurveIsogeny__init_odd_kernel_polynomial(E, x+6)
-            (x^3 + 5*x^2 + 3*x + 2, x^3*y - 3*x^2*y + x*y, 2, 6, 1, 3)
+            (x^3 + 5*x^2 + 3*x + 2, (x^3 + 4*x^2 + x)*y, 2, 6, 1, 3)
 
             sage: F = GF(2^4, 'alpha'); R.<x> = F[]
             sage: alpha = F.gen()
@@ -2319,7 +2321,7 @@ def __init_odd_kernel_polynomial(self, E, psi):
             sage: R.<x> = F[]
             sage: f = x + alpha^2 + 1
             sage: phi._EllipticCurveIsogeny__init_odd_kernel_polynomial(E, f)
-            (x^3 + (alpha^2 + 1)*x + alpha^3 + alpha^2 + alpha, x^3*y + (alpha^2 + 1)*x^2*y + (alpha^2 + alpha + 1)*x^2 + (alpha^2 + 1)*x*y + (alpha^2 + alpha)*x + alpha*y + alpha, alpha^2 + alpha + 1, alpha^3 + alpha^2 + alpha, 1, 3)
+            (x^3 + (alpha^2 + 1)*x + alpha^3 + alpha^2 + alpha, (x^3 + (alpha^2 + 1)*x^2 + (alpha^2 + 1)*x + alpha)*y + (alpha^2 + alpha + 1)*x^2 + (alpha^2 + alpha)*x + alpha, alpha^2 + alpha + 1, alpha^3 + alpha^2 + alpha, 1, 3)
 
             sage: E = EllipticCurve(j=-262537412640768000)
             sage: f = E.isogenies_prime_degree()[0].kernel_polynomial()
@@ -2407,11 +2409,12 @@ def __compute_omega_fast(self, E, psi, psi_pr, phi, phi_pr):
             sage: fi = phi._EllipticCurveIsogeny__phi
             sage: fi_pr = fi.derivative()
             sage: phi._EllipticCurveIsogeny__compute_omega_fast(E, psi, psi_pr, fi, fi_pr)
-            x^3*y - 3*x^2*y + x*y
+            (x^3 + 4*x^2 + x)*y
         """
         a1 = E.a1()
         a3 = E.a3()
-        x, y = self.__mpoly_ring.gens()
+        x = self.__poly_ring.gen()
+        y = self.__mpoly_ring.gen()
 
         psi_2 = 2*y + a1*x + a3
 
@@ -2456,7 +2459,7 @@ def __compute_omega_general(self, E, psi, psi_pr, phi, phi_pr):
             sage: fi = phi._EllipticCurveIsogeny__phi
             sage: fi_pr = fi.derivative()
             sage: phi._EllipticCurveIsogeny__compute_omega_general(E, psi, psi_pr, fi, fi_pr)
-            x^3*y + (alpha^2 + 1)*x^2*y + (alpha^2 + alpha + 1)*x^2 + (alpha^2 + 1)*x*y + (alpha^2 + alpha)*x + alpha*y + alpha
+            (x^3 + (alpha^2 + 1)*x^2 + (alpha^2 + 1)*x + alpha)*y + (alpha^2 + alpha + 1)*x^2 + (alpha^2 + alpha)*x + alpha
 
         A bug fixed in :trac:`7907`::
 
@@ -2470,7 +2473,8 @@ def __compute_omega_general(self, E, psi, psi_pr, phi, phi_pr):
         """
         a1, a2, a3, a4, a6 = E.a_invariants()
         b2, b4, _, _ = E.b_invariants()
-        x, y = self.__mpoly_ring.gens()
+        x = self.__poly_ring.gen()
+        y = self.__mpoly_ring.gen()
 
         n = psi.degree()
         d = 2 * n + 1
@@ -2559,7 +2563,9 @@ def __compute_via_kohel(self, xP, yP):
             ((x^3 - 2*x^2 + 3*x + 2)/(x^2 - 2*x + 1), (x^3*y - 3*x^2*y + x*y)/(x^3 - 3*x^2 + 3*x - 1))
         """
         a = self.__phi(xP)
-        b = self.__omega(xP, yP)
+        omega0 = self.__omega[0]
+        omega1 = self.__omega[1]
+        b = omega0(xP) + omega1(xP)*yP
         c = self.__psi(xP)
         return a/c**2, b/c**3