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>

7 files changed, 1354 insertions, 0 deletions

diff --git a/theta_lib/theta_structures/Theta_dim1.py b/theta_lib/theta_structures/Theta_dim1.py
new file mode 100644
index 0000000..e8e94dd
--- /dev/null
+++ b/theta_lib/theta_structures/Theta_dim1.py
@@ -0,0 +1,98 @@
+from sage.all import *
+from ..utilities.supersingular import compute_linearly_independent_point
+from .montgomery_theta import (
+    montgomery_point_to_theta_point,
+    theta_point_to_montgomery_point,
+    torsion_to_theta_null_point,
+)
+
+
+class ThetaStructureDim1:
+	def __init__(self,E,P=None,Q=None):
+		self.E=E
+
+		a_inv=E.a_invariants()
+		
+		A =a_inv[1]
+		if a_inv != (0,A,0,1,0):
+			raise ValueError("The elliptic curve E is not in the Montgomery model.")
+
+		if Q==None:
+			y2=A-2
+			y=y2.sqrt()
+			Q=E([-1,y,1])
+			P=compute_linearly_independent_point(E,Q,4)
+		else:
+			if Q[0]!=-1:
+				raise ValueError("You should enter a canonical 4-torsion basis. Q[0] should be -1.")
+
+		self.P=P
+		self.Q=Q
+
+		self._base_ring=E.base_ring()
+
+		self._point=ThetaPointDim1
+		self._null_point=self._point(self,torsion_to_theta_null_point(P))
+
+	def null_point(self):
+		"""
+		"""
+		return self._null_point
+
+	def base_ring(self):
+		"""
+		"""
+		return self._base_ring
+
+	def zero(self):
+		"""
+		"""
+		return self.null_point()
+
+	def elliptic_curve(self):
+		return self.E
+
+	def torsion_basis(self):
+		return (self.P,self.Q)
+
+	def __call__(self,coords):
+		r"""
+		Input: either a tuple or list of 2 coordinates or an elliptic curve point.
+
+		Output: the corresponding theta point for the self theta structure.
+		"""
+		if isinstance(coords,tuple):
+			return self._point(self,coords)
+		elif isinstance(coords,list):
+			return self._point(self,coords)
+		else:
+			return self._point(self,montgomery_point_to_theta_point(self.null_point().coords(),coords))
+
+	def __repr__(self):
+		return f"Theta structure on {self.elliptic_curve()} with null point: {self.null_point()} induced by the 4-torsion basis {self.torsion_basis()}"
+
+	def to_montgomery_point(self,P):
+		return theta_point_to_montgomery_point(self.E,self.null_point().coords(),P.coords())
+
+
+class ThetaPointDim1:
+    def __init__(self, parent, coords):
+        """
+        """
+        if not isinstance(parent, ThetaStructureDim1):
+            raise ValueError("Entry parent should be a ThetaStructureDim1 object.")
+
+        self._parent = parent
+        self._coords = tuple(coords)
+
+
+    def coords(self):
+    	return self._coords
+
+    def __repr__(self):
+    	return f"Theta point with coordinates: {self.coords()}"
+
+
+
+
+
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()}"
diff --git a/theta_lib/theta_structures/Theta_dim4.py b/theta_lib/theta_structures/Theta_dim4.py
new file mode 100644
index 0000000..7851c56
--- /dev/null
+++ b/theta_lib/theta_structures/Theta_dim4.py
@@ -0,0 +1,351 @@
+from sage.all import *
+from sage.structure.element import get_coercion_model
+
+from .theta_helpers_dim4 import (
+    hadamard,
+    batch_inversion,
+    product_theta_point_dim4,
+    product_to_theta_points_dim4,
+    product_theta_point_dim2_dim4,
+    product_to_theta_points_dim4_dim2,
+    act_point,
+    squared,
+)
+from .Theta_dim1 import ThetaStructureDim1
+from .Tuple_point import TuplePoint
+from ..basis_change.base_change_dim4 import (
+    apply_base_change_theta_dim4,
+    random_symplectic_matrix,
+    base_change_theta_dim4,
+)
+
+cm = get_coercion_model()
+
+
+class ThetaStructureDim4:
+	def __init__(self,null_point,null_point_dual=None,inv_null_point_dual_sq=None):
+		r"""
+		INPUT:
+		- null_point: theta-constants.
+		- inv_null_point_dual_sq: inverse of the squares of dual theta-constants, if provided 
+		(meant to prevent duplicate computation, since this data is already computed when the 
+		codomain of an isogeny is computed).
+		"""
+		if not len(null_point) == 16:
+			raise ValueError("Entry null_point should have 16 coordinates.")
+
+		self._base_ring = cm.common_parent(*(c.parent() for c in null_point))
+		self._point = ThetaPointDim4
+		self._null_point = self._point(self, null_point)
+		self._null_point_dual=null_point_dual
+		self._inv_null_point=None
+		self._inv_null_point_dual_sq=inv_null_point_dual_sq
+
+	def null_point(self):
+		"""
+		"""
+		return self._null_point
+
+	def null_point_dual(self):
+		if self._null_point_dual==None:
+			self._null_point_dual=hadamard(self._null_point.coords())
+		return self._null_point_dual
+
+	def base_ring(self):
+		"""
+		"""
+		return self._base_ring
+
+	def zero(self):
+		"""
+		"""
+		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 __call__(self,coords):
+		return self._point(self,coords)
+
+	def act_null(self,I,J):
+		r"""
+		Point of 2-torsion.
+
+		INPUT:
+		- I, J: two 4-tuples of indices in {0,1}.
+
+		OUTPUT: the action of (I,\chi_J) on the theta null point given by:
+		(I,\chi_J)*P=(\chi_J(I+K)^{-1}P[I+K])_K
+		"""
+		return self.null_point().act_point(I,J)
+
+	def is_K2(self,B):
+		r"""
+		Given a symplectic decomposition A[2]=K_1\oplus K_2 canonically 
+		induced by the theta-null point, determines if B is the canonical
+		basis of K_2 given by act_nul(0,\delta_i)_{1\leq i\leq 4}.
+
+		INPUT:
+		- B: Basis of 4 points of 2-torsion.
+
+		OUTPUT: Boolean True if and only if B is the canonical basis of K_2.
+		"""
+		I0=(0,0,0,0)
+		if B[0]!=self.act_null(I0,(1,0,0,0)):
+			return False
+		if B[1]!=self.act_null(I0,(0,1,0,0)):
+			return False
+		if B[2]!=self.act_null(I0,(0,0,1,0)):
+			return False
+		if B[3]!=self.act_null(I0,(0,0,0,1)):
+			return False
+		return True
+
+	def base_change_struct(self,N):
+		null_coords=self.null_point().coords()
+		new_null_coords=apply_base_change_theta_dim4(N,null_coords)
+		return ThetaStructureDim4(new_null_coords)
+
+	def base_change_coords(self,N,P):
+		coords=P.coords()
+		new_coords=apply_base_change_theta_dim4(N,coords)
+		return self.__call__(new_coords)
+
+	#@cached_method
+	def _arithmetic_precomputation(self):
+		r"""
+		Initializes the precomputation containing the inverse of the theta-constants in standard
+		and dual (Hadamard transformed) theta-coordinates. Assumes no theta-constant is zero.
+		"""
+		O=self.null_point()
+		if all([O[k]!=0 for k in range(16)]):
+			self._inv_null_point=batch_inversion(O.coords())
+			if self._inv_null_point_dual_sq==None:
+				U_chi_0_sq=hadamard(squared(O.coords()))
+				if all([U_chi_0_sq[k]!=0 for k in range(16)]):	
+					self._inv_null_point_dual_sq=batch_inversion(U_chi_0_sq)
+					self._arith_base_change=False
+				else:
+					self._arith_base_change=True
+					self._arithmetic_base_change()
+					#print("Zero dual theta constants.\nRandom symplectic base change is being used for duplication.\nDoublings are more costly than expected.")
+		else:
+			self._arith_base_change=True
+			self._arithmetic_base_change()
+			#print("Zero theta constants.\nRandom symplectic base change is being used for duplication.\nDoublings are more costly than expected.")
+
+	def _arithmetic_base_change(self,max_iter=50):
+		F=self._base_ring
+		if F.degree() == 2:
+			i=self._base_ring.gen()
+		else:
+			assert(F.degree() == 1)
+			Fp2 = GF((F.characteristic(), 2), name='i', modulus=var('x')**2 + 1)
+			i=Fp2.gen()
+
+		count=0
+		O=self.null_point()
+		while count<max_iter:
+			count+=1
+			M=random_symplectic_matrix(4)
+			N=base_change_theta_dim4(M,i)
+
+			NO=apply_base_change_theta_dim4(N,O)
+			NU_chi_0_sq=hadamard(squared(NO))
+
+			if all([NU_chi_0_sq[k]!=0 for k in range(16)]) and all([NO[k]!=0 for k in range(16)]):
+				self._arith_base_change_matrix=N
+				self._arith_base_change_matrix_inv=N.inverse()
+				self._inv_null_point=batch_inversion(NO)
+				self._inv_null_point_dual_sq=batch_inversion(NU_chi_0_sq)
+				break
+
+	def has_suitable_doubling(self):
+		O=self.null_point()
+		UO=hadamard(O.coords())
+		if all([O[k]!=0 for k in range(16)]) and all([UO[k]!=0 for k in range(16)]):
+			return True
+		else:
+			return False
+
+	def hadamard(self):
+		return ThetaStructureDim4(self.null_point_dual())
+
+	def hadamard_change_coords(self,P):
+		new_coords=hadamard(P)
+		return self.__call__(new_coords)
+
+
+class ProductThetaStructureDim1To4(ThetaStructureDim4):
+	def __init__(self,*args):
+		r"""Defines the product theta structure at level 2 of 4 elliptic curves.
+
+		Input: Either
+		- 4 theta structures of dimension 1: T0, T1, T2, T3;
+		- 4 elliptic curves: E0, E1, E2, E3.
+		- 4 elliptic curves E0, E1, E2, E3 and their respective canonical 4-torsion basis B0, B1, B2, B3.
+		"""
+		if len(args)==4:
+			theta_structures=list(args)
+			for k in range(4):
+				if not isinstance(theta_structures[k],ThetaStructureDim1):
+					try:
+						theta_structures[k]=ThetaStructureDim1(theta_structures[k])
+					except:
+						pass
+		elif len(args)==8:
+			theta_structures=[ThetaStructureDim1(args[k],args[4+k][0],args[4+k][1]) for k in range(4)]
+		else:
+			raise ValueError("4 or 8 arguments expected but {} were given.\nYou should enter a list of 4 elliptic curves or ThetaStructureDim1\nor a list of 4 elliptic curves with a 4-torsion basis for each of them.".format(len(args)))
+
+		self._theta_structures=theta_structures
+
+		null_point=product_theta_point_dim4(theta_structures[0].zero().coords(),theta_structures[1].zero().coords(),
+			theta_structures[2].zero().coords(),theta_structures[3].zero().coords())
+
+		ThetaStructureDim4.__init__(self,null_point)
+
+	def product_theta_point(self,theta_points):
+		return self._point(self,product_theta_point_dim4(theta_points[0].coords(),theta_points[1].coords(),
+			theta_points[2].coords(),theta_points[3].coords()))
+
+	def __call__(self,point):
+		if isinstance(point,TuplePoint):
+			theta_points=[]
+			theta_structures=self._theta_structures
+			for i in range(4):
+				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):
+		theta_coords=product_to_theta_points_dim4(P)
+		theta_points=[self._theta_structures[i](theta_coords[i]) for i in range(4)]
+		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(4)]
+		return TuplePoint(montgomery_points)
+
+class ProductThetaStructureDim2To4(ThetaStructureDim4):
+	def __init__(self,theta1,theta2):
+		self._theta_structures=(theta1,theta2)
+
+		null_point=product_theta_point_dim2_dim4(theta1.zero().coords(),theta2.zero().coords())
+
+		ThetaStructureDim4.__init__(self,null_point)
+
+	def product_theta_point(self,P1,P2):
+		return self._point(self,product_theta_point_dim2_dim4(P1.coords(),P2.coords()))
+
+	def __call__(self,point):
+		return self._point(self,point)
+
+	def to_theta_points(self,P):
+		theta_coords=product_to_theta_points_dim4_dim2(P)
+		theta_points=[self._theta_structures[i](theta_coords[i]) for i in range(2)]
+		return theta_points
+
+class ThetaPointDim4:
+    def __init__(self, parent, coords):
+        """
+        """
+        if not isinstance(parent, ThetaStructureDim4):
+            raise ValueError("Entry parent should be a ThetaStructureDim4 object.")
+
+        self._parent = parent
+        self._coords = tuple(coords)
+
+    def parent(self):
+        """
+        """
+        return self._parent
+    
+    def theta(self):
+        """
+        """
+        return self.parent()
+    
+    def coords(self):
+        """
+        """
+        return self._coords
+    
+    def is_zero(self):
+        """
+        """
+        return self == self.parent().zero()
+
+    def __eq__(self, other):
+    	P=self.coords()
+    	Q=other.coords()
+    	
+    	k0=0
+    	while k0<15 and P[k0]==0:
+    		k0+=1
+
+    	for l in range(16):
+    		if P[l]*Q[k0]!=Q[l]*P[k0]:
+    			return False
+    	return True
+
+    def __repr__(self):
+    	return f"Theta point with coordinates: {self.coords()}"
+
+    def __getitem__(self,i):
+    	return self._coords[i]
+
+    def scale(self,lamb):
+    	if lamb==0:
+    		raise ValueError("Entry lamb should be non-zero.")
+
+    	P=self.coords()
+    	return self._parent([lamb*x for x in P])
+
+    def act_point(self,I,J):
+    	r"""
+		Translation by a point of 2-torsion.
+
+		Input:
+		- I, J: two 4-tuples of indices in {0,1}.
+
+		Output: the action of (I,\chi_J) on P given by:
+		(I,\chi_J)*P=(\chi_J(I+K)^{-1}P[I+K])_K
+		"""
+    	return self._parent(act_point(self._coords,I,J))
+
+    def double(self):
+    	## This formula is projective.
+    	## Works only when theta constants are non-zero.
+    	P=self.coords()
+    	if self.parent()._inv_null_point==None or self.parent()._inv_null_point_dual_sq==None:
+    		self.parent()._arithmetic_precomputation()
+    	inv_O,inv_U_chi_0_sq=self.parent()._inv_null_point,self.parent()._inv_null_point_dual_sq
+
+    	if self.parent()._arith_base_change:
+    		P=apply_base_change_theta_dim4(self.parent()._arith_base_change_matrix,P)
+
+    	U_chi_P=squared(hadamard(squared(P)))
+    	for chi in range(16):
+    		U_chi_P[chi]*=inv_U_chi_0_sq[chi]
+
+    	theta_2P = list(hadamard(U_chi_P))
+    	for i in range(16):
+    		theta_2P[i] *= inv_O[i]
+
+    	if self.parent()._arith_base_change:
+    		theta_2P=apply_base_change_theta_dim4(self.parent()._arith_base_change_matrix_inv,theta_2P)
+
+    	return self._parent(theta_2P)
+
+    def double_iter(self,n):
+    	## Computes 2**n*self
+    	Q=self
+    	for i in range(n):
+    		Q=Q.double()
+    	return Q
diff --git a/theta_lib/theta_structures/Tuple_point.py b/theta_lib/theta_structures/Tuple_point.py
new file mode 100644
index 0000000..aac248b
--- /dev/null
+++ b/theta_lib/theta_structures/Tuple_point.py
@@ -0,0 +1,106 @@
+from sage.all import *
+from ..utilities.discrete_log import weil_pairing_pari
+
+class TuplePoint:
+	def __init__(self,*args):
+		if len(args)==1:
+			self._points=list(args[0])
+		else:
+			self._points=list(args)
+
+	def points(self):
+		return self._points
+
+	def parent_curves(self):
+		return [x.curve() for x in self._points]
+
+	def parent_curve(self,i):
+		return self._points[i].curve()
+
+	def n_points(self):
+		return len(self._points)
+
+	def is_zero(self):
+		return all([self._points[i]==0 for i in range(self.n_points())])
+
+	def __repr__(self):
+		return str(self._points)
+
+	def __getitem__(self,i):
+		return self._points[i]
+
+	def __setitem__(self,i,P):
+		self._points[i]=P
+
+	def __eq__(self,other):
+		n_self=self.n_points()
+		n_other=self.n_points()
+		return n_self==n_other and all([self._points[i]==other._points[i] for i in range(n_self)])
+
+	def __add__(self,other):
+		n_self=self.n_points()
+		n_other=self.n_points()
+		
+		if n_self!=n_other:
+			raise ValueError("Cannot add TuplePoint of distinct lengths {} and {}.".format(n_self,n_other))
+
+		points=[]
+		for i in range(n_self):
+			points.append(self._points[i]+other._points[i])
+		return self.__class__(points)
+
+	def __sub__(self,other):
+		n_self=self.n_points()
+		n_other=self.n_points()
+		
+		if n_self!=n_other:
+			raise ValueError("Cannot substract TuplePoint of distinct lengths {} and {}.".format(n_self,n_other))
+
+		points=[]
+		for i in range(n_self):
+			points.append(self._points[i]-other._points[i])
+		return self.__class__(points)
+
+	def __neg__(self):
+		n_self=self.n_points()
+		points=[]
+		for i in range(n_self):
+			points.append(-self._points[i])
+		return self.__class__(points)
+
+	def __mul__(self,m):
+		n_self=self.n_points()
+		points=[]
+		for i in range(n_self):
+			points.append(m*self._points[i])
+		return self.__class__(points)
+
+	def __rmul__(self,m):
+		return self*m
+
+	def double_iter(self,n):
+		result=self
+		for i in range(n):
+			result=2*result
+		return result
+
+	def weil_pairing(self,other,n):
+		n_self=self.n_points()
+		n_other=self.n_points()
+		
+		if n_self!=n_other:
+			raise ValueError("Cannot compute the Weil pairing of TuplePoint of distinct lengths {} and {}.".format(n_self,n_other))
+
+		zeta=1
+		for i in range(n_self):
+			zeta*=weil_pairing_pari(self._points[i],other._points[i],n)
+
+		return zeta
+
+
+
+
+
+
+
+		
diff --git a/theta_lib/theta_structures/montgomery_theta.py b/theta_lib/theta_structures/montgomery_theta.py
new file mode 100644
index 0000000..6dddb2b
--- /dev/null
+++ b/theta_lib/theta_structures/montgomery_theta.py
@@ -0,0 +1,68 @@
+from sage.all import *
+
+def torsion_to_theta_null_point(P):
+	r=P[0]
+	s=P[2]
+	return (r+s,r-s)
+
+def montgomery_point_to_theta_point(O,P):
+	if P[0]==P[2]==0:
+		return O
+	else:
+		a,b=O
+		return (a*(P[0]-P[2]),b*(P[0]+P[2]))
+
+def theta_point_to_montgomery_point(E,O,P,twist=False):
+	a,b=O
+
+	x=a*P[1]+b*P[0]
+	z=a*P[1]-b*P[0]
+
+	if twist:
+		x=-x
+
+	if z==0:
+		return E(0)
+	else:
+		x=x/z
+
+		a_inv=E.a_invariants()
+		
+		A =a_inv[1]
+		if a_inv != (0,A,0,1,0):
+			raise ValueError("The elliptic curve E is not in the Montgomery model.")
+
+		y2=x*(x**2+A*x+1)
+		if not is_square(y2):
+			raise ValueError("The Montgomery point is not on the base field.")
+		else:
+			y=y2.sqrt()
+			return E([x,y,1])
+
+def lift_kummer_montgomery_point(E,x,z=1):
+	if z==0:
+		return E(0)
+	elif z!=1:
+		x=x/z
+
+	a_inv=E.a_invariants()
+		
+	A =a_inv[1]
+	if a_inv != (0,A,0,1,0):
+		raise ValueError("The elliptic curve E is not in the Montgomery model.")
+
+	y2=x*(x**2+A*x+1)
+	if not is_square(y2):
+		raise ValueError("The Montgomery point is not on the base field.")
+	else:
+		y=y2.sqrt()
+		return E([x,y,1])
+
+def null_point_to_montgomery_coeff(a,b):
+	return -2*(a**4+b**4)/(a**4-b**4)
+
+
+
+
+
+
diff --git a/theta_lib/theta_structures/theta_helpers_dim2.py b/theta_lib/theta_structures/theta_helpers_dim2.py
new file mode 100644
index 0000000..a57789d
--- /dev/null
+++ b/theta_lib/theta_structures/theta_helpers_dim2.py
@@ -0,0 +1,32 @@
+from sage.all import *
+
+def batch_inversion(L):
+	r"""Does n inversions in 3(n-1)M+1I.
+
+	Input:
+	- L: list of elements to invert.
+
+	Output:
+	- [1/x for x in L]
+	"""
+	# Given L=[a0,...,an]
+	# Computes multiples=[a0, a0.a1, ..., a0...an]
+	multiples=[L[0]]
+	for ai in L[1:]:
+		multiples.append(multiples[-1]*ai)
+
+	# Computes inverses=[1/(a0...an),...,1/a0]
+	inverses=[1/multiples[-1]]
+	for i in range(1,len(L)):
+		inverses.append(inverses[-1]*L[-i])
+
+	# Finally computes [1/a0,...,1/an]
+	result=[inverses[-1]]
+	for i in range(2,len(L)+1):
+		result.append(inverses[-i]*multiples[i-2])
+	return result
+
+def product_theta_point(theta1,theta2):
+    a,b=theta1
+    c,d=theta2
+    return [a*c,b*c,a*d,b*d]
diff --git a/theta_lib/theta_structures/theta_helpers_dim4.py b/theta_lib/theta_structures/theta_helpers_dim4.py
new file mode 100644
index 0000000..29599b0
--- /dev/null
+++ b/theta_lib/theta_structures/theta_helpers_dim4.py
@@ -0,0 +1,259 @@
+from sage.all import *
+import itertools
+
+## Index management
+@cached_function
+def index_to_multindex(k):
+	r"""
+	Input: 
+	- k: integer between 0 and 15.
+
+	Output: binary decomposition of k.
+	"""
+	L_ind=[]
+	l=k
+	for i in range(4):
+		L_ind.append(l%2)
+		l=l//2
+	return tuple(L_ind)
+
+@cached_function
+def multindex_to_index(*args):
+	r"""
+	Input: 4 elements i0,i1,i2,i3 in {0,1}.
+
+	Output: k=i0+2*i1+4*i2+8*i3.
+	"""
+	if len(args)==4:
+		i0,i1,i2,i3=args
+	else:
+		i0,i1,i2,i3=args[0]
+	return i0+2*i1+4*i2+8*i3
+
+@cached_function
+def scal_prod(i,j):
+	r"""
+	Input: Two integers i and j in {0,...,15}.
+
+	Output: Scalar product of the bits of i and j mod 2.
+	"""
+	return (int(i)&int(j)).bit_count()%2
+
+def act_point(P,I,J):
+	r"""
+	Input:
+	- P: a point with 16 coordinates.
+	- I, J: two 4-tuples of indices in {0,1}.
+
+	Output: the action of (I,\chi_J) on P given by:
+	(I,\chi_J)*P=(\chi_J(I+K)^{-1}P[I+K])_K
+	"""
+	Q=[]
+	i=multindex_to_index(I)
+	j=multindex_to_index(J)
+	for k in range(16):
+		ipk=i^k
+		Q.append((-1)**scal_prod(ipk,j)*P[ipk])
+	return Q
+
+## Product of theta points
+def product_theta_point_dim4(P0,P1,P2,P3):
+	# Computes the product theta coordinates of a product of 4 elliptic curves.
+	P=[0 for k in range(16)]
+	for i0, i1, i2, i3 in itertools.product([0,1],repeat=4):
+		P[multindex_to_index(i0,i1,i2,i3)]=P0[i0]*P1[i1]*P2[i2]*P3[i3]
+	return P
+
+def product_theta_point_dim2_dim4(P0,P1):
+	# Computes the product theta coordinates of a product of 2 abelian surfaces.
+	P=[0 for k in range(16)]
+	for i0, i1, i2, i3 in itertools.product([0,1],repeat=4):
+		P[multindex_to_index(i0,i1,i2,i3)]=P0[i0+2*i1]*P1[i2+2*i3]
+	return P
+
+## 4-dimensional product Theta point to 1-dimensional Theta points
+def product_to_theta_points_dim4(P):
+	Fp2=P[0].parent()
+	d_Pi={0:None,1:None,2:None,3:None}
+	d_index_ratios={0:None,1:None,2:None,3:None}# Index of numertors/denominators to compute the theta points.
+	for k in range(4):
+		is_zero=True# theta_1(Pk)=0
+		for I in itertools.product([0,1],repeat=3):
+			J=list(I)
+			J.insert(k,1)
+			j=multindex_to_index(*J)
+			if P[j]!=0:
+				is_zero=False
+				d_index_ratios[k]=[j^(2**k),j]
+				break
+		if is_zero:
+			d_Pi[k]=(Fp2(1),Fp2(0))
+	L_num=[]
+	L_denom=[]
+	d_index_num_denom={}
+	for k in range(4):
+		if d_Pi[k]==None:# Point has non-zero coordinate theta_1
+			d_index_num_denom[k]=len(L_num)
+			L_num.append(P[d_index_ratios[k][0]])
+			L_denom.append(P[d_index_ratios[k][1]])
+	L_denom=batch_inversion(L_denom)
+	for k in range(4):
+		if d_Pi[k]==None:
+			d_Pi[k]=(L_num[d_index_num_denom[k]]*L_denom[d_index_num_denom[k]],Fp2(1))
+	return d_Pi[0],d_Pi[1],d_Pi[2],d_Pi[3]
+
+## 4-dimensional product Theta point to 2-dimensional Theta points
+def product_to_theta_points_dim4_dim2(P):
+	Fp2=P[0].parent()
+	k0=0# Index of the denominator k0=multindex_to_index(I0,J0) and
+	# we divide by theta_{I0,J0}=theta_{I0}*theta_{J0}!=0
+	while k0<=15 and P[k0]==0:
+		k0+=1
+	i0, j0 = k0%4, k0//4
+	inv_theta_k0=1/P[k0]
+	theta_P1=[]
+	theta_P2=[]
+	for i in range(4):
+		if i==i0:
+			theta_P1.append(1)
+		else:
+			theta_P1.append(P[i+4*j0]*inv_theta_k0)
+	for j in range(4):
+		if j==j0:
+			theta_P2.append(1)
+		else:
+			theta_P2.append(P[i0+4*j]*inv_theta_k0)
+	return theta_P1, theta_P2
+
+
+
+
+
+## Usual theta transformations
+@cached_function
+def inv_16(F):
+	return 1/F(16)
+
+def hadamard2(x,y):
+    return (x+y, x-y)
+
+def hadamard4(x,y,z,t):
+    x,y=hadamard2(x,y)
+    z,t=hadamard2(z,t)
+    return (x+z, y+t, x-z, y-t)
+
+def hadamard8(a,b,c,d,e,f,g,h):
+    a,b,c,d=hadamard4(a,b,c,d)
+    e,f,g,h=hadamard4(e,f,g,h)
+    return (a+e, b+f, c+g, d+h, a-e, b-f, c-g, d-h)
+
+def hadamard16(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p):
+    a,b,c,d,e,f,g,h=hadamard8(a,b,c,d,e,f,g,h)
+    i,j,k,l,m,n,o,p=hadamard8(i,j,k,l,m,n,o,p)
+    return (a+i, b+j, c+k, d+l, e+m, f+n, g+o, h+p, a-i, b-j, c-k, d-l, e-m, f-n, g-o, h-p)
+
+def hadamard_inline(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p):
+    a,b=a+b,a-b
+    c,d=c+d,c-d
+    e,f=e+f,e-f
+    g,h=g+h,g-h
+    i,j=i+j,i-j
+    k,l=k+l,k-l
+    m,n=m+n,m-n
+    o,p=o+p,o-p
+    a,b,c,d=a+c,b+d,a-c,b-d
+    e,f,g,h=e+g,f+h,e-g,f-h
+    i,j,k,l=i+k,j+l,i-k,j-l
+    m,n,o,p=m+o,n+p,m-o,n-p
+    a,b,c,d,e,f,g,h=a+e, b+f, c+g, d+h, a-e, b-f, c-g, d-h
+    i,j,k,l,m,n,o,p=i+m,j+n,k+o,l+p,i-m,j-n,k-o,l-p
+    return (a+i, b+j, c+k, d+l, e+m, f+n, g+o, h+p, a-i, b-j, c-k, d-l, e-m, f-n, g-o, h-p)
+
+def hadamard(P):
+    return hadamard16(*P)
+    #return hadamard_inline(*P)
+
+def hadamard_inverse(P):
+	H_inv_P=[]
+	C=inv_16(P[0].parent())
+	for j in range(16):
+		HP.append(0)
+		for k in range(16):
+			if scal_prod(k,j)==0:
+				H_inv_P[j]+=P[k]
+			else:
+				H_inv_P[j]-=P[k]
+		H_inv_P[j]=H_inv_P[j]*C
+	return H_inv_P
+
+def squared(P):
+	return [x**2 for x in P]
+
+def batch_inversion(L):
+	r"""Does n inversions in 3(n-1)M+1I.
+
+	Input:
+	- L: list of elements to invert.
+
+	Output:
+	- [1/x for x in L]
+	"""
+	# Given L=[a0,...,an]
+	# Computes multiples=[a0, a0.a1, ..., a0...an]
+	multiples=[L[0]]
+	for ai in L[1:]:
+		multiples.append(multiples[-1]*ai)
+
+	# Computes inverses=[1/(a0...an),...,1/a0]
+	inverses=[1/multiples[-1]]
+	for i in range(1,len(L)):
+		inverses.append(inverses[-1]*L[-i])
+
+	# Finally computes [1/a0,...,1/an]
+	result=[inverses[-1]]
+	for i in range(2,len(L)+1):
+		result.append(inverses[-i]*multiples[i-2])
+	return result
+
+
+## Functions to handle zero theta coordinates
+def find_zeros(P):
+	L_ind_zeros=[]
+	for i in range(16):
+		if P[i]==0:
+			L_ind_zeros.append(i)
+	return L_ind_zeros
+
+def find_translates(L_ind_zeros):
+	L_ind_non_zero=[]
+	L_ind_origin=L_ind_zeros.copy()
+	
+	for i in range(16):
+		if i not in L_ind_zeros:
+			L_ind_non_zero.append(i)
+
+	L_trans=[]
+	while L_ind_origin!=[]:
+		n_target_max=0
+		ind_trans_max=0
+		for k in range(16):
+			trans=[x^k for x in L_ind_origin]
+			n_target=0
+			for y in trans:
+				if y in L_ind_non_zero:
+					n_target+=1
+			if n_target>n_target_max:
+				n_target_max=n_target
+				ind_trans_max=k
+		L_trans.append(ind_trans_max)
+		L_ind_remove=[]
+		for x in L_ind_origin:
+			if x^ind_trans_max in L_ind_non_zero:
+				L_ind_remove.append(x)
+		for x in L_ind_remove:
+			L_ind_origin.remove(x)
+	return L_trans
+
+
+
+