Clean up PEGASIS submodule inclusion

7 files changed, 0 insertions, 1354 deletions

diff --git a/theta_lib/theta_structures/Theta_dim1.py b/theta_lib/theta_structures/Theta_dim1.py
deleted file mode 100644
index e8e94dd..0000000
--- a/theta_lib/theta_structures/Theta_dim1.py
+++ /dev/null
@@ -1,98 +0,0 @@
-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
deleted file mode 100644
index 0d9955a..0000000
--- a/theta_lib/theta_structures/Theta_dim2.py
+++ /dev/null
@@ -1,440 +0,0 @@
-# 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
deleted file mode 100644
index 7851c56..0000000
--- a/theta_lib/theta_structures/Theta_dim4.py
+++ /dev/null
@@ -1,351 +0,0 @@
-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
deleted file mode 100644
index aac248b..0000000
--- a/theta_lib/theta_structures/Tuple_point.py
+++ /dev/null
@@ -1,106 +0,0 @@
-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
deleted file mode 100644
index 6dddb2b..0000000
--- a/theta_lib/theta_structures/montgomery_theta.py
+++ /dev/null
@@ -1,68 +0,0 @@
-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
deleted file mode 100644
index a57789d..0000000
--- a/theta_lib/theta_structures/theta_helpers_dim2.py
+++ /dev/null
@@ -1,32 +0,0 @@
-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
deleted file mode 100644
index 29599b0..0000000
--- a/theta_lib/theta_structures/theta_helpers_dim4.py
+++ /dev/null
@@ -1,259 +0,0 @@
-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
-
-
-
-