Initial Commit

Co-Authored-By: Damien Robert <Damien.Olivier.Robert+git@gmail.com>
Co-Authored-By: Frederik Vercauteren <frederik.vercauteren@gmail.com>
Co-Authored-By: Jonathan Komada Eriksen <jonathan.eriksen97@gmail.com>
Co-Authored-By: Pierrick Dartois <pierrickdartois@icloud.com>
Co-Authored-By: Riccardo Invernizzi <nidadoni@gmail.com>
Co-Authored-By: Ryan Rueger <git@rueg.re>                                                                                                                                           [0.01s]
Co-Authored-By: Benjamin Wesolowski <benjamin@pasch.umpa.ens-lyon.fr>
Co-Authored-By: Arthur Herlédan Le Merdy <ahlm@riseup.net>
Co-Authored-By: Boris Fouotsa <tako.fouotsa@epfl.ch>

1 files changed, 440 insertions, 0 deletions

diff --git a/theta_lib/theta_structures/Theta_dim2.py b/theta_lib/theta_structures/Theta_dim2.py
new file mode 100644
index 0000000..0d9955a
--- /dev/null
+++ b/theta_lib/theta_structures/Theta_dim2.py
@@ -0,0 +1,440 @@
+# Sage Imports
+from sage.all import Integer
+from sage.structure.element import get_coercion_model, RingElement
+
+cm = get_coercion_model()
+
+from .theta_helpers_dim2 import batch_inversion, product_theta_point
+from .Theta_dim1 import ThetaStructureDim1
+from .Tuple_point import TuplePoint
+from ..basis_change.base_change_dim2 import apply_base_change_theta_dim2
+
+# =========================================== #
+#     Class for Theta Structure (level-2)     #
+# =========================================== #
+
+
+class ThetaStructureDim2:
+    """
+    Class for the Theta Structure in dimension 2, defined by its theta null point. This type
+    represents the generic domain/codomain of the (2,2)-isogeny in the theta model.
+    """
+
+    def __init__(self, null_point, null_point_dual=None):
+        if not len(null_point) == 4:
+            raise ValueError
+
+        self._base_ring = cm.common_parent(*(c.parent() for c in null_point))
+        self._point = ThetaPointDim2
+        self._precomputation = None
+
+        self._null_point = self._point(self, null_point)
+        self._null_point_dual = null_point_dual
+
+    def null_point(self):
+        """
+        Return the null point of the given theta structure
+        """
+        return self._null_point
+
+    def null_point_dual(self):
+        if self._null_point_dual==None:
+            self._null_point_dual = self._point.to_hadamard(*self.coords())
+        return self._null_point_dual
+
+    def base_ring(self):
+        """
+        Return the base ring of the common parent of the coordinates of the null point
+        """
+        return self._base_ring
+
+    def zero(self):
+        """
+        The additive identity is the theta null point
+        """
+        return self.null_point()
+
+    def zero_dual(self):
+        return self.null_point_dual()
+
+    def __repr__(self):
+        return f"Theta structure over {self.base_ring()} with null point: {self.null_point()}"
+
+    def coords(self):
+        """
+        Return the coordinates of the theta null point of the theta structure
+        """
+        return self.null_point().coords()
+
+    def hadamard(self):
+        """
+        Compute the Hadamard transformation of the theta structure
+        """
+        return ThetaStructureDim2(self.null_point_dual(),null_point_dual=self.coords())
+
+    def squared_theta(self):
+        """
+        Square the coefficients and then compute the Hadamard transformation of
+        the theta null point of the theta structure
+        """
+        return self.null_point().squared_theta()
+
+    def _arithmetic_precomputation(self):
+        """
+        Precompute 6 field elements used in arithmetic and isogeny computations
+        """
+        if self._precomputation is None:
+            a, b, c, d = self.null_point().coords()
+
+            # Technically this computes 4A^2, 4B^2, ...
+            # but as we take quotients this doesnt matter
+            # Cost: 4S
+            AA, BB, CC, DD = self.squared_theta()
+
+            # Precomputed constants for addition and doubling
+            b_inv, c_inv, d_inv, BB_inv, CC_inv, DD_inv = batch_inversion([
+                b, c, d, BB, CC, DD]
+            )
+
+            y0 = a * b_inv
+            z0 = a * c_inv
+            t0 = a * d_inv
+
+            Y0 = AA * BB_inv
+            Z0 = AA * CC_inv
+            T0 = AA * DD_inv
+
+            self._precomputation = (y0, z0, t0, Y0, Z0, T0)
+        return self._precomputation
+
+    def __call__(self, coords):
+        return self._point(self, coords)
+
+    def base_change_struct(self,N):
+        null_coords=self.null_point().coords()
+        new_null_coords=apply_base_change_theta_dim2(N,null_coords)
+        return ThetaStructure(new_null_coords)
+
+    def base_change_coords(self,N,P):
+        coords=P.coords()
+        new_coords=apply_base_change_theta_dim2(N,coords)
+        return self.__call__(new_coords)
+
+
+# =================================================== #
+#     Class for Product Theta Structure (level-2)     #
+# =================================================== #
+
+
+class ProductThetaStructureDim2(ThetaStructureDim2):
+    def __init__(self,*args):
+        r"""Defines the product theta structure at level 2 of 2 elliptic curves.
+
+        Input: Either
+        - 2 theta structures of dimension 1: T0, T1;
+        - 2 elliptic curves: E0, E1.
+        - 2 elliptic curves E0, E1 and their respective canonical 4-torsion basis B0, B1.
+        """
+        if len(args)==2:
+            theta_structures=list(args)
+            for k in range(2):
+                if not isinstance(theta_structures[k],ThetaStructureDim1):
+                    theta_structures[k]=ThetaStructureDim1(theta_structures[k])
+        elif len(args)==4:
+            theta_structures=[ThetaStructureDim1(args[k],args[2+k][0],args[2+k][1]) for k in range(2)]
+        else:
+            raise ValueError("2 or 4 arguments expected but {} were given.\nYou should enter a list of 2 elliptic curves or ThetaStructureDim1\nor a list of 2 elliptic curves with a 4-torsion basis for each of them.".format(len(args)))
+
+        self._theta_structures=theta_structures
+
+        null_point=product_theta_point(theta_structures[0].zero().coords(),theta_structures[1].zero().coords())
+
+        ThetaStructureDim2.__init__(self,null_point)
+
+    def product_theta_point(self,theta_points):
+        t0,t1=theta_points[0].coords()
+        u0,u1=theta_points[1].coords()
+        return self._point(self,[t0*u0,t1*u0,t0*u1,t1*u1])
+
+    def __call__(self,point):
+        if isinstance(point,TuplePoint):
+            theta_points=[]
+            theta_structures=self._theta_structures
+            for i in range(2):
+                theta_points.append(theta_structures[i](point[i]))
+            return self.product_theta_point(theta_points)
+        else:
+            return self._point(self,point)
+
+    def to_theta_points(self,P):
+        coords=P.coords()
+        theta_coords=[(coords[0],coords[1]),(coords[1],coords[3])]
+        theta_points=[self._theta_structures[i](theta_coords[i]) for i in range(2)]
+        return theta_points
+
+    def to_tuple_point(self,P):
+        theta_points=self.to_theta_points(P)
+        montgomery_points=[self._theta_structures[i].to_montgomery_point(theta_points[i]) for i in range(2)]
+        return TuplePoint(montgomery_points)
+
+
+# ======================================= #
+#     Class for Theta Point (level-2)     #
+# ======================================= #
+
+
+class ThetaPointDim2:
+    """
+    A Theta Point in the level-2 Theta Structure is defined with four projective
+    coordinates
+
+    We cannot perform arbitrary arithmetic, but we can compute doubles and
+    differential addition, which like x-only points on the Kummer line, allows
+    for scalar multiplication
+    """
+
+    def __init__(self, parent, coords):
+        if not isinstance(parent, ThetaStructureDim2) and not isinstance(parent, ProductThetaStructureDim2):
+            raise ValueError
+
+        self._parent = parent
+        self._coords = tuple(coords)
+
+        self._hadamard = None
+        self._squared_theta = None
+
+    def parent(self):
+        """
+        Return the parent of the element, of type ThetaStructureDim2
+        """
+        return self._parent
+
+    def theta(self):
+        """
+        Return the parent theta structure of this ThetaPointDim2"""
+        return self.parent()
+
+    def coords(self):
+        """
+        Return the projective coordinates of the ThetaPointDim2
+        """
+        return self._coords
+
+    def is_zero(self):
+        """
+        An element is zero if it is equivalent to the null point of the parent
+        ThetaStrcuture
+        """
+        return self == self.parent().zero()
+
+    @staticmethod
+    def to_hadamard(x_00, x_10, x_01, x_11):
+        """
+        Compute the Hadamard transformation of four coordinates, using recursive
+        formula.
+        """
+        x_00, x_10 = (x_00 + x_10, x_00 - x_10)
+        x_01, x_11 = (x_01 + x_11, x_01 - x_11)
+        return x_00 + x_01, x_10 + x_11, x_00 - x_01, x_10 - x_11
+
+    def hadamard(self):
+        """
+        Compute the Hadamard transformation of this element
+        """
+        if self._hadamard is None:
+            self._hadamard = self.to_hadamard(*self.coords())
+        return self._hadamard
+
+    @staticmethod
+    def to_squared_theta(x, y, z, t):
+        """
+        Square the coordinates and then compute the Hadamard transform of the
+        input
+        """
+        return ThetaPointDim2.to_hadamard(x * x, y * y, z * z, t * t)
+
+    def squared_theta(self):
+        """
+        Compute the Squared Theta transformation of this element
+        which is the square operator followed by Hadamard.
+        """
+        if self._squared_theta is None:
+            self._squared_theta = self.to_squared_theta(*self.coords())
+        return self._squared_theta
+
+    def double(self):
+        """
+        Computes [2]*self
+
+        NOTE: Assumes that no coordinate is zero
+
+        Cost: 8S 6M
+        """
+        # If a,b,c,d = 0, then the codomain of A->A/K_2 is a product of
+        # elliptic curves with a non product theta structure.
+        # Unless we are very unlucky, A/K_1 will not be in this case, so we
+        # just need to Hadamard, double, and Hadamard inverse
+        # If A,B,C,D=0 then the domain itself is a product of elliptic
+        # curves with a non product theta structure. The Hadamard transform
+        # will not change this, we need a symplectic change of variable
+        # that puts us back in a product theta structure
+        y0, z0, t0, Y0, Z0, T0 = self.parent()._arithmetic_precomputation()
+
+        # Temp coordinates
+        # Cost 8S 3M
+        xp, yp, zp, tp = self.squared_theta()
+        xp = xp**2
+        yp = Y0 * yp**2
+        zp = Z0 * zp**2
+        tp = T0 * tp**2
+
+        # Final coordinates
+        # Cost 3M
+        X, Y, Z, T = self.to_hadamard(xp, yp, zp, tp)
+        X = X
+        Y = y0 * Y
+        Z = z0 * Z
+        T = t0 * T
+
+        coords = (X, Y, Z, T)
+        return self._parent(coords)
+
+    def diff_addition(P, Q, PQ):
+        """
+        Given the theta points of P, Q and P-Q computes the theta point of
+        P + Q.
+
+        NOTE: Assumes that no coordinate is zero
+
+        Cost: 8S 17M
+        """
+        # Extract out the precomputations
+        Y0, Z0, T0 = P.parent()._arithmetic_precomputation()[-3:]
+
+        # Transform with the Hadamard matrix and multiply
+        # Cost: 8S 7M
+        p1, p2, p3, p4 = P.squared_theta()
+        q1, q2, q3, q4 = Q.squared_theta()
+
+        xp = p1 * q1
+        yp = Y0 * p2 * q2
+        zp = Z0 * p3 * q3
+        tp = T0 * p4 * q4
+
+        # Final coordinates
+        PQx, PQy, PQz, PQt = PQ.coords()
+
+        # Note:
+        # We replace the four divisions by
+        # PQx, PQy, PQz, PQt by 10 multiplications
+        # Cost: 10M
+        PQxy = PQx * PQy
+        PQzt = PQz * PQt
+
+        X, Y, Z, T = P.to_hadamard(xp, yp, zp, tp)
+        X = X * PQzt * PQy
+        Y = Y * PQzt * PQx
+        Z = Z * PQxy * PQt
+        T = T * PQxy * PQz
+
+        coords = (X, Y, Z, T)
+        return P.parent()(coords)
+
+    def scale(self, n):
+        """
+        Scale all coordinates of the ThetaPointDim2 by `n`
+        """
+        x, y, z, t = self.coords()
+        if not isinstance(n, RingElement):
+            raise ValueError(f"Cannot scale by element {n} of type {type(n)}")
+        scaled_coords = (n * x, n * y, n * z, n * t)
+        return self._parent(scaled_coords)
+
+    def double_iter(self, m):
+        """
+        Compute [2^n] Self
+
+        NOTE: Assumes that no coordinate is zero at any point during the doubling
+        """
+        if not isinstance(m, Integer):
+            try:
+                m = Integer(m)
+            except:
+                raise TypeError(f"Cannot coerce input scalar {m = } to an integer")
+
+        if m.is_zero():
+            return self.parent().zero()
+
+        P1 = self
+        for _ in range(m):
+            P1 = P1.double()
+        return P1
+
+    def __mul__(self, m):
+        """
+        Uses Montgomery ladder to compute [m] Self
+
+        NOTE: Assumes that no coordinate is zero at any point during the doubling
+        """
+        # Make sure we're multiplying by something value
+        if not isinstance(m, (int, Integer)):
+            try:
+                m = Integer(m)
+            except:
+                raise TypeError(f"Cannot coerce input scalar {m = } to an integer")
+
+        # If m is zero, return the null point
+        if not m:
+            return self.parent().zero()
+
+        # We are with ±1 identified, so we take the absolute value of m
+        m = abs(m)
+
+        P0, P1 = self, self
+        P2 = P1.double()
+        # If we are multiplying by two, the chain stops here
+        if m == 2:
+            return P2
+
+        # Montgomery double and add.
+        for bit in bin(m)[3:]:
+            Q = P2.diff_addition(P1, P0)
+            if bit == "1":
+                P2 = P2.double()
+                P1 = Q
+            else:
+                P1 = P1.double()
+                P2 = Q
+
+        return P1
+
+    def __rmul__(self, m):
+        return self * m
+
+    def __imul__(self, m):
+        self = self * m
+        return self
+
+    def __eq__(self, other):
+        """
+        Check the quality of two ThetaPoints. Note that as this is a
+        projective equality, we must be careful for when certain coefficients may
+        be zero.
+        """
+        if not isinstance(other, ThetaPointDim2):
+            return False
+
+        a1, b1, c1, d1 = self.coords()
+        a2, b2, c2, d2 = other.coords()
+
+        if d1 != 0 or d2 != 0:
+            return all([a1 * d2 == a2 * d1, b1 * d2 == b2 * d1, c1 * d2 == c2 * d1])
+        elif c1 != 0 or c2 != 0:
+            return all([a1 * c2 == a2 * c1, b1 * c2 == b2 * c1])
+        elif b1 != 0 or b2 != 0:
+            return a1 * b2 == a2 * b1
+        else:
+            return True
+
+    def __repr__(self):
+        return f"Theta point with coordinates: {self.coords()}"