Clean up PEGASIS submodule inclusion

31 files changed, 24 insertions, 5826 deletions

diff --git a/.gitmodules b/.gitmodules
new file mode 100644
index 0000000..cfe8997
--- /dev/null
+++ b/.gitmodules
@@ -0,0 +1,3 @@
+[submodule "Theta_dim4"]
+	path = Theta_dim4
+	url = git@github.com:Pierrick-Dartois/Theta_dim4.git
diff --git a/LICENSE b/LICENSE
index 4ad0dcc..9ffc20c 100644
--- a/LICENSE
+++ b/LICENSE
@@ -12,7 +12,8 @@ the Software, and to permit persons to whom the Software is furnished to do so,
 subject to the following conditions:
 
 The above copyright notice and this permission notice shall be included in all
-copies or substantial portions of the Software.
+copies or substantial portions of the Software along with relevant license and 
+copyright information from the Theta_dim4 submodule.
 
 THE SOFTWARE IS PROVIDED “AS IS”, WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
 IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY, FITNESS
diff --git a/README.md b/README.md
index 29bb845..f316781 100644
--- a/README.md
+++ b/README.md
@@ -95,7 +95,7 @@ logging.getLogger("lcsidh").setLevel(logging.INFO)
 
 ## Project structure
 
-- `theta_lib`: code for computing with 4D isogenies;
+- `Theta_dim4`: code for computing with 4D isogenies imported as a submodule from <https://github.com/Pierrick-Dartois/Theta_dim4>
 - `pegasis.py`: contains `PEGASIS`, the main class running the algorithm
 - `coin.py`: utilities to solve the norm equation with u and v sums of squares
 - `const_precomp.py`: precomputed constants for trial division
@@ -106,3 +106,17 @@ logging.getLogger("lcsidh").setLevel(logging.INFO)
 - `uv_params.py`: contains the parameters for solving the norm equation
   ("finding uv")
 - `test_vectors.txt`: contains test vectors for verification
+
+## License
+
+See [LICENSE](LICENSE)
+
+Copyright 2025
+
+Pierrick Dartois, Jonathan Komada Eriksen, Tako Boris Fouotsa, Arthur Herlédan
+Le Merdy, Riccardo Invernizzi, Damien Robert, Ryan Rueger, Frederik Vercauteren,
+Benjamin Wesolowski
+
+Third party code is used in this library: 
+
+- `Theta_dim4`; Apache 2.0: "Copyright (c) 2025 Pierrick Dartois"
diff --git a/Theta_dim4 b/Theta_dim4
new file mode 160000
+Subproject cc67feef0152c13a45af185d6666ad1404bc787
diff --git a/pegasis.py b/pegasis.py
index 6664e0d..38384e6 100644
--- a/pegasis.py
+++ b/pegasis.py
@@ -9,7 +9,8 @@ from uv_params import UV_params
 from ideals import ideal_to_sage, PrincipalGenerator, Conjugate, RandomDegreeOnePrimeIdeal
 from elkies import Elkies
 
-from theta_lib.isogenies.Kani_clapoti import KaniClapotiIsog
+#from theta_lib.isogenies.Kani_clapoti import KaniClapotiIsog
+from Theta_dim4.Theta_dim4_sage.pkg.isogenies.Kani_clapoti import KaniClapotiIsog
 
 logger = logging.getLogger(__name__)
 logger.setLevel(logging.WARNING)
diff --git a/theta_lib/basis_change/base_change_dim2.py b/theta_lib/basis_change/base_change_dim2.py
deleted file mode 100644
index 4994529..0000000
--- a/theta_lib/basis_change/base_change_dim2.py
+++ /dev/null
@@ -1,150 +0,0 @@
-from sage.all import *
-
-import itertools
-
-def complete_symplectic_matrix_dim2(C, D, n=4):
-	Zn = Integers(n)
-
-	# Compute first row
-	F = D.stack(-C).transpose()
-	xi = F.solve_right(vector(Zn, [1, 0]))
-
-	# Rearrange TODO: why?
-	v = vector(Zn, [xi[2], xi[3], -xi[0], -xi[1]])
-
-	# Compute second row
-	F2 = F.stack(v)
-	yi = F2.solve_right(vector(Zn, [0, 1, 0]))
-	M = Matrix(Zn,
-		[
-			[xi[0], yi[0], C[0, 0], C[0, 1]],
-			[xi[1], yi[1], C[1, 0], C[1, 1]],
-			[xi[2], yi[2], D[0, 0], D[0, 1]],
-			[xi[3], yi[3], D[1, 0], D[1, 1]],
-		],
-	)
-
-	return M
-
-def block_decomposition(M):
-	"""
-	Given a 4x4 matrix, return the four 2x2 matrices as blocks
-	"""
-	A = M.matrix_from_rows_and_columns([0, 1], [0, 1])
-	B = M.matrix_from_rows_and_columns([2, 3], [0, 1])
-	C = M.matrix_from_rows_and_columns([0, 1], [2, 3])
-	D = M.matrix_from_rows_and_columns([2, 3], [2, 3])
-	return A, B, C, D
-
-def is_symplectic_matrix_dim2(M):
-	A, B, C, D = block_decomposition(M)
-	if B.transpose() * A != A.transpose() * B:
-		return False
-	if C.transpose() * D != D.transpose() * C:
-		return False
-	if A.transpose() * D - B.transpose() * C != identity_matrix(2):
-		return False
-	return True
-
-def base_change_theta_dim2(M, zeta):
-    r"""
-    Computes the matrix N (in row convention) of the new level 2
-    theta coordinates after a symplectic base change of A[4] given
-    by M\in Sp_4(Z/4Z). We have:
-    (theta'_i)=N*(theta_i)
-    where the theta'_i are the new theta coordinates and theta_i,
-    the theta coordinates induced by the product theta structure.
-    N depends on the fourth root of unity zeta=e_4(Si,Ti) (where
-    (S0,S1,T0,T1) is a symplectic basis if A[4]).
-
-    Inputs:
-    - M: symplectic base change matrix.
-    - zeta: a primitive 4-th root of unity induced by the Weil-pairings
-    of the symplectic basis of A[4].
-
-    Output: Matrix N of base change of theta-coordinates.
-    """
-    # Split 4x4 matrix into 2x2 blocks
-    A, B, C, D = block_decomposition(M)
-
-    # Initialise N to 4x4 zero matrix
-    N = [[0 for _ in range(4)] for _ in range(4)]
-
-    def choose_non_vanishing_index(C, D, zeta):
-        """
-        Choice of reference non-vanishing index (ir0,ir1)
-        """
-        for ir0, ir1 in itertools.product([0, 1], repeat=2):
-            L = [0, 0, 0, 0]
-            for j0, j1 in itertools.product([0, 1], repeat=2):
-                k0 = C[0, 0] * j0 + C[0, 1] * j1
-                k1 = C[1, 0] * j0 + C[1, 1] * j1
-
-                l0 = D[0, 0] * j0 + D[0, 1] * j1
-                l1 = D[1, 0] * j0 + D[1, 1] * j1
-
-                e = -(k0 + 2 * ir0) * l0 - (k1 + 2 * ir1) * l1
-                L[ZZ(k0 + ir0) % 2 + 2 * (ZZ(k1 + ir1) % 2)] += zeta ** (ZZ(e))
-
-            # Search if any L value in L is not zero
-            if any([x != 0 for x in L]):
-                return ir0, ir1
-
-    ir0, ir1 = choose_non_vanishing_index(C, D, zeta)
-
-    for i0, i1, j0, j1 in itertools.product([0, 1], repeat=4):
-        k0 = A[0, 0] * i0 + A[0, 1] * i1 + C[0, 0] * j0 + C[0, 1] * j1
-        k1 = A[1, 0] * i0 + A[1, 1] * i1 + C[1, 0] * j0 + C[1, 1] * j1
-
-        l0 = B[0, 0] * i0 + B[0, 1] * i1 + D[0, 0] * j0 + D[0, 1] * j1
-        l1 = B[1, 0] * i0 + B[1, 1] * i1 + D[1, 0] * j0 + D[1, 1] * j1
-
-        e = i0 * j0 + i1 * j1 - (k0 + 2 * ir0) * l0 - (k1 + 2 * ir1) * l1
-        N[i0 + 2 * i1][ZZ(k0 + ir0) % 2 + 2 * (ZZ(k1 + ir1) % 2)] += zeta ** (ZZ(e))
-
-    return Matrix(N)
-
-def montgomery_to_theta_matrix_dim2(zero12,N=identity_matrix(4),return_null_point=False):
-	r"""
-	Computes the matrix that transforms Montgomery coordinates on E1*E2 into theta coordinates 
-	with respect to a certain theta structure given by a base change matrix N.
-
-	Input: 
-	- zero12: theta-null point for the product theta structure on E1*E2[2]:
-	zero12=[zero1[0]*zero2[0],zero1[1]*zero2[0],zero1[0]*zero2[1],zero1[1]*zero2[1]],
-	where the zeroi are the theta-null points of Ei for i=1,2.
-	- N: base change matrix from the product theta-structure.
-	- return_null_point: True if the theta-null point obtained after applying N to zero12 is returned.
-
-	Output:
-	- Matrix for the change of coordinates:
-	(X1*X2,Z1*X2,X1*Z2,Z1*Z2)-->(theta_00,theta_10,theta_01,theta_11).
-	- If return_null_point, the theta-null point obtained after applying N to zero12 is returned.
-	"""
-
-	Fp2=zero12[0].parent()
-
-	M=zero_matrix(Fp2,4,4)
-
-	for i in range(4):
-		for j in range(4):
-			M[i,j]=N[i,j]*zero12[j]
-
-	M2=[]
-	for i in range(4):
-		M2.append([M[i,0]+M[i,1]+M[i,2]+M[i,3],-M[i,0]-M[i,1]+M[i,2]+M[i,3],-M[i,0]+M[i,1]-M[i,2]+M[i,3],M[i,0]-M[i,1]-M[i,2]+M[i,3]])
-
-	if return_null_point:
-		null_point=[M[i,0]+M[i,1]+M[i,2]+M[i,3] for i in range(4)]
-		return Matrix(M2), null_point
-	else:
-		return Matrix(M2)
-
-def apply_base_change_theta_dim2(N,P):
-	Q=[]
-	for i in range(4):
-		Q.append(0)
-		for j in range(4):
-			Q[i]+=N[i,j]*P[j]
-	return Q
-
diff --git a/theta_lib/basis_change/base_change_dim4.py b/theta_lib/basis_change/base_change_dim4.py
deleted file mode 100644
index aaf685f..0000000
--- a/theta_lib/basis_change/base_change_dim4.py
+++ /dev/null
@@ -1,152 +0,0 @@
-from sage.all import *
-import itertools
-
-from ..theta_structures.theta_helpers_dim4 import multindex_to_index
-
-def bloc_decomposition(M):
-	I1=[0,1,2,3]
-	I2=[4,5,6,7]
-
-	A=M[I1,I1]
-	B=M[I2,I1]
-	C=M[I1,I2]
-	D=M[I2,I2]
-
-	return A,B,C,D
-
-def mat_prod_vect(A,I):
-	J=[]
-	for i in range(4):
-		J.append(0)
-		for k in range(4):
-			J[i]+=A[i,k]*I[k]
-	return tuple(J)
-
-def add_tuple(I,J):
-	K=[]
-	for k in range(4):
-		K.append(I[k]+J[k])
-	return tuple(K)
-
-def scal_prod_tuple(I,J):
-	s=0
-	for k in range(4):
-		s+=I[k]*J[k]
-	return s
-
-def red_mod_2(I):
-	J=[]
-	for x in I:
-		J.append(ZZ(x)%2)
-	return tuple(J)
-
-def choose_non_vanishing_index(C,D,zeta):
-	for I0 in itertools.product([0,1],repeat=4):
-		L=[0 for k in range(16)]
-		for J in itertools.product([0,1],repeat=4):
-			CJ=mat_prod_vect(C,J)
-			DJ=mat_prod_vect(D,J)
-			e=-scal_prod_tuple(CJ,DJ)-2*scal_prod_tuple(I0,DJ)
-
-			I0pDJ=add_tuple(I0,CJ)
-
-			L[multindex_to_index(red_mod_2(I0pDJ))]+=zeta**(ZZ(e))
-		for k in range(16):
-			if L[k]!=0:
-				return I0,L
-
-
-def base_change_theta_dim4(M,zeta):
-	
-	Z4=Integers(4)
-	A,B,C,D=bloc_decomposition(M.change_ring(Z4))
-
-	I0,L0=choose_non_vanishing_index(C,D,zeta)
-
-
-	N=[L0]+[[0 for j in range(16)] for i in range(15)]
-	for I in itertools.product([0,1],repeat=4):
-		if I!=(0,0,0,0):
-			AI=mat_prod_vect(A,I)
-			BI=mat_prod_vect(B,I)
-			for J in itertools.product([0,1],repeat=4):
-				CJ=mat_prod_vect(C,J)
-				DJ=mat_prod_vect(D,J)
-
-				AIpCJ=add_tuple(AI,CJ)
-				BIpDJ=add_tuple(BI,DJ)
-
-				e=scal_prod_tuple(I,J)-scal_prod_tuple(AIpCJ,BIpDJ)-2*scal_prod_tuple(I0,BIpDJ)
-				N[multindex_to_index(I)][multindex_to_index(red_mod_2(add_tuple(AIpCJ,I0)))]+=zeta**(ZZ(e))
-
-	Fp2=zeta.parent()
-	return matrix(Fp2,N)
-
-def apply_base_change_theta_dim4(N,P):
-	Q=[]
-	for i in range(16):
-		Q.append(0)
-		for j in range(16):
-			Q[i]+=N[i,j]*P[j]
-	return Q
-
-def complete_symplectic_matrix_dim4(C,D,n=4):
-	Zn=Integers(n)
-
-	Col_I4=[matrix(Zn,[[1],[0],[0],[0]]),matrix(Zn,[[0],[1],[0],[0],[0]]),
-	matrix(Zn,[[0],[0],[1],[0],[0],[0]]),matrix(Zn,[[0],[0],[0],[1],[0],[0],[0]])]
-
-	L_DC=block_matrix([[D.transpose(),-C.transpose()]])
-	Col_AB_i=L_DC.solve_right(Col_I4[0])
-	A_t=Col_AB_i[[0,1,2,3],0].transpose()
-	B_t=Col_AB_i[[4,5,6,7],0].transpose()
-
-	for i in range(1,4):
-		F=block_matrix(2,1,[L_DC,block_matrix(1,2,[B_t,-A_t])])
-		Col_AB_i=F.solve_right(Col_I4[i])
-		A_t=block_matrix(2,1,[A_t,Col_AB_i[[0,1,2,3],0].transpose()])
-		B_t=block_matrix(2,1,[B_t,Col_AB_i[[4,5,6,7],0].transpose()])
-
-	A=A_t.transpose()
-	B=B_t.transpose()
-
-	M=block_matrix([[A,C],[B,D]])
-
-	return M
-
-def random_symmetric_matrix(n,r):
-	Zn=Integers(n)
-	M=zero_matrix(Zn,r,r)
-	for i in range(r):
-		M[i,i]=randint(0,n-1)
-	for i in range(r):
-		for j in range(i+1,r):
-			M[i,j]=randint(0,n-1)
-			M[j,i]=M[i,j]
-	return M
-
-
-def random_symplectic_matrix(n):
-	Zn=Integers(n)
-
-	C=random_symmetric_matrix(n,4)
-	Cp=random_symmetric_matrix(n,4)
-
-	A=random_matrix(Zn,4,4)
-	while not A.is_invertible():
-		A=random_matrix(Zn,4,4)
-
-	tA_inv=A.inverse().transpose()
-	ACp=A*Cp
-	return block_matrix([[ACp,A],[C*ACp-tA_inv,C*A]])
-
-def is_symplectic_matrix_dim4(M):
-	A,B,C,D=bloc_decomposition(M)
-	if B.transpose()*A!=A.transpose()*B:
-		return False
-	if C.transpose()*D!=D.transpose()*C:
-		return False
-	if A.transpose()*D-B.transpose()*C!=identity_matrix(4):
-		return False
-	return True
-
diff --git a/theta_lib/basis_change/canonical_basis_dim1.py b/theta_lib/basis_change/canonical_basis_dim1.py
deleted file mode 100644
index e1c3d1f..0000000
--- a/theta_lib/basis_change/canonical_basis_dim1.py
+++ /dev/null
@@ -1,76 +0,0 @@
-from sage.all import *
-from ..utilities.discrete_log import weil_pairing_pari, discrete_log_pari
-
-def last_four_torsion(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.")
-	y2=A-2
-	y=y2.sqrt()
-	return E([-1,y,1])
-
-
-def make_canonical(P,Q,A,preserve_pairing=False):
-	r"""
-	Input:
-	- P,Q: a basis of E[A].
-	- A: an integer divisible by 4.
-	- preserve_pairing: boolean indicating if we want to preserve pairing at level 4.
-	
-	Output:
-	- P1,Q1: basis of E[A].
-	- U1,U2: basis of E[4] induced by (P1,Q1) ((A//4)*P1=U1, (A//4)*Q1=U2) such that U2[0]=-1 
-	and e_4(U1,U2)=i if not preserve_pairing and e_4(U1,U2)=e_4((A//4)*P,(A//4)*Q) if preserve_pairing.
-	- M: base change matrix (in row convention) from (P1,Q1) to (P,Q).
-	
-	We say that (U1,U2) is canonical and that (P1,Q1) induces or lies above a canonical basis. 
-	"""
-	E=P.curve()
-	Fp2=E.base_ring()
-	i=Fp2.gen()
-
-	assert i**2==-1
-
-	T2=last_four_torsion(E)
-	V1=(A//4)*P
-	V2=(A//4)*Q
-	U1=V1
-	U2=V2
-	
-	a1=discrete_log_pari(weil_pairing_pari(U1,T2,4),i,4)
-	b1=discrete_log_pari(weil_pairing_pari(U2,T2,4),i,4)
-
-	if a1%2!=0:
-		c1=inverse_mod(a1,4)
-		d1=c1*b1
-		P1=P
-		Q1=Q-d1*P
-		U1,U2=U1,U2-d1*U1
-		M=matrix(ZZ,[[1,0],[d1,1]])
-	else:
-		c1=inverse_mod(b1,4)
-		d1=c1*a1
-		P1=Q
-		Q1=P-d1*Q
-		U1,U2=U2,U1-d1*U2
-		M=matrix(ZZ,[[d1,1],[1,0]])
-
-	if preserve_pairing:
-		e4=weil_pairing_pari(V1,V2,4)
-	else:
-		e4=i
-	
-	if weil_pairing_pari(U1,U2,4)!=e4:
-		U2=-U2
-		Q1=-Q1
-		M[0,1]=-M[0,1]
-		M[1,1]=-M[1,1]
-	
-	assert (A//4)*P1==U1
-	assert (A//4)*Q1==U2
-	assert weil_pairing_pari(U1,U2,4)==e4
-	assert M[0,0]*P1+M[0,1]*Q1==P
-	assert M[1,0]*P1+M[1,1]*Q1==Q
-	
-	return P1,Q1,U1,U2,M
diff --git a/theta_lib/basis_change/kani_base_change.py b/theta_lib/basis_change/kani_base_change.py
deleted file mode 100644
index e6de2e2..0000000
--- a/theta_lib/basis_change/kani_base_change.py
+++ /dev/null
@@ -1,975 +0,0 @@
-from sage.all import *
-from ..basis_change.canonical_basis_dim1 import make_canonical
-from ..basis_change.base_change_dim2 import is_symplectic_matrix_dim2
-from ..basis_change.base_change_dim4 import (
-    complete_symplectic_matrix_dim4,
-    is_symplectic_matrix_dim4,
-    bloc_decomposition,
-)
-from ..theta_structures.Tuple_point import TuplePoint
-
-
-def base_change_canonical_dim2(P1,P2,R1,R2,q,f):
-	r"""
-
-	Input:
-	- P1, P2: basis of E1[2**f].
-	- R1, R2: images of P1, P2 by \sigma: E1 --> E2.
-	- q: degree of \sigma.
-	- f: log_2(order of P1 and P2).
-
-	Output:
-	- P1_doubles: list of 2**i*P1 for i in {0,...,f-2}.
-	- P2_doubles: list of 2**i*P2 for i in {0,...,f-2}.
-	- R1_doubles: list of 2**i*R1 for i in {0,...,f-2}.
-	- R2_doubles: list of 2**i*R2 for i in {0,...,f-2},
-	- T1, T2: canonical basis of E1[4].
-	- U1, U2: canonical basis of E2[4].
-	- M0: base change matrix of the symplectic basis 2**(f-2)*B1 of E1*E2[4] given by:
-	B1:=[[(P1,0),(0,R1)],[(P2,0),(0,lamb*R2)]]
-	where lamb is the modular inverse of q mod 2**f, so that:
-	e_{2**f}(P1,P2)=e_{2**f}(R1,lamb*R2).
-	in the canonical symplectic basis:
-	B0:=[[(T1,0),(0,U1)],[(T2,0),(0,U2)]].
-	"""
-	lamb=inverse_mod(q,4)
-
-	P1_doubles=[P1]
-	P2_doubles=[P2]
-	R1_doubles=[R1]
-	R2_doubles=[R2]
-
-	for i in range(f-2):
-		P1_doubles.append(2*P1_doubles[-1])
-		P2_doubles.append(2*P2_doubles[-1])
-		R1_doubles.append(2*R1_doubles[-1])
-		R2_doubles.append(2*R2_doubles[-1])
-
-	# Constructing canonical basis of E1[4] and E2[4].
-	_,_,T1,T2,MT=make_canonical(P1_doubles[-1],P2_doubles[-1],4,preserve_pairing=True)
-	_,_,U1,U2,MU=make_canonical(R1_doubles[-1],lamb*R2_doubles[-1],4,preserve_pairing=True)
-
-	Z4=Integers(4)
-	M0=matrix(Z4,[[MT[0,0],0,MT[1,0],0],
-		[0,MU[0,0],0,MU[1,0]],
-		[MT[0,1],0,MT[1,1],0],
-		[0,MU[0,1],0,MU[1,1]]])
-
-	return P1_doubles,P2_doubles,R1_doubles,R2_doubles,T1,T2,U1,U2,M0
-
-def gluing_base_change_matrix_dim2(a1,a2,q):
-	r"""Computes the symplectic base change matrix of a symplectic basis (*,B_K4) of E1*E2[4]
-	given by the kernel of the dimension 2 gluing isogeny:
-	B_K4=2**(f-2)[([a1]P1-[a2]P2,R1),([a1]P2+[a2]P1,R2)]
-	in the basis $2**(f-2)*B1$ given by:
-	B1:=[[(P1,0),(0,R1)],[(P2,0),(0,lamb*R2)]]
-	where:
-	- lamb is the inverse of q modulo 2**f.
-	- (P1,P2) is the canonical basis of E1[2**f].
-	- (R1,R2) is the image of (P1,P2) by sigma.
-
-	Input:
-	- a1, q: integers.
-
-	Output:
-	- M: symplectic base change matrix of (*,B_K4) in 2**(f-2)*B1.
-	"""
-
-	Z4=Integers(4)
-
-	mu=inverse_mod(a1,4)
-
-	A=matrix(Z4,[[0,mu],
-		[0,0]])
-	B=matrix(Z4,[[0,0],
-		[-1,-ZZ(mu*a2)]])
-
-	C=matrix(Z4,[[ZZ(a1),ZZ(a2)],
-		[1,0]])
-	D=matrix(Z4,[[-ZZ(a2),ZZ(a1)],
-		[0,ZZ(q)]])
-
-	#M=complete_symplectic_matrix_dim2(C, D, 4)
-	M=block_matrix([[A,C],[B,D]])
-
-	assert is_symplectic_matrix_dim2(M)
-
-	return M
-
-# ============================================== #
-#     Functions for the class KaniClapotiIsog    #
-# ============================================== #
-
-def clapoti_cob_matrix_dim2(integers):
-	gu,xu,yu,gv,xv,yv,Nbk,Nck,e,m = integers
-
-	xu = ZZ(xu)
-	xv = ZZ(xv)
-	Nbk = ZZ(Nbk)
-	Nck = ZZ(Nck)
-	u = ZZ(gu*(xu**2+yu**2))
-	v = ZZ(gv*(xv**2+yv**2))
-	mu = inverse_mod(u,4)
-	suv = xu*xv+yu*yv
-	inv_Nbk = inverse_mod(Nbk,4)
-	inv_gugvNcksuv = inverse_mod(gu*gv*Nck*suv,4)
-
-	Z4=Integers(4)
-
-	M=matrix(Z4,[[0,0,u*Nbk,0],
-		[0,inv_Nbk*inv_gugvNcksuv,gu*suv,0],
-		[-inv_Nbk*mu,0,0,gu*Nbk*u],
-		[0,0,0,gu*gv*Nbk*Nck*suv]])
-
-	assert is_symplectic_matrix_dim2(M)
-
-	return M
-
-def clapoti_cob_matrix_dim2_dim4(integers):
-	gu,xu,yu,gv,xv,yv,Nbk,Nck,e,m = integers
-
-	xu = ZZ(xu)
-	yu = ZZ(yu)
-	xv = ZZ(xv)
-	yv = ZZ(yv)
-	gu = ZZ(gu)
-	gv = ZZ(gv)
-	Nbk = ZZ(Nbk)
-	Nck = ZZ(Nck)
-	u = ZZ(gu*(xu**2+yu**2))
-	v = ZZ(gv*(xv**2+yv**2))
-	suv = xu*xv+yu*yv
-	duv = xv*yu-xu*yv
-	duv_2m = duv//2**m
-	mu = inverse_mod(u,4)
-	nu = inverse_mod(v,4)
-	sigmauv = inverse_mod(suv,4)
-	inv_guNbk = inverse_mod(gu*Nbk,4)
-	lamb = nu*gu*gv*Nbk*suv
-	mu1 = ZZ(mu*gu**2*gv*suv*Nbk*Nck*duv_2m)
-	mu2 = ZZ(duv_2m*gu*sigmauv*(Nbk*u*yu+gv*xv*Nck*duv))
-	mu3 = ZZ(duv_2m*gu*sigmauv*(Nbk*u*xu-gv*yv*Nck*duv))
-
-	Z4=Integers(4)
-
-	M=matrix(Z4,[[gu*xu,-gu*yu,0,0,0,0,mu2,mu3],
-		[0,0,lamb*xv,-lamb*yv,mu1*yu,mu1*xu,0,0],
-		[gu*yu,gu*xu,0,0,0,0,-mu3,mu2],
-		[0,0,lamb*yv,lamb*xv,-mu1*xu,mu1*yu,0,0],
-		[0,0,0,0,mu*xu,-mu*yu,0,0],
-		[0,0,0,0,0,0,inv_guNbk*xv*sigmauv,-inv_guNbk*yv*sigmauv],
-		[0,0,0,0,mu*yu,mu*xu,0,0],
-		[0,0,0,0,0,0,inv_guNbk*yv*sigmauv,inv_guNbk*xv*sigmauv]])
-
-	assert is_symplectic_matrix_dim4(M)
-
-	return M
-
-def clapoti_cob_splitting_matrix(integers):
-	gu,xu,yu,gv,xv,yv,Nbk,Nck,e,m = integers
-
-	v=ZZ(gv*(xv**2+yv**2))
-	vNck=ZZ(v*Nck)
-	inv_vNck=inverse_mod(vNck,4)
-
-	Z4=Integers(4)
-
-	M=matrix(Z4,[[0,0,0,0,-1,0,0,0],
-		[0,0,0,0,0,-1,0,0],
-		[0,0,vNck,0,0,0,0,0],
-		[0,0,0,vNck,0,0,0,0],
-		[1,0,-vNck,0,0,0,0,0],
-		[0,1,0,-vNck,0,0,0,0],
-		[0,0,0,0,1,0,inv_vNck,0],
-		[0,0,0,0,0,1,0,inv_vNck]])
-
-	assert is_symplectic_matrix_dim4(M)
-
-	return M
-
-# =============================================== #
-#     Functions for the class KaniFixedDegDim2    #
-# =============================================== #
-
-def fixed_deg_gluing_matrix_Phi1(u,a,b,c,d):
-	u,a,b,c,d = ZZ(u),ZZ(a),ZZ(b),ZZ(c),ZZ(d)
-
-	mu = inverse_mod(u,4)
-	inv_cmd = inverse_mod(c-d,4)
-
-	Z4 = Integers(4)
-
-	M = matrix(Z4,[[0,0,u,0],
-		[0,inv_cmd,c+d,0],
-		[-mu,0,0,(d**2-c**2)*mu],
-		[0,0,0,c-d]])
-
-	assert is_symplectic_matrix_dim2(M)
-
-	return M
-
-def fixed_deg_gluing_matrix_Phi2(u,a,b,c,d):
-	u,a,b,c,d = ZZ(u),ZZ(a),ZZ(b),ZZ(c),ZZ(d)
-
-	mu = inverse_mod(u,4)
-	inv_cpd = inverse_mod(c+d,4)
-
-	Z4 = Integers(4)
-
-	M = matrix(Z4,[[0,0,u,0],
-		[0,-inv_cpd,d-c,0],
-		[-mu,0,0,(d**2-c**2)*mu],
-		[0,0,0,-(c+d)]])
-
-	assert is_symplectic_matrix_dim2(M)
-
-	return M
-
-def fixed_deg_gluing_matrix_dim4(u,a,b,c,d,m):
-	u,a,b,c,d = ZZ(u),ZZ(a),ZZ(b),ZZ(c),ZZ(d)
-
-	mu = inverse_mod(u,4)
-	nu = ZZ((-mu**2)%4)
-	amb_2m = ZZ((a-b)//2**m)
-	apb_2m = ZZ((a+b)//2**m)
-	u2pc2md2_2m = ZZ((u**2+c**2-d**2)//2**m)
-	inv_cmd = inverse_mod(c-d,4)
-	inv_cpd = inverse_mod(c+d,4)
-
-
-	Z4 = Integers(4)
-
-	M = matrix(Z4,[[1,0,0,0,0,0,-u2pc2md2_2m,-apb_2m*(c+d)],
-		[0,0,(c+d)*(c-d)*nu,(a-b)*(c-d)*nu,0,amb_2m*(c-d),0,0],
-		[0,1,0,0,0,0,amb_2m*(c-d),-u2pc2md2_2m],
-		[0,0,-(a+b)*(c+d)*nu,(c+d)*(c-d)*nu,-apb_2m*(c+d),0,0,0],
-		[0,0,0,0,1,0,0,0],
-		[0,0,0,0,0,0,1,(a+b)*inv_cmd],
-		[0,0,0,0,0,1,0,0],
-		[0,0,0,0,0,0,(b-a)*inv_cpd,1]])
-
-	assert is_symplectic_matrix_dim4(M)
-
-	return M
-
-def fixed_deg_gluing_matrix(u,a,b,c,d):
-	r"""
-	Deprecated.
-	"""
-
-	mu = inverse_mod(u,4)
-	nu = (-mu**2)%4
-
-	Z4 = Integers(4)
-
-	M = matrix(Z4,[[0,0,0,0,ZZ(u),0,0,0],
-		[0,0,0,0,0,ZZ(u),0,0],
-		[0,0,ZZ(nu*(a+b)),ZZ(nu*(d-c)),ZZ(a+b),ZZ(d-c),0,0],
-		[0,0,ZZ(nu*(c+d)),ZZ(nu*(a-b)),ZZ(c+d),ZZ(a-b),0,0],
-		[ZZ(-mu),0,0,0,0,0,ZZ(u),0],
-		[0,ZZ(-mu),0,0,0,0,0,ZZ(u)],
-		[0,0,0,0,0,0,ZZ(a-b),ZZ(-c-d)],
-		[0,0,0,0,0,0,ZZ(c-d),ZZ(a+b)]])
-
-	assert is_symplectic_matrix_dim4(M)
-
-	return M
-
-def fixed_deg_splitting_matrix(u):
-
-	mu = inverse_mod(u,4)
-
-	Z4 = Integers(4)
-
-	M = matrix(Z4,[[0,0,0,0,-1,0,0,0],
-		[0,0,0,0,0,-1,0,0],
-		[0,0,ZZ(-u),0,0,0,0,0],
-		[0,0,0,ZZ(-u),0,0,0,0],
-		[1,0,ZZ(-mu),0,0,0,0,0],
-		[0,1,0,ZZ(-mu),0,0,0,0],
-		[0,0,0,0,ZZ(mu),0,ZZ(mu),0],
-		[0,0,0,0,0,ZZ(mu),0,ZZ(mu)]])
-
-	assert is_symplectic_matrix_dim4(M)
-
-	return M
-	
-
-# ========================================================== #
-#     Functions for the class KaniEndo (one isogeny chain)   #
-# ========================================================== #
-
-def gluing_base_change_matrix_dim2_dim4(a1,a2,m,mua2):
-	r"""Computes the symplectic base change matrix of a symplectic basis (*,B_K4) of Am*Am[4]
-	given by the kernel of the dimension 4 gluing isogeny Am*Am-->B:
-	
-	B_K4=2**(e-m)[(Phi([a1]P1,sigma(P1)),Phi([a2]P1,0)),(Phi([a1]P2,sigma(P2)),Phi([a2]P2,0)),
-	(Phi(-[a2]P1,0),Phi([a1]P1,sigma(P1))),(Phi(-[a2]P2,0),Phi([a2]P2,sigma(P2)))]
-	
-	in the basis associated to the product theta-structure of level 2 of Am*Am:
-	
-	B:=[(S1,0),(S2,0),(0,S1),(0,S2),(T1,0),(T2,0),(0,T1),(0,T2)]
-	
-	where:
-	- (P1,P2) is the canonical basis of E1[2**f].
-	- (R1,R2) is the image of (P1,P2) by sigma.
-	- Phi is the 2**m-isogeny E1*E2-->Am (m first steps of the chain in dimension 2).
-	- S1=[2**e]Phi([lamb]P2,[a]sigma(P1)+[b]sigma(P2)).
-	- S2=[2**e]Phi([mu]P1,[c]sigma(P1)+[d]sigma(P2)).
-	- T1=[2**(e-m)]Phi([a1]P1-[a2]P2,sigma(P1)).
-	- T2=[2**(e-m)]Phi([a1]P2+[a2]P1,sigma(P2)).
-	- (S1,S2,T1,T2) is induced by the image by Phi of a symplectic basis of E1*E2[2**(m+2)] lying
-	above the symplectic basis of E1*E2[4] outputted by gluing_base_change_matrix_dim2.
-
-	INPUT:
-	- a1, a2: integers.
-	- m: integer (number of steps in dimension 2).
-	- mua2: product mu*a2.
-
-	OUTPUT:
-	- M: symplectic base change matrix of (*,B_K4) in B.
-	"""
-	a1a2_2m=ZZ(a1*a2//2**m)
-	a22_2m=ZZ(a2**2//2**m)
-
-	Z4=Integers(4)
-
-	C=matrix(Z4,[[-a1a2_2m,a22_2m,a22_2m,a1a2_2m],
-		[-a22_2m,-a1a2_2m,-a1a2_2m,a22_2m],
-		[-a22_2m,-a1a2_2m,-a1a2_2m,a22_2m],
-		[a1a2_2m,-a22_2m,-a22_2m,-a1a2_2m]])
-
-	D=matrix(Z4,[[1,0,0,0],
-		[mua2,1,0,-mua2],
-		[0,0,1,0],
-		[0,mua2,mua2,1]])
-
-	M=complete_symplectic_matrix_dim4(C,D,4)
-
-	assert is_symplectic_matrix_dim4(M)
-
-	return M
-
-def splitting_base_change_matrix_dim4(a1,a2,q,m,M0,A_B,mu=None):
-	r"""
-	Let F be the endomorphism of E1^2*E2^2 given by Kani's lemma. Write: 
-			E1^2*E2^2 -- Phi x Phi --> Am^2 -- G --> E1^2*E2^2,
-	where Phi: E1 x E1 --> Am is a 2**m-isogeny in dimension 2. 
-	Let (U_1,...,U_4,V_1,...,V_4) be a symplectic basis of Am^2[2**(e-m+2)] 
-	such that V_i=Phi x Phi(W_i), where W_1,...,W_4 have order 2**(e+2), lie over ker(F) 
-	and generate an isotropic subgroup:
-	W_1=([a1]P1-[2^e/a1]P1,[a2]P1,R2,0)
-	W_2=([a1]Q1,[a2]Q1,S2,0)
-	W_3=(-[a2]P1,[a1]P1,0,R2)
-	W_4=(-[a2]Q1,[a1]Q1,0,S2),
-	with (P1,Q1), a basis of E1[2**(e+2)] and (R2,S2) its image via 
-	sigma: E1 --> E2. Then B:=([2^(e-m)]G(U_1),...,[2^(e-m)]G(U_4),G(V_1),...,G(V_4))
-	is a symplectic basis of E1^2*E2^2[4].
-
-	We assume that ([2^(e-m)]U_1,...,[2^(e-m)]U_4) is the symplectic complement of 
-	([2^(e-m)]V_1,...,[2^(e-m)]V_4) that has been outputted by 
-	gluing_base_change_matrix_dim2_dim4 for the gluing isogeny on Am^2 
-	(first 2-isogeny of G). This function computes the base change matrix of B
-	in the symplectic basis of E1^2*E2^2[4]:
-	B0=[(T1,0,0,0),(0,T1,0,0),(0,0,T2,0),(0,0,0,T2),(U1,0,0,0),(0,U1,0,0),
-	(0,0,U2,0),(0,0,0,U2)]
-	associated to the product Theta structure on E1^2*E2^2.
-
-	INPUT:
-	- a1,a2,q: integers defining F (q=deg(sigma)).
-	- m: 2-adic valuation of a2.
-	- M0: base change matrix of the symplectic basis 2**e*B1 of E1*E2[4] 
-	given by:
-	B1:=[[(P1,0),(0,R2)],[(Q1,0),(0,lamb*S2)]]
-	in the canonical symplectic basis:
-	B0:=[[(T1,0),(0,T2)],[(U1,0),(0,U2)]],
-	where lamb is the modular inverse of q mod 2**(e+2), so that:
-	e_{2**(e+2)}(P1,P2)=e_{2**(e+2)}(R1,lamb*R2).
-	- A_B: 4 first columns (left part) of the symplectic matrix outputted by
-	gluing_base_change_matrix_dim2_dim4.
-	- mu, a, b, c, d: integers defining the product Theta structure of Am^2
-	given by the four torsion basis [2**(e-m)]*B1 of Am, where:
-	B1=[[2**m]Phi([2**(m+1)]P2,[a]sigma(P1)+[b]sigma(P2)),
-	[2**m]Phi([mu]P1,[2**(m+1)]sigma(P1)+[d]sigma(P2)),
-	Phi([a1]P1-[a2]P2,sigma(P1)),
-	Phi([a1]P2+[a2]P1,sigma(P2))].
-	Only mu is given.
-
-	OUTPUT: The desired base change matrix.
-	"""
-	Z4=Integers(4)
-
-	a2_2m=ZZ(a2//2**m)
-	a12_q_2m=ZZ((a1**2+q)//2**m)
-
-	inv_q=inverse_mod(q,4)
-	inv_a1=inverse_mod(a1,4)
-
-	lamb=ZZ(2**(m+1))
-	if mu==None:
-		mu=ZZ((1-2**(m+1)*q)*inv_a1)
-	a=ZZ(2**(m+1)*a2*inv_q)
-	b=ZZ(-(1+2**(m+1)*a1)*inv_q)
-	c=ZZ(2**(m+1))
-	d=ZZ(-mu*a2*inv_q)
-
-	# Matrix of the four torsion basis of E1^2*E2^2[4] given by
-	# ([2^(e-m)]G(B1[0],0),[2^(e-m)]G(B1[1],0),[2^(e-m)]G(0,B1[0]),[2^(e-m)]G(0,B1[1]),
-	# G(B1[2],0),G(B1[3],0),G(0,B1[2]),G(0,B1[3])) in the basis induced by 
-	# [2**e](P1,Q1,R2,[1/q]S2)
-	M1=matrix(Z4,[[a*q,mu*a1+c*q,0,mu*a2,a12_q_2m,a1*a2_2m,a1*a2_2m,a2*a2_2m],
-		[0,-mu*a2,a*q,mu*a1+c*q,-a1*a2_2m,-a2*a2_2m,a12_q_2m,a1*a2_2m],
-		[a1*a,a1*c-mu,-a*a2,-c*a2,0,-a2_2m,-a2_2m,0],
-		[a2*a,a2*c,a*a1,c*a1-mu,a2_2m,0,0,-a2_2m],
-		[lamb*a1+b*q,d*q,lamb*a2,0,-a1*a2_2m,a12_q_2m,-a2*a2_2m,a1*a2_2m],
-		[-lamb*a2,0,lamb*a1+b*q,d*q,a2*a2_2m,-a1*a2_2m,-a1*a2_2m,a12_q_2m],
-		[(a1*b-lamb)*q,a1*d*q,-b*a2*q,-a2*d*q,a2_2m*q,0,0,-a2_2m*q],
-		[a2*b*q,a2*d*q,(b*a1-lamb)*q,a1*d*q,0,a2_2m*q,a2_2m*q,0]])
-	#A,B,C,D=bloc_decomposition(M1)
-	#if B.transpose()*A!=A.transpose()*B:
-		#print("B^T*A!=A^T*B")
-	#if C.transpose()*D!=D.transpose()*C:
-		#print("C^T*D!=D^T*C")
-	#if A.transpose()*D-B.transpose()*C!=identity_matrix(4):
-		#print(A.transpose()*D-B.transpose()*C)
-		#print("A^T*D-B^T*C!=I")
-	#print(M1)
-	#print(M1.det())
-
-	# Matrix of ([2^e]G(U_1),...,[2^e]G(U_4)) in the basis induced by
-	# [2**e](P1,Q1,R2,[1/q]S2)
-	M_left=M1*A_B
-	#print(A_B)
-	#print(M_left)
-
-	# Matrix of (G(V_1),...,G(V_4)) in the basis induced by [2**e](P1,Q1,R2,[1/q]S2)
-	M_right=matrix(Z4,[[0,0,0,0],
-		[a2*inv_a1,0,1,0],
-		[inv_a1,0,0,0],
-		[0,0,0,0],
-		[0,1,0,-a2*inv_a1],
-		[0,0,0,0],
-		[0,0,0,0],
-		[0,0,0,q*inv_a1]])
-
-	# Matrix of the basis induced by [2**e](P1,Q1,R2,[1/q]S2) in the basis 
-	# B0 (induced by T1, U1, T2, U2)
-	MM0=matrix(Z4,[[M0[0,0],0,M0[0,1],0,M0[0,2],0,M0[0,3],0],
-				  [0,M0[0,0],0,M0[0,1],0,M0[0,2],0,M0[0,3]],
-				  [M0[1,0],0,M0[1,1],0,M0[1,2],0,M0[1,3],0],
-				  [0,M0[1,0],0,M0[1,1],0,M0[1,2],0,M0[1,3]],
-				  [M0[2,0],0,M0[2,1],0,M0[2,2],0,M0[2,3],0],
-				  [0,M0[2,0],0,M0[2,1],0,M0[2,2],0,M0[2,3]],
-				  [M0[3,0],0,M0[3,1],0,M0[3,2],0,M0[3,3],0],
-				  [0,M0[3,0],0,M0[3,1],0,M0[3,2],0,M0[3,3]]])
-
-	M=MM0*block_matrix(1,2,[M_left,M_right])
-
-	#A,B,C,D=bloc_decomposition(M)
-
-	#M=complete_symplectic_matrix_dim4(C,D)
-
-	#print(M.det())
-	#print(M)
-
-	A,B,C,D=bloc_decomposition(M)
-	if B.transpose()*A!=A.transpose()*B:
-		print("B^T*A!=A^T*B")
-	if C.transpose()*D!=D.transpose()*C:
-		print("C^T*D!=D^T*C")
-	if A.transpose()*D-B.transpose()*C!=identity_matrix(4):
-		print("A^T*D-B^T*C!=I")
-	assert is_symplectic_matrix_dim4(M)
-
-	return M
-
-# ============================================================================ #
-#     Functions for the class KaniEndoHalf (isogeny chain decomposed in two)   #
-# ============================================================================ #
-
-def complete_kernel_matrix_F1(a1,a2,q,f):
-	r"""Computes the symplectic base change matrix of a symplectic basis of the form (*,B_Kp1) 
-	in the symplectic basis of E1^2*E2^2[2**f] given by:
-	B1:=[[(P1,0,0,0),(0,P1,0,0),(0,0,R1,0),(0,0,0,R1)],
-	[(P2,0,0,0),(0,P2,0,0),(0,0,lamb*R2,0),(0,0,0,lamb*R2)]]
-	where:
-	- B_Kp1 is a basis of an isotropic subgroup of E1^2*E2^2[2**f] lying above ker(F1).
-	By convention B_Kp1=[(\tilde{\alpha}_1(P1,0),\Sigma(P1,0)),
-	(\tilde{\alpha}_1(P2,0),\Sigma(P2,0)),
-	(\tilde{\alpha}_1(0,P1),\Sigma(0,P1)),
-	(\tilde{\alpha}_1(0,P2),\Sigma(0,P2))]
-	- lamb is the inverse of q modulo 2**f.
-	- (P1,P2) is the canonical basis of E1[2**f].
-	- (R1,R2) is the image of (P1,P2) by sigma.
-
-	Input:
-	- a1, a2, q: Integers such that q+a1**2+a2**2=2**e.
-	- f: integer determining the accessible 2-torsion in E1 (E1[2**f]).
-
-	Output:
-	- M: symplectic base change matrix of (*,B_Kp1) in B1.
-	"""
-	N=2**f
-	ZN=Integers(N)
-
-	C=matrix(ZN,[[a1,0,-a2,0],
-		[a2,0,a1,0],
-		[1,0,0,0],
-		[0,0,1,0]])
-
-	D=matrix(ZN,[[0,a1,0,-a2],
-		[0,a2,0,a1],
-		[0,q,0,0],
-		[0,0,0,q]])
-
-	assert C.transpose()*D==D.transpose()*C
-
-	M=complete_symplectic_matrix_dim4(C,D,N)
-
-	assert is_symplectic_matrix_dim4(M)
-
-	return M
-
-def complete_kernel_matrix_F2_dual(a1,a2,q,f):
-	r"""Computes the symplectic base change matrix of a symplectic basis of the form (*,B_Kp2) 
-	in the symplectic basis of E1^2*E2^2[2**f] given by:
-	B1:=[[(P1,0,0,0),(0,P1,0,0),(0,0,R1,0),(0,0,0,R1)],
-	[(P2,0,0,0),(0,P2,0,0),(0,0,lamb*R2,0),(0,0,0,lamb*R2)]]
-	where:
-	- B_Kp2 is a basis of an isotropic subgroup of E1^2*E2^2[2**f] lying above ker(\tilde{F2}).
-	By convention B_Kp2=[(\alpha_1(P1,0),-\Sigma(P1,0)),
-	(\alpha_1(P2,0),-\Sigma(P2,0)),
-	(\alpha_1(0,P1),-\Sigma(0,P1)),
-	(\alpha_1(0,P2),-\Sigma(0,P2))]. 
-	- lamb is the inverse of q modulo 2**f.
-	- (P1,P2) is the canonical basis of E1[2**f].
-	- (R1,R2) is the image of (P1,P2) by sigma.
-
-	Input:
-	- a1, a2, q: Integers such that q+a1**2+a2**2=2**e.
-	- f: integer determining the accessible 2-torsion in E1 (E1[2**f]).
-
-	Output:
-	- M: symplectic base change matrix of (*,B_Kp2) in B1.
-	"""
-	N=2**f
-	ZN=Integers(N)
-
-	C=matrix(ZN,[[a1,0,a2,0],
-		[-a2,0,a1,0],
-		[-1,0,0,0],
-		[0,0,-1,0]])
-
-	D=matrix(ZN,[[0,a1,0,a2],
-		[0,-a2,0,a1],
-		[0,-q,0,0],
-		[0,0,0,-q]])
-
-
-
-	M=complete_symplectic_matrix_dim4(C,D,N)
-
-	assert is_symplectic_matrix_dim4(M)
-
-	return M
-
-def matrix_F_dual(a1,a2,q,f):
-	r""" Computes the matrix of \tilde{F}(B1) in B1, where:
-	B1:=[[(P1,0,0,0),(0,P1,0,0),(0,0,R1,0),(0,0,0,R1)],
-	[(P2,0,0,0),(0,P2,0,0),(0,0,lamb*R2,0),(0,0,0,lamb*R2)]]
-	as defined in complete_kernel_matrix_F2_dual.
-
-	Input:
-	- a1, a2, q: Integers such that q+a1**2+a2**2=2**e.
-	- f: integer determining the accessible 2-torsion in E1 (E1[2**f]).
-
-	Output:
-	- M: symplectic base change matrix of \tilde{F}(B1) in B1.
-	"""
-	N=2**f
-	ZN=Integers(N)
-
-	M=matrix(ZN,[[a1,-a2,-q,0,0,0,0,0],
-		[a2,a1,0,-q,0,0,0,0],
-		[1,0,a1,a2,0,0,0,0],
-		[0,1,-a2,a1,0,0,0,0],
-		[0,0,0,0,a1,-a2,-1,0],
-		[0,0,0,0,a2,a1,0,-1],
-		[0,0,0,0,q,0,a1,a2],
-		[0,0,0,0,0,q,-a2,a1]])
-
-	return M
-
-def matrix_F(a1,a2,q,f):
-	r""" Computes the matrix of F(B1) in B1, where:
-	B1:=[[(P1,0,0,0),(0,P1,0,0),(0,0,R1,0),(0,0,0,R1)],
-	[(P2,0,0,0),(0,P2,0,0),(0,0,lamb*R2,0),(0,0,0,lamb*R2)]]
-	as defined in complete_kernel_matrix_F1.
-
-	Input:
-	- a1, a2, q: Integers such that q+a1**2+a2**2=2**e.
-	- f: integer determining the accessible 2-torsion in E1 (E1[2**f]).
-
-	Output:
-	- M: symplectic base change matrix of \tilde{F}(B1) in B1.
-	"""
-	N=2**f
-	ZN=Integers(N)
-
-	M=matrix(ZN,[[a1,a2,q,0,0,0,0,0],
-		[-a2,a1,0,q,0,0,0,0],
-		[-1,0,a1,-a2,0,0,0,0],
-		[0,-1,a2,a1,0,0,0,0],
-		[0,0,0,0,a1,a2,1,0],
-		[0,0,0,0,-a2,a1,0,1],
-		[0,0,0,0,-q,0,a1,-a2],
-		[0,0,0,0,0,-q,a2,a1]])
-
-	return M
-
-def starting_two_symplectic_matrices(a1,a2,q,f):
-	r"""
-	Computes the matrices of two symplectic basis of E1^2*E2^2[2**f] given 
-	by (*,B_Kp1) and (*,B_Kp2) in the basis
-	B1:=[[(P1,0,0,0),(0,P1,0,0),(0,0,R1,0),(0,0,0,R1)],
-	[(P2,0,0,0),(0,P2,0,0),(0,0,lamb*R2,0),(0,0,0,lamb*R2)]]
-	as defined in complete_kernel_matrix_F1.
-
-	Input:
-	- a1, a2, q: Integers such that q+a1**2+a2**2=2**e.
-	- f: integer determining the accessible 2-torsion in E1 (E1[2**f]).
-
-	Output:
-	- M1, M2: the symplectic base change matrices of (*,B_Kp1) and (*,B_Kp2) in B1.
-	"""
-	M1_0=complete_kernel_matrix_F1(a1,a2,q,f)
-	MatF=matrix_F(a1,a2,q,f)
-
-	# Matrix of an isotropic subgroup of E1^2*E2^2[2**f] lying above ker(\tilde{F2}).
-	Block_right2=MatF*M1_0[:,[0,1,2,3]]
-
-	N=ZZ(2**f)
-
-	C=Block_right2[[0,1,2,3],:]
-	D=Block_right2[[4,5,6,7],:]
-
-	assert C.transpose()*D==D.transpose()*C
-
-	# Matrix of the resulting symplectic basis (*,B_Kp2)
-	M2=complete_symplectic_matrix_dim4(C,D,N)
-
-	MatF_dual=matrix_F_dual(a1,a2,q,f)
-
-	Block_right1=MatF_dual*M2[:,[0,1,2,3]]
-
-	C=Block_right1[[0,1,2,3],:]
-	D=Block_right1[[4,5,6,7],:]
-
-	A=M1_0[[0,1,2,3],[0,1,2,3]]
-	B=M1_0[[4,5,6,7],[0,1,2,3]]
-
-	assert C.transpose()*D==D.transpose()*C
-	assert B.transpose()*A==A.transpose()*B
-
-	# Matrix of the resulting symplectic basis (*,B_Kp1)
-	M1=block_matrix(1,2,[M1_0[:,[0,1,2,3]],-Block_right1])
-
-	assert is_symplectic_matrix_dim4(M1)
-
-	A,B,C,D=bloc_decomposition(M1)
-	a2_div=a2
-	m=0
-	while a2_div%2==0:
-		m+=1
-		a2_div=a2_div//2
-	for j in range(4):
-		assert (-D[0,j]*a1-C[0,j]*a2-D[2,j])%2**m==0
-		assert (C[0,j]*a1-D[0,j]*a2+C[2,j]*q)%2**m==0
-		assert (-D[1,j]*a1-C[1,j]*a2-D[3,j])%2**m==0
-		assert (C[1,j]*a1-D[1,j]*a2+C[3,j]*q)%2**m==0
-
-	return M1, M2
-
-def gluing_base_change_matrix_dim2_F1(a1,a2,q):
-	r"""Computes the symplectic base change matrix of a symplectic basis (*,B_K4) of E1*E2[4]
-	given by the kernel of the dimension 2 gluing isogeny:
-	B_K4=2**(f-2)[([a1]P1-[a2]P2,R1),([a1]P2+[a2]P1,R2)]
-	in the basis $2**(f-2)*B1$ given by:
-	B1:=[[(P1,0),(0,R1)],[(P2,0),(0,[1/q]*R2)]]
-	where:
-	- lamb is the inverse of q modulo 2**f.
-	- (P1,P2) is the canonical basis of E1[2**f].
-	- (R1,R2) is the image of (P1,P2) by sigma.
-
-	Input:
-	- a1, q: integers.
-
-	Output:
-	- M: symplectic base change matrix of (*,B_K4) in 2**(f-2)*B1.
-	"""
-	return gluing_base_change_matrix_dim2(a1,a2,q)
-
-def gluing_base_change_matrix_dim2_dim4_F1(a1,a2,q,m,M1):
-	r"""Computes the symplectic base change matrix of the symplectic basis Bp of Am*Am[4] induced
-	by the image of the symplectic basis (x_1, ..., x_4, y_1, ..., y_4) of E1^2*E2^2[2**(e1+2)]
-	adapted to ker(F1)=[4]<y_1, ..., y_4> in the basis associated to the product theta-structure 
-	of level 2 of Am*Am:
-	
-	B:=[(S1,0),(S2,0),(0,S1),(0,S2),(T1,0),(T2,0),(0,T1),(0,T2)]
-	
-	where:
-	- (P1,Q1) is the canonical basis of E1[2**f].
-	- (R2,S2) is the image of (P1,P2) by sigma.
-	- Phi is the 2**m-isogeny E1*E2-->Am (m first steps of the chain in dimension 2).
-	- S1=[2**e]Phi([lamb]Q1,[a]sigma(P1)+[b]sigma(Q1)).
-	- S2=[2**e]Phi([mu]P1,[c]sigma(P1)+[d]sigma(Q1)).
-	- T1=[2**(e-m)]Phi([a1]P1-[a2]Q1,sigma(P1)).
-	- T2=[2**(e-m)]Phi([a1]Q1+[a2]P1,sigma(Q1)).
-	- (S1,S2,T1,T2) is induced by the image by Phi of a symplectic basis of E1*E2[2**(m+2)] lying
-	above the symplectic basis of E1*E2[4] outputted by gluing_base_change_matrix_dim2.
-
-	INPUT:
-	- a1, a2, q: integers.
-	- m: integer (number of steps in dimension 2 and 2-adic valuation of a2).
-	- M1: matrix of (x_1, ..., x_4, y_1, ..., y_4) in the symplectic basis of E1^2*E2^2[2**(e1+2)]
-	given by:
-
-	B1:=[[(P1,0,0,0),(0,P1,0,0),(0,0,R2,0),(0,0,0,R2)],
-	[(Q1,0,0,0),(0,Q1,0,0),(0,0,[1/q]*S2,0),(0,0,0,[1/q]*S2)]]
-
-	OUTPUT:
-	- M: symplectic base change matrix of Bp in B.
-	"""
-
-	inv_a1=inverse_mod(a1,2**(m+2))
-	inv_q=inverse_mod(q,2**(m+2))
-	lamb=ZZ(2**(m+1))
-	mu=ZZ((1-2**(m+1)*q)*inv_a1)
-	a=ZZ(2**(m+1)*a2*inv_q)
-	bq=ZZ((-1-2**(m+1)*a1))
-	c=ZZ(2**(m+1))
-	dq=-ZZ(mu*a2)
-
-	Z4=Integers(4)
-
-	A,B,C,D=bloc_decomposition(M1)
-
-	Ap=matrix(Z4,[[ZZ(-B[0,j]*a1-A[0,j]*a2-B[2,j]) for j in range(4)],
-		[ZZ(A[0,j]*a1-B[0,j]*a2+A[2,j]*q) for j in range(4)],
-		[ZZ(-B[1,j]*a1-A[1,j]*a2-B[3,j]) for j in range(4)],
-		[ZZ(A[1,j]*a1-B[1,j]*a2+A[3,j]*q) for j in range(4)]])
-
-	Bp=matrix(Z4,[[ZZ(2**m*(B[2,j]*a-A[0,j]*lamb-A[2,j]*bq)) for j in range(4)],
-		[ZZ(2**m*(B[2,j]*c+B[0,j]*mu-A[2,j]*dq)) for j in range(4)],
-		[ZZ(2**m*(B[3,j]*a-A[1,j]*lamb-A[3,j]*bq)) for j in range(4)],
-		[ZZ(2**m*(B[3,j]*c+B[1,j]*mu-A[3,j]*dq)) for j in range(4)]])
-
-	Cp=matrix(Z4,[[ZZ(ZZ(-D[0,j]*a1-C[0,j]*a2-D[2,j])//(2**m)) for j in range(4)],
-		[ZZ(ZZ(C[0,j]*a1-D[0,j]*a2+C[2,j]*q)//2**m) for j in range(4)],
-		[ZZ(ZZ(-D[1,j]*a1-C[1,j]*a2-D[3,j])//2**m) for j in range(4)],
-		[ZZ(ZZ(C[1,j]*a1-D[1,j]*a2+C[3,j]*q)//2**m) for j in range(4)]])
-
-	Dp=matrix(Z4,[[ZZ(D[2,j]*a-C[0,j]*lamb-C[2,j]*bq) for j in range(4)],
-		[ZZ(D[0,j]*mu+D[2,j]*c-C[2,j]*dq) for j in range(4)],
-		[ZZ(D[3,j]*a-C[1,j]*lamb-C[3,j]*bq) for j in range(4)],
-		[ZZ(D[1,j]*mu+D[3,j]*c-C[3,j]*dq) for j in range(4)]])
-
-	M=block_matrix(2,2,[[Ap,Cp],[Bp,Dp]])
-
-	assert is_symplectic_matrix_dim4(M)
-
-	return M
-
-def gluing_base_change_matrix_dim2_F2(a1,a2,q):
-	r"""Computes the symplectic base change matrix of a symplectic basis (*,B_K4) of E1*E2[4]
-	given by the kernel of the dimension 2 gluing isogeny:
-	B_K4=2**(f-2)[([a1]P1-[a2]P2,R1),([a1]P2+[a2]P1,R2)]
-	in the basis $2**(f-2)*B1$ given by:
-	B1:=[[(P1,0),(0,R1)],[(P2,0),(0,[1/q]*R2)]]
-	where:
-	- lamb is the inverse of q modulo 2**f.
-	- (P1,P2) is the canonical basis of E1[2**f].
-	- (R1,R2) is the image of (P1,P2) by sigma.
-
-	Input:
-	- a1, q: integers.
-
-	Output:
-	- M: symplectic base change matrix of (*,B_K4) in 2**(f-2)*B1.
-	"""
-
-	Z4=Integers(4)
-
-	mu=inverse_mod(a1,4)
-
-	A=matrix(Z4,[[0,mu],
-		[0,0]])
-	B=matrix(Z4,[[0,0],
-		[1,-ZZ(mu*a2)]])
-
-	C=matrix(Z4,[[ZZ(a1),-ZZ(a2)],
-		[-1,0]])
-	D=matrix(Z4,[[ZZ(a2),ZZ(a1)],
-		[0,-ZZ(q)]])
-
-	M=block_matrix([[A,C],[B,D]])
-
-	assert is_symplectic_matrix_dim2(M)
-
-	return M
-
-def gluing_base_change_matrix_dim2_dim4_F2(a1,a2,q,m,M2):
-	r"""Computes the symplectic base change matrix of the symplectic basis Bp of Am*Am[4] induced
-	by the image of the symplectic basis (x_1, ..., x_4, y_1, ..., y_4) of E1^2*E2^2[2**(e1+2)]
-	adapted to ker(F1)=[4]<y_1, ..., y_4> in the basis associated to the product theta-structure 
-	of level 2 of Am*Am:
-	
-	B:=[(S1,0),(S2,0),(0,S1),(0,S2),(T1,0),(T2,0),(0,T1),(0,T2)]
-	
-	where:
-	- (P1,Q1) is the canonical basis of E1[2**f].
-	- (R2,S2) is the image of (P1,P2) by sigma.
-	- Phi is the 2**m-isogeny E1*E2-->Am (m first steps of the chain in dimension 2).
-	- S1=[2**e]Phi([lamb]Q1,[a]sigma(P1)+[b]sigma(Q1)).
-	- S2=[2**e]Phi([mu]P1,[c]sigma(P1)+[d]sigma(Q1)).
-	- T1=[2**(e-m)]Phi([a1]P1-[a2]Q1,sigma(P1)).
-	- T2=[2**(e-m)]Phi([a1]Q1+[a2]P1,sigma(Q1)).
-	- (S1,S2,T1,T2) is induced by the image by Phi of a symplectic basis of E1*E2[2**(m+2)] lying
-	above the symplectic basis of E1*E2[4] outputted by gluing_base_change_matrix_dim2.
-
-	INPUT:
-	- a1, a2, q: integers.
-	- m: integer (number of steps in dimension 2 and 2-adic valuation of a2).
-	- M2: matrix of (x_1, ..., x_4, y_1, ..., y_4) in the symplectic basis of E1^2*E2^2[2**(e1+2)]
-	given by:
-
-	B1:=[[(P1,0,0,0),(0,P1,0,0),(0,0,R2,0),(0,0,0,R2)],
-	[(Q1,0,0,0),(0,Q1,0,0),(0,0,[1/q]*S2,0),(0,0,0,[1/q]*S2)]]
-
-	OUTPUT:
-	- M: symplectic base change matrix of Bp in B.
-	"""
-
-	inv_a1=inverse_mod(a1,2**(m+2))
-	inv_q=inverse_mod(q,2**(m+2))
-	lamb=ZZ(2**(m+1))
-	mu=ZZ((1+2**(m+1)*q)*inv_a1)
-	a=ZZ(2**(m+1)*a2*inv_q)
-	bq=ZZ((1+2**(m+1)*a1))
-	c=ZZ(2**(m+1))
-	dq=-ZZ(mu*a2)
-
-	Z4=Integers(4)
-
-	A,B,C,D=bloc_decomposition(M2)
-
-	Ap=matrix(Z4,[[ZZ(-B[0,j]*a1+A[0,j]*a2+B[2,j]) for j in range(4)],
-		[ZZ(A[0,j]*a1+B[0,j]*a2-A[2,j]*q) for j in range(4)],
-		[ZZ(-B[1,j]*a1+A[1,j]*a2+B[3,j]) for j in range(4)],
-		[ZZ(A[1,j]*a1+B[1,j]*a2-A[3,j]*q) for j in range(4)]])
-
-	Bp=matrix(Z4,[[ZZ(2**m*(B[2,j]*a-A[0,j]*lamb-A[2,j]*bq)) for j in range(4)],
-		[ZZ(2**m*(B[2,j]*c+B[0,j]*mu-A[2,j]*dq)) for j in range(4)],
-		[ZZ(2**m*(B[3,j]*a-A[1,j]*lamb-A[3,j]*bq)) for j in range(4)],
-		[ZZ(2**m*(B[3,j]*c+B[1,j]*mu-A[3,j]*dq)) for j in range(4)]])
-
-	Cp=matrix(Z4,[[ZZ(ZZ(-D[0,j]*a1+C[0,j]*a2+D[2,j])//(2**m)) for j in range(4)],
-		[ZZ(ZZ(C[0,j]*a1+D[0,j]*a2-C[2,j]*q)//2**m) for j in range(4)],
-		[ZZ(ZZ(-D[1,j]*a1+C[1,j]*a2+D[3,j])//2**m) for j in range(4)],
-		[ZZ(ZZ(C[1,j]*a1+D[1,j]*a2-C[3,j]*q)//2**m) for j in range(4)]])
-
-	Dp=matrix(Z4,[[ZZ(D[2,j]*a-C[0,j]*lamb-C[2,j]*bq) for j in range(4)],
-		[ZZ(D[0,j]*mu+D[2,j]*c-C[2,j]*dq) for j in range(4)],
-		[ZZ(D[3,j]*a-C[1,j]*lamb-C[3,j]*bq) for j in range(4)],
-		[ZZ(D[1,j]*mu+D[3,j]*c-C[3,j]*dq) for j in range(4)]])
-
-	M=block_matrix(2,2,[[Ap,Cp],[Bp,Dp]])
-
-	assert is_symplectic_matrix_dim4(M)
-
-	return M
-
-def point_matrix_product(M,L_P,J=None,modulus=None):
-	r"""
-	Input:
-	- M: matrix with (modular) integer values.
-	- L_P: list of elliptic curve points [P1,P2,R1,R2] such that the rows of M correspond to the vectors
-	(P1,0,0,0),(0,P1,0,0),(0,0,R1,0),(0,0,0,R1),(P2,0,0,0),(0,P2,0,0),(0,0,R2,0),(0,0,0,R2).
-	- J: list of column indices (default, all the columns).
-	- modulus: order of points in L_P (default, None).
-
-	Output:
-	- L_ret: list of points corresponding to the columns of M with indices in J.
-	"""
-	if modulus==None:
-		M1=M
-	else:
-		Zmod=Integers(modulus)
-		M1=matrix(Zmod,M)
-
-	if J==None:
-		J=range(M1.ncols())
-
-	L_ret=[]
-	for j in J:
-		L_ret.append(TuplePoint(M1[0,j]*L_P[0]+M1[4,j]*L_P[1],M1[1,j]*L_P[0]+M1[5,j]*L_P[1],
-			M1[2,j]*L_P[2]+M1[6,j]*L_P[3],M1[3,j]*L_P[2]+M1[7,j]*L_P[3]))
-
-	return L_ret
-
-
-def kernel_basis(M,ei,mP1,mP2,mR1,mlambR2):
-	r"""
-	Input:
-	- M: matrix of a symplectic basis in the basis
-	B1:=[[(P1,0,0,0),(0,P1,0,0),(0,0,R1,0),(0,0,0,R1)],
-	[(P2,0,0,0),(0,P2,0,0),(0,0,lamb*R2,0),(0,0,0,lamb*R2)]]
-	as defined in complete_kernel_matrix_F1.
-	- ei: length of F1 or F2.
-	- mP1,mP2: canonical basis (P1,P2) of E1[2**f] multiplied by m:=2**(f-ei-2).
-	- mR1,mlambR2: (mR1,mlambR2)=(m*sigma(P1),m*sigma(P2)), where lamb is the
-	inverse of q=deg(sigma) modulo 2**f.
-
-	Output:
-	- Basis of the second symplectic subgroup basis of E1^2*E2^2[2**(ei+2)] induced by M. 
-	"""
-	modulus=2**(ei+2)
-
-	return point_matrix_product(M,[mP1,mP2,mR1,mlambR2],[4,5,6,7],modulus)
-
-def base_change_canonical_dim4(P1,P2,R1,R2,q,f,e1,e2):
-	lamb=inverse_mod(q,2**f)
-	
-	lambR2=lamb*R2
-
-	P1_doubles=[P1]
-	P2_doubles=[P2]
-	R1_doubles=[R1]
-	lambR2_doubles=[lambR2]
-
-	for i in range(f-2):
-		P1_doubles.append(2*P1_doubles[-1])
-		P2_doubles.append(2*P2_doubles[-1])
-		R1_doubles.append(2*R1_doubles[-1])
-		lambR2_doubles.append(2*lambR2_doubles[-1])
-
-	# Constructing canonical basis of E1[4] and E2[4].
-	_,_,T1,T2,MT=make_canonical(P1_doubles[-1],P2_doubles[-1],4,preserve_pairing=True)
-	_,_,U1,U2,MU=make_canonical(R1_doubles[-1],lambR2_doubles[-1],4,preserve_pairing=True)
-
-	# Base change matrix of the symplectic basis 2**(f-2)*B1 of E1^2*E2^2[4] in the basis:
-	# B0:=[[(T1,0,0,0),(0,T1,0,0),(0,0,U1,0),(0,0,0,U1)],
-	#[(T2,0,0,0),(0,T2,0,0),(0,0,U2,0),(0,0,0,U2)]]
-	# where B1:=[[(P1,0,0,0),(0,P1,0,0),(0,0,R1,0),(0,0,0,R1)],
-	#[(P2,0,0,0),(0,P2,0,0),(0,0,lamb*R2,0),(0,0,0,lamb*R2)]]
-	Z4=Integers(4)
-	M0=matrix(Z4,[[MT[0,0],0,0,0,MT[1,0],0,0,0],
-		[0,MT[0,0],0,0,0,MT[1,0],0,0],
-		[0,0,MU[0,0],0,0,0,MU[1,0],0],
-		[0,0,0,MU[0,0],0,0,0,MU[1,0]],
-		[MT[0,1],0,0,0,MT[1,1],0,0,0],
-		[0,MT[0,1],0,0,0,MT[1,1],0,0],
-		[0,0,MU[0,1],0,0,0,MU[1,1],0],
-		[0,0,0,MU[0,1],0,0,0,MU[1,1]]])
-
-	return P1_doubles,P2_doubles,R1_doubles,lambR2_doubles,T1,T2,U1,U2,MT,MU,M0
diff --git a/theta_lib/isogenies/Kani_clapoti.py b/theta_lib/isogenies/Kani_clapoti.py
deleted file mode 100644
index 66030e2..0000000
--- a/theta_lib/isogenies/Kani_clapoti.py
+++ /dev/null
@@ -1,258 +0,0 @@
-from sage.all import *
-import itertools
-
-from ..basis_change.kani_base_change import clapoti_cob_splitting_matrix
-from ..basis_change.base_change_dim4 import base_change_theta_dim4
-from ..theta_structures.Tuple_point import TuplePoint
-from ..theta_structures.montgomery_theta import null_point_to_montgomery_coeff, theta_point_to_montgomery_point
-from ..theta_structures.theta_helpers_dim4 import product_to_theta_points_dim4
-from ..utilities.supersingular import torsion_basis_to_Fp_rational_point
-from .Kani_gluing_isogeny_chain_dim4 import KaniClapotiGluing
-from .isogeny_chain_dim4 import IsogenyChainDim4
-
-class KaniClapotiIsog(IsogenyChainDim4):
-	r"""Class representing the 4-dimensional isogeny obtained via Kani's lemma F: Eu^2*Ev^2 --> Ea^2*A 
-	where Ea=[\mf{a}]*E is the result of the ideal class group action by \mf{a} when given relevant 
-	constants and torsion point information.
-
-	INPUT: 
-	- Pu, Qu = phi_u(P, Q)\in Eu;
-	- Pv, Qv = phi_v*phi_{ck}*\hat{\phi}_{bk}(P, Q)\in Ev;
-	- gu, xu, yu, gv, xv, yv, Nbk, Nck, e: positive integers;
-	where:
-		* gu(xu^2+yu^2)Nbk+gv(xv^2+yv^2)Nck=2^e;
-		* gcd(u*Nbk,v*Nck)=1 with u:=gu(xu^2+yu^2) and v:=gv(xv^2+yv^2);
-		* xu and xv are odd and yu and yv are even;
-		* \mf{b}=\mf{be}*\mf{bk} is a product of ideals of norms Nbe and Nbk respectively,
-		where Nbe is a product of small Elkies primes;
-		* \mf{c}=\mf{ce}*\mf{ck} is a product of ideals of norms Nce and Nck respectively,
-		where Nbe is a product of small Elkies primes;
-		* phi_{bk}: E --> E1 and phi_{ck}: E --> E2 are induced by the action of 
-		ideals \mf{bk} and \mf{ck} respectively;
-		* <P,Q>=E_1[2^{e+2}];
-		* phi_u: E1 --> Eu and phi_v: E2 --> Ev are gu and gv-isogenies respectively.
-
-	OUTPUT: F: Eu^2*Ev^2 --> Ea^2*A is the isogeny:
-	
-	F := [[Phi_{bp}*\tilde{Phi}_u, Phi_{cp}*\tilde{Phi}_v],
-		 [-Psi, \tilde{Phi}]]
-	
-	obtained from the Kani isogeny diamond:
-
-	A --------------------Phi------------------> Ev^2
-	^                                             ^
-	|                                             |
-	|                                           Phi_v
-	|                                             |
-	Psi                                          E2^2
-	|                                             ^
-	|                                             |
-	|                                     \tilde{Phi}_{ck}
-	|                                             |
-	Eu^2 --\tilde{Phi}_{u}--> E1^2 --Phi_{bk}--> Ea^2 
-
-	where Phi_{bk}:=Diag(phi_{bk},phi_{bk}), Phi_{ck}:=Diag(phi_{ck},phi_{ck}),
-
-	Phi_u := [[xu, -yu],
-			[yu, xu]] * Diag(phi_u,phi_u)
-	
-	Phi_v := [[xv, -yv],
-			[yv, xv]] * Diag(phi_v,phi_v)
-	"""
-
-	def __init__(self,points,integers,strategy=None):
-		gu,xu,yu,gv,xv,yv,Nbk,Nck,e = integers
-		Pu,Qu,Pv,Qv = points
-		if gu*(xu**2+yu**2)*Nbk+gv*(xv**2+yv**2)*Nck!=2**e:
-			raise ValueError("Wrong parameters: gu(xu^2+yu^2)Nbk + gv(xv^2+yv^2)Nck != 2^e")
-		if gcd(ZZ(gu*(xu**2+yu**2)*Nbk),ZZ(gv*(xv**2+yv**2)*Nck))!=1:
-			raise ValueError("Non coprime parameters: gcd(gu(xu^2+yu^2)Nbk, gv(xv^2+yv^2)Nck) != 1")
-		if xu%2==0:
-			xu,yu=yu,xu
-		if xv%2==0:
-			xv,yv=yv,xv
-
-		self.Eu = Pu.curve()
-		self.Ev = Pv.curve()
-		Fp2 = parent(Pu[0])
-		Fp = Fp2.base_ring()
-
-		# Number of dimension 2 steps before dimension 4 gluing Am^2-->B 
-		m=valuation(xv*yu-xu*yv,2)
-		integers=[gu,xu,yu,gv,xv,yv,Nbk,Nck,e,m]
-
-		points_mp3=[(2**(e-m-1))*P for P in points]
-		points_mp2=[2*P for P in points_mp3]
-		points_4=[(2**m)*P for P in points_mp2]
-
-		self.Ru_Fp = torsion_basis_to_Fp_rational_point(self.Eu,points_4[0],points_4[1],4)
-
-		self.gluing_isogeny_chain = KaniClapotiGluing(points_mp3,points_mp2,points_4,integers, coerce=Fp)
-
-		xuNbk = xu*Nbk
-		yuNbk = yu*Nbk
-		two_ep2 = 2**(e+2)
-		inv_Nbk = inverse_mod(Nbk,two_ep2)
-		u = gu*(xu**2+yu**2)
-		inv_u = inverse_mod(u,4)
-		lambxu = ((1-2**e*inv_u*inv_Nbk)*xu)%two_ep2
-		lambyu = ((1-2**e*inv_u*inv_Nbk)*yu)%two_ep2
-		xv_Nbk = (xv*inv_Nbk)%two_ep2
-		yv_Nbk = (yv*inv_Nbk)%two_ep2
-
-
-		B_Kpp = [TuplePoint(xuNbk*Pu,yuNbk*Pu,xv*Pv,yv*Pv),
-		TuplePoint(-yuNbk*Pu,xuNbk*Pu,-yv*Pv,xv*Pv),
-		TuplePoint(lambxu*Qu,lambyu*Qu,xv_Nbk*Qv,yv_Nbk*Qv),
-		TuplePoint(-lambyu*Qu,lambxu*Qu,-yv_Nbk*Qv,xv_Nbk*Qv)]
-
-		IsogenyChainDim4.__init__(self, B_Kpp, self.gluing_isogeny_chain, e, m, splitting=True, strategy=strategy)
-
-		# Splitting
-		M_split = clapoti_cob_splitting_matrix(integers)
-
-		self.N_split = base_change_theta_dim4(M_split,self.gluing_isogeny_chain.e4)
-
-		self.codomain_product = self._isogenies[-1]._codomain.base_change_struct(self.N_split)
-
-		# Extracting the group action image Ea=[\mathfrak{a}]*E from the codomain Ea^2*E'^2
-		self.theta_null_Ea, self.theta_null_Ep, self.Ea, self.Ep = self.extract_montgomery_curve()
-
-
-
-	def extract_montgomery_curve(self):
-
-		# Computing the theta null point of Ea
-		null_point=self.codomain_product.zero()
-		Fp2=parent(null_point[0])
-		Fp = Fp2.base_ring()
-		for i3, i4 in itertools.product([0,1],repeat=2):
-			if null_point[4*i3+8*i4]!=0:
-				i30=i3
-				i40=i4
-				theta_Ea_0=Fp(null_point[4*i3+8*i4])
-				theta_Ea_1=Fp(null_point[1+4*i3+8*i4])
-				break
-		for i1, i2 in itertools.product([0,1],repeat=2):
-			if null_point[i1+2*i2]!=0:
-				i10=i1
-				i20=i2
-				theta_Ep_0=Fp(null_point[i1+2*i2])
-				theta_Ep_1=Fp(null_point[i1+2*i2+4])
-				break
-
-		# Sanity check: is the codomain of F a product of the form Ea^2*E'^2 ?
-		theta_Ea=[Fp(theta_Ea_0),Fp(theta_Ea_1)]
-		theta_Ep=[Fp(theta_Ep_0),Fp(theta_Ep_1)]
-
-		theta_Ea2Ep2=[0 for i in range(16)]
-		for i1,i2,i3,i4 in itertools.product([0,1],repeat=4):
-			theta_Ea2Ep2[i1+2*i2+4*i3+8*i4]=theta_Ea[i1]*theta_Ea[i2]*theta_Ep[i3]*theta_Ep[i4]
-		theta_Ea2Ep2=self.codomain_product(theta_Ea2Ep2)
-
-		assert theta_Ea2Ep2.is_zero()
-
-		A_Ep = null_point_to_montgomery_coeff(theta_Ep_0,theta_Ep_1)
-		Ep = EllipticCurve([0,A_Ep,0,1,0])
-
-		## ## Recovering Ea over Fp and not Fp2
-		## self.find_Fp_rational_theta_struct_Ea()
-
-		## theta_Ea = self.iso_Ea(theta_Ea)
-		## A_Ea = null_point_to_montgomery_coeff(theta_Ea[0],theta_Ea[1])
-
-		## # Sanity check : the curve Ea should be defined over Fp
-		## # assert A_Ea[1] == 0
-
-		## # Twisting Ea if necessary: if A_Ea+2 is not a square in Fp, then we take the twist (A_Ea --> -A_Ea)
-		## p=self.Eu.base_field().characteristic()
-		## self.twist = False
-		## if (A_Ea+2)**((p-1)//2)==-1:
-		## 	A_Ea = -A_Ea
-		## 	self.twist = True
-
-		## Ea = EllipticCurve([0,A_Ea,0,1,0])
-
-		A = null_point_to_montgomery_coeff(theta_Ea_0, theta_Ea_1)
-		Ab = null_point_to_montgomery_coeff(theta_Ea_0+theta_Ea_1, theta_Ea_0-theta_Ea_1)
-		Acan = min([A, -A, Ab, -Ab])
-		Acan = A
-		if (Acan == A or Acan == -A):
-			# 'Id' corresponds to the point on the twist
-			self.iso_type = 'Id'
-		else:
-			# 'Hadamard' corresponds to the point on the curve
-			self.iso_type = 'Hadamard'
-		if ((self.iso_type == 'Hadamard' and not (Acan+2).is_square()) or (self.iso_type == 'Id' and (Acan+2).is_square())):
-			Acan=-Acan
-		if (Acan == A or Acan == Ab):
-			self.twist = False
-		else:
-			self.twist = True
-		Ea = EllipticCurve([0,Acan,0,1,0])
-
-		# Find the dual null point
-		return theta_Ea, theta_Ep, Ea, Ep
-
-	def eval_rational_point_4_torsion(self):
-		T = TuplePoint(self.Ru_Fp,self.Eu(0),self.Ev(0),self.Ev(0))
-
-		FPu_4 = self.evaluate_isogeny(T)
-		FPu_4=self.codomain_product.base_change_coords(self.N_split,FPu_4)
-		FPu_4=product_to_theta_points_dim4(FPu_4)
-
-		return FPu_4[0]
-
-	def find_Fp_rational_theta_struct_Ea(self):
-		Pa = self.eval_rational_point_4_torsion()
-
-		HPa = (Pa[0]+Pa[1],Pa[0]-Pa[1])
-		i = self.Eu.base_field().gen()
-		self.i = i
-		iHPa = (Pa[0]+i*Pa[1],Pa[0]-i*Pa[1])
-
-		if Pa[0]==0 or Pa[1]==0:
-			self.iso_type='Id'
-		elif HPa[0]==0 or HPa[1]==0:
-			self.iso_type='Hadamard'
-		elif iHPa[0]==0 or iHPa[1]==0:
-			self.iso_type='iHadamard'
-		else:
-			raise ValueError("A rational theta point should be mapped to (0:1) or (1:0) after change of theta coordinates on Ea.")
-
-	def iso_Ea(self,P):
-		# Change of theta coordinates to obtain Fp-rational theta coordinates on Ea
-
-		if self.iso_type == 'Id':
-			return P
-		elif self.iso_type == 'Hadamard':
-			return (P[0]+P[1],P[0]-P[1])
-		else:
-			return (P[0]+self.i*P[1],P[0]-self.i*P[1])
-
-
-	def evaluate(self,P):
-		FP=self.evaluate_isogeny(P)
-		FP=self.codomain_product.base_change_coords(self.N_split,FP)
-		
-		FP=product_to_theta_points_dim4(FP)
-		FP=TuplePoint([theta_point_to_montgomery_point(self.Ea,self.theta_null_Ea,self.iso_Ea(FP[0]),self.twist),
-			theta_point_to_montgomery_point(self.Ea,self.theta_null_Ea,self.iso_Ea(FP[1]),self.twist),
-			theta_point_to_montgomery_point(self.Ep,self.theta_null_Ep,FP[2]),
-			theta_point_to_montgomery_point(self.Ep,self.theta_null_Ep,FP[3])])
-
-		return FP
-
-	def __call__(self,P):
-		return self.evaluate(P)
-
-
-
-
-
-
-
-
-
-
-
diff --git a/theta_lib/isogenies/Kani_gluing_isogeny_chain_dim4.py b/theta_lib/isogenies/Kani_gluing_isogeny_chain_dim4.py
deleted file mode 100644
index 282219c..0000000
--- a/theta_lib/isogenies/Kani_gluing_isogeny_chain_dim4.py
+++ /dev/null
@@ -1,567 +0,0 @@
-from sage.all import *
-from ..utilities.discrete_log import weil_pairing_pari
-from ..basis_change.canonical_basis_dim1 import make_canonical
-from ..basis_change.kani_base_change import (
-    fixed_deg_gluing_matrix_Phi1,
-    fixed_deg_gluing_matrix_Phi2,
-    fixed_deg_gluing_matrix_dim4,
-    clapoti_cob_matrix_dim2,
-    clapoti_cob_matrix_dim2_dim4,
-    gluing_base_change_matrix_dim2,
-    gluing_base_change_matrix_dim2_dim4,
-    gluing_base_change_matrix_dim2_F1,
-    gluing_base_change_matrix_dim2_F2,
-    kernel_basis,
-)
-from ..basis_change.base_change_dim2 import base_change_theta_dim2
-from ..basis_change.base_change_dim4 import base_change_theta_dim4
-from ..theta_structures.Theta_dim1 import ThetaStructureDim1
-from ..theta_structures.Theta_dim2 import ProductThetaStructureDim2
-from ..theta_structures.Tuple_point import TuplePoint
-from ..theta_structures.Theta_dim4 import ProductThetaStructureDim2To4, ThetaPointDim4
-from ..isogenies_dim2.isogeny_chain_dim2 import IsogenyChainDim2
-from .gluing_isogeny_dim4 import GluingIsogenyDim4
-
-class KaniFixedDegDim2Gluing:
-	def __init__(self,P_mp3,Q_mp3,a,b,c,d,u,f,m,strategy_dim2=None):
-		r"""
-		INPUT:
-		- P_mp3, Q_mp3: basis of E[2^(m+3)] such that pi(P_mp3)=P_mp3 and pi(Q_mp3)=-Q_mp3.
-		- a,b,c,d,u,f: integers such that a**2+c**2+p*(b**2+d**2)=u*(2**f-u), where p is
-		ths characteristic of the base field.
-		- m: integer such that m=min(v_2(a-b),v_2(a+b)).
-
-		OUTPUT: Gluing isogeny chain F_{m+1}\circ...\circ F_1 containing the first m+1 steps of
-		the isogeny F: E^4 --> A*A' representing a u-isogeny in dimension 2.
-		"""
-
-		P_mp2 = 2*P_mp3
-		Q_mp2 = 2*Q_mp3
-		P_4 = 2**m*P_mp2
-		Q_4 = 2**m*Q_mp2
-
-		E = P_mp3.curve()
-
-		# Canonical basis with S_4=(1,0)
-		_, _, R_4, S_4, M_dim1 = make_canonical(P_4,Q_4,4,preserve_pairing=True)
-
-		Z4 = Integers(4)
-		M0 = matrix(Z4,[[M_dim1[0,0],0,M_dim1[0,1],0],
-			[0,M_dim1[0,0],0,M_dim1[0,1]],
-			[M_dim1[1,0],0,M_dim1[1,1],0],
-			[0,M_dim1[1,0],0,M_dim1[1,1]]])
-
-		# Theta structures
-		Theta_E = ThetaStructureDim1(E,R_4,S_4)
-		Theta_EE = ProductThetaStructureDim2(Theta_E,Theta_E)
-
-		# Gluing change of basis in dimension 2
-		M1 = fixed_deg_gluing_matrix_Phi1(u,a,b,c,d)
-		M2 = fixed_deg_gluing_matrix_Phi2(u,a,b,c,d)
-
-		M10 = M0*M1
-		M20 = M0*M2
-
-		Fp2 = E.base_field()
-		e4 = Fp2(weil_pairing_pari(R_4,S_4,4))
-
-		N_Phi1 = base_change_theta_dim2(M10,e4)
-		N_Phi2 = base_change_theta_dim2(M20,e4)
-
-		# Gluing change of basis dimension 2 * dimension 2 --> dimension 4
-		M_dim4 = fixed_deg_gluing_matrix_dim4(u,a,b,c,d,m)
-
-		self.N_dim4 = base_change_theta_dim4(M_dim4,e4)
-
-		# Kernel of Phi1 : E^2 --> A_m1 and Phi2 : E^2 --> A_m2
-		two_mp2 = 2**(m+2)
-		two_mp3 = 2*two_mp2
-		mu = inverse_mod(u,two_mp3)
-		
-		B_K_Phi1 = [TuplePoint((u%two_mp2)*P_mp2,((c+d)%two_mp2)*P_mp2),
-		TuplePoint((((d**2-c**2)*mu)%two_mp2)*Q_mp2,((c-d)%two_mp2)*Q_mp2)]
-
-		B_K_Phi2 = [TuplePoint((u%two_mp2)*P_mp2,((d-c)%two_mp2)*P_mp2),
-		TuplePoint((((d**2-c**2)*mu)%two_mp2)*Q_mp2,(-(c+d)%two_mp2)*Q_mp2)]
-
-		# Computation of the 2**m-isogenies Phi1 and Phi2
-		self._Phi1=IsogenyChainDim2(B_K_Phi1,Theta_EE,N_Phi1,m,strategy_dim2)
-		self._Phi2=IsogenyChainDim2(B_K_Phi2,Theta_EE,N_Phi2,m,strategy_dim2)
-
-		# Kernel of the (m+1)-th isogeny in dimension 4 F_{m+1}: A_m1*A_m2 --> B (gluing isogeny)
-
-		B_K_dim4 =[TuplePoint((u%two_mp3)*P_mp3,E(0),((a+b)%two_mp3)*P_mp3,((c+d)%two_mp3)*P_mp3),
-		TuplePoint(E(0),(u%two_mp3)*P_mp3,((d-c)%two_mp3)*P_mp3,((a-b)%two_mp3)*P_mp3),
-		TuplePoint(((u-2**f)%two_mp3)*Q_mp3,E(0),((a-b)%two_mp3)*Q_mp3,((c-d)%two_mp3)*Q_mp3),
-		TuplePoint(E(0),((u-2**f)%two_mp3)*Q_mp3,((-c-d)%two_mp3)*Q_mp3,((a+b)%two_mp3)*Q_mp3)]
-
-		L_K_dim4=B_K_dim4+[B_K_dim4[0]+B_K_dim4[1]]
-
-		L_K_dim4=[[self._Phi1(TuplePoint(T[0],T[3])),self._Phi2(TuplePoint(T[1],T[2]))] for T in L_K_dim4]
-
-		# Product Theta structure on A_m1*A_m2
-		self.domain_product=ProductThetaStructureDim2To4(self._Phi1._codomain,self._Phi2._codomain)
-		
-		# Theta structure on A_m1*A_m2 after change of theta coordinates
-		self.domain_base_change=self.domain_product.base_change_struct(self.N_dim4)
-		
-		# Converting the kernel to the Theta structure domain_base_change
-		L_K_dim4=[self.domain_product.product_theta_point(T[0],T[1]) for T in L_K_dim4]
-		L_K_dim4=[self.domain_base_change.base_change_coords(self.N_dim4,T) for T in L_K_dim4]
-
-		# Computing the gluing isogeny in dimension 4
-		self._gluing_isogeny_dim4=GluingIsogenyDim4(self.domain_base_change,L_K_dim4,[(1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1),(1,1,0,0)])
-
-		# Translates for the evaluation of the gluing isogeny in dimension 4
-		self.L_trans=[2*B_K_dim4[k] for k in range(2)]
-		self.L_trans_ind=[1,2] # Corresponds to multi indices (1,0,0,0) and (0,1,0,0)
-
-		self._codomain=self._gluing_isogeny_dim4._codomain
-
-	def evaluate(self,P):
-		if not isinstance(P, TuplePoint):
-			raise TypeError("KaniGluingIsogenyChainDim4 isogeny expects as input a TuplePoint on the domain product E^4")
-
-		# Translating P
-		L_P_trans=[P+T for T in self.L_trans]
-		
-		# dim4 --> dim2 x dim2
-		eval_P=[TuplePoint(P[0],P[3]),TuplePoint(P[1],P[2])]
-		eval_L_P_trans=[[TuplePoint(Q[0],Q[3]),TuplePoint(Q[1],Q[2])] for Q in L_P_trans]
-
-		# evaluating through the dimension 2 isogenies
-		eval_P=[self._Phi1(eval_P[0]),self._Phi2(eval_P[1])]
-		eval_L_P_trans=[[self._Phi1(Q[0]),self._Phi2(Q[1])] for Q in eval_L_P_trans]
-
-		# Product Theta structure and base change
-		eval_P=self.domain_product.product_theta_point(eval_P[0],eval_P[1])
-		eval_P=self.domain_base_change.base_change_coords(self.N_dim4,eval_P)
-
-		eval_L_P_trans=[self.domain_product.product_theta_point(Q[0],Q[1]) for Q in eval_L_P_trans]
-		eval_L_P_trans=[self.domain_base_change.base_change_coords(self.N_dim4,Q) for Q in eval_L_P_trans]
-
-		return self._gluing_isogeny_dim4.special_image(eval_P,eval_L_P_trans,self.L_trans_ind)
-
-	def __call__(self,P):
-		return self.evaluate(P)
-
-
-class KaniClapotiGluing:
-	def __init__(self,points_mp3,points_mp2,points_4,integers,strategy_dim2=None,coerce=None):
-		self._coerce=coerce
-		Pu_mp3,Qu_mp3,Pv_mp3,Qv_mp3 = points_mp3
-		Pu_mp2,Qu_mp2,Pv_mp2,Qv_mp2 = points_mp2
-		Pu_4,Qu_4,Pv_4,Qv_4 = points_4
-		gu,xu,yu,gv,xv,yv,Nbk,Nck,e,m = integers
-
-		Eu=Pu_4.curve()
-		Ev=Pv_4.curve()
-
-		lamb_u = inverse_mod(ZZ(gu),4)
-		lamb_v = inverse_mod(ZZ(gv*Nbk*Nck),4)
-
-
-		# 4-torsion canonical change of basis in Eu and Ev (Su=(1,0), Sv=(1,0))
-		_,_,Ru,Su,Mu=make_canonical(Pu_4,lamb_u*Qu_4,4,preserve_pairing=True)
-		_,_,Rv,Sv,Mv=make_canonical(Pv_4,lamb_v*Qv_4,4,preserve_pairing=True)
-
-		Z4 = Integers(4)
-		M0=matrix(Z4,[[Mu[0,0],0,Mu[1,0],0],
-			[0,Mv[0,0],0,Mv[1,0]],
-			[Mu[0,1],0,Mu[1,1],0],
-			[0,Mv[0,1],0,Mv[1,1]]])
-
-		self.M_product_dim2=M0
-
-		# Theta structures in dimension 1 and 2
-		Theta_u=ThetaStructureDim1(Eu,Ru,Su)
-		Theta_v=ThetaStructureDim1(Ev,Rv,Sv)
-
-		Theta_uv=ProductThetaStructureDim2(Theta_u,Theta_v)
-
-		# Gluing change of basis in dimension 2
-		M1 = clapoti_cob_matrix_dim2(integers)
-		M10 = M0*M1
-
-		Fp2 = Eu.base_field()
-		e4 = Fp2(weil_pairing_pari(Ru,Su,4))
-		self.e4 = e4
-
-		N_dim2 = base_change_theta_dim2(M10,e4)
-
-		# Gluing change of basis dimension 2 * dimension 2 --> dimension 4
-		M2 = clapoti_cob_matrix_dim2_dim4(integers)
-
-		self.N_dim4 = base_change_theta_dim4(M2,e4)
-
-		# Kernel of the 2**m-isogeny chain in dimension 2
-		two_mp2=2**(m+2)
-		two_mp3=2*two_mp2
-		u=ZZ(gu*(xu**2+yu**2))
-		mu=inverse_mod(u,two_mp2)
-		suv=ZZ(xu*xv+yu*yv)
-		duv=ZZ(xv*yu-xu*yv)
-		uNbk=(u*Nbk)%two_mp2
-		gusuv=(gu*suv)%two_mp2
-		xK2=(uNbk+gu*gv*mu*Nck*duv**2)%two_mp2
-		B_K_dim2 = [TuplePoint(uNbk*Pu_mp2,gusuv*Pv_mp2),TuplePoint(xK2*Qu_mp2,gusuv*Qv_mp2)]
-
-		# Computation of the 2**m-isogeny chain in dimension 2
-		self._isogenies_dim2=IsogenyChainDim2(B_K_dim2,Theta_uv,N_dim2,m,strategy_dim2)
-
-		# Kernel of the (m+1)-th isogeny in dimension 4 f_{m+1}: A_m^2 --> B (gluing isogeny)
-		xuNbk = (xu*Nbk)%two_mp3
-		yuNbk = (yu*Nbk)%two_mp3
-		inv_Nbk = inverse_mod(Nbk,two_mp3)
-		lambxu = ((1-2**e)*xu)%two_mp3 # extreme case m=e-2, 2^e = 2^(m+2) so 2^e/(u*Nbk) = 2^e mod 2^(m+3).
-		lambyu = ((1-2**e)*yu)%two_mp3
-		xv_Nbk = (xv*inv_Nbk)%two_mp3
-		yv_Nbk = (yv*inv_Nbk)%two_mp3
-
-		B_K_dim4 = [TuplePoint(xuNbk*Pu_mp3,yuNbk*Pu_mp3,xv*Pv_mp3,yv*Pv_mp3),
-		TuplePoint(-yuNbk*Pu_mp3,xuNbk*Pu_mp3,-yv*Pv_mp3,xv*Pv_mp3),
-		TuplePoint(lambxu*Qu_mp3,lambyu*Qu_mp3,xv_Nbk*Qv_mp3,yv_Nbk*Qv_mp3),
-		TuplePoint(-lambyu*Qu_mp3,lambxu*Qu_mp3,-yv_Nbk*Qv_mp3,xv_Nbk*Qv_mp3)]
-
-		L_K_dim4=B_K_dim4+[B_K_dim4[0]+B_K_dim4[1]]
-
-		L_K_dim4=[[self._isogenies_dim2(TuplePoint(L_K_dim4[k][0],L_K_dim4[k][2])),self._isogenies_dim2(TuplePoint(L_K_dim4[k][1],L_K_dim4[k][3]))] for k in range(5)]
-
-		# Product Theta structure on A_m^2
-		self.domain_product=ProductThetaStructureDim2To4(self._isogenies_dim2._codomain,self._isogenies_dim2._codomain)
-		
-		# Theta structure on A_m^2 after change of theta coordinates
-		self.domain_base_change=self.domain_product.base_change_struct(self.N_dim4)
-		
-		# Converting the kernel to the Theta structure domain_base_change
-		L_K_dim4=[self.domain_product.product_theta_point(L_K_dim4[k][0],L_K_dim4[k][1]) for k in range(5)]
-		L_K_dim4=[self.domain_base_change.base_change_coords(self.N_dim4,L_K_dim4[k]) for k in range(5)]
-
-		# Computing the gluing isogeny in dimension 4
-		self._gluing_isogeny_dim4=GluingIsogenyDim4(self.domain_base_change,L_K_dim4,[(1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1),(1,1,0,0)], coerce=self._coerce)
-
-		# Translates for the evaluation of the gluing isogeny in dimension 4
-		self.L_trans=[2*B_K_dim4[k] for k in range(2)]
-		self.L_trans_ind=[1,2] # Corresponds to multi indices (1,0,0,0) and (0,1,0,0)
-
-		self._codomain=self._gluing_isogeny_dim4._codomain
-
-	def evaluate(self,P):
-		if not isinstance(P, TuplePoint):
-			raise TypeError("KaniGluingIsogenyChainDim4 isogeny expects as input a TuplePoint on the domain product Eu^2 x Ev^2")
-
-		# Translating P
-		L_P_trans=[P+T for T in self.L_trans]
-		
-		# dim4 --> dim2 x dim2
-		eval_P=[TuplePoint(P[0],P[2]),TuplePoint(P[1],P[3])]
-		eval_L_P_trans=[[TuplePoint(Q[0],Q[2]),TuplePoint(Q[1],Q[3])] for Q in L_P_trans]
-
-		# evaluating through the dimension 2 isogenies
-		eval_P=[self._isogenies_dim2(eval_P[0]),self._isogenies_dim2(eval_P[1])]
-		eval_L_P_trans=[[self._isogenies_dim2(Q[0]),self._isogenies_dim2(Q[1])] for Q in eval_L_P_trans]
-
-		# Product Theta structure and base change
-		eval_P=self.domain_product.product_theta_point(eval_P[0],eval_P[1])
-		eval_P=self.domain_base_change.base_change_coords(self.N_dim4,eval_P)
-
-		eval_L_P_trans=[self.domain_product.product_theta_point(Q[0],Q[1]) for Q in eval_L_P_trans]
-		eval_L_P_trans=[self.domain_base_change.base_change_coords(self.N_dim4,Q) for Q in eval_L_P_trans]
-
-		return self._gluing_isogeny_dim4.special_image(eval_P,eval_L_P_trans,self.L_trans_ind)
-
-	def __call__(self,P):
-		return self.evaluate(P)
-
-		
-		
-class KaniGluingIsogenyChainDim4:
-	def __init__(self,points_m,points_4,a1,a2,q,m,strategy_dim2=None):
-		r"""
-
-		INPUT: 
-		- points_m: list of 4 points P1_m, Q1_m, R2_m, S2_m of order 2**(m+3)
-		such that (P1_m,Q1_m) generates E1[2**(m+3)] and (R2_m,S2_m) is
-		its image by sigma: E1 --> E2.
-		- points_4: list of 4 points P1_4, Q1_4, R2_4, S2_4 of order 4 obtained by
-		multiplying P1_m, Q1_m, R2_m, S2_m by 2**(m+1).
-		- a1, a2, q: integers such that a1**2+a2**2+q=2**e.
-		- m: 2-adic valuation of a2.
-
-		OUTPUT: Composition of the m+1 first isogenies in the isogeny chained
-		E1^2*E1^2 --> E1^2*E2^2 parametrized by a1, a2 and sigma via Kani's lemma.
-		"""
-
-		P1_m, Q1_m, R2_m, S2_m = points_m
-		P1_4, Q1_4, R2_4, S2_4 = points_4
-
-		E1=P1_m.curve()
-		E2=R2_m.curve()
-
-		Fp2=E1.base_field()
-
-		lamb=inverse_mod(q,4)
-
-		_,_,T1,T2,MT=make_canonical(P1_4,Q1_4,4,preserve_pairing=True)
-		_,_,U1,U2,MU=make_canonical(R2_4,lamb*S2_4,4,preserve_pairing=True)
-
-		Z4=Integers(4)
-		M0=matrix(Z4,[[MT[0,0],0,MT[1,0],0],
-			[0,MU[0,0],0,MU[1,0]],
-			[MT[0,1],0,MT[1,1],0],
-			[0,MU[0,1],0,MU[1,1]]])
-
-		self.M_product_dim2=M0
-
-		# Theta structures in dimension 1 and 2
-		Theta1=ThetaStructureDim1(E1,T1,T2)
-		Theta2=ThetaStructureDim1(E2,U1,U2)
-
-		Theta12=ProductThetaStructureDim2(Theta1,Theta2)
-
-		self.Theta1=Theta1
-		self.Theta2=Theta2
-		self.Theta12=Theta12
-
-		# Gluing base change in dimension 2
-		M1=gluing_base_change_matrix_dim2(a1,a2,q)
-		M10=M0*M1
-
-		self.M_gluing_dim2=M1
-
-		e4=Fp2(weil_pairing_pari(T1,T2,4))
-
-		self.e4=e4
-
-		N_dim2=base_change_theta_dim2(M10,e4)
-		#N_dim2=montgomery_to_theta_matrix_dim2(Theta12.zero().coords(),N1)
-
-		# Gluing base change in dimension 4
-		mua2=-M1[3,1]
-		M2=gluing_base_change_matrix_dim2_dim4(a1,a2,m,mua2)
-
-		self.M_gluing_dim4=M2
-
-		self.N_dim4=base_change_theta_dim4(M2,e4)
-
-		# Kernel of the 2**m-isogeny chain in dimension 2
-		a1_red=a1%(2**(m+2))
-		a2_red=a2%(2**(m+2))
-		B_K_dim2=[TuplePoint(2*a1_red*P1_m-2*a2_red*Q1_m,2*R2_m),TuplePoint(2*a1_red*Q1_m+2*a2_red*P1_m,2*S2_m)]
-
-		# Computation of the 2**m-isogeny chain in dimension 2
-		self._isogenies_dim2=IsogenyChainDim2(B_K_dim2,Theta12,N_dim2,m,strategy_dim2)
-
-		# Kernel of the (m+1)-th isogeny in dimension 4 f_{m+1}: A_m^2 --> B (gluing isogeny)
-		a1_red=a1%(2**(m+3))
-		a2_red=a2%(2**(m+3))
-		
-		a1P1_m=(a1_red)*P1_m
-		a2P1_m=(a2_red)*P1_m
-		a1Q1_m=(a1_red)*Q1_m
-		a2Q1_m=(a2_red)*Q1_m
-
-		OE2=E2(0)
-
-		B_K_dim4=[TuplePoint(a1P1_m,a2P1_m,R2_m,OE2),TuplePoint(a1Q1_m,a2Q1_m,S2_m,OE2),
-				TuplePoint(-a2P1_m,a1P1_m,OE2,R2_m),TuplePoint(-a2Q1_m,a1Q1_m,OE2,S2_m)]
-		L_K_dim4=B_K_dim4+[B_K_dim4[0]+B_K_dim4[1]]
-
-		L_K_dim4=[[self._isogenies_dim2(TuplePoint(L_K_dim4[k][0],L_K_dim4[k][2])),self._isogenies_dim2(TuplePoint(L_K_dim4[k][1],L_K_dim4[k][3]))] for k in range(5)]
-
-		# Product Theta structure on A_m^2
-		self.domain_product=ProductThetaStructureDim2To4(self._isogenies_dim2._codomain,self._isogenies_dim2._codomain)
-		
-		# Theta structure on A_m^2 after base change
-		self.domain_base_change=self.domain_product.base_change_struct(self.N_dim4)
-		
-		# Converting the kernel to the Theta structure domain_base_change
-		L_K_dim4=[self.domain_product.product_theta_point(L_K_dim4[k][0],L_K_dim4[k][1]) for k in range(5)]
-		L_K_dim4=[self.domain_base_change.base_change_coords(self.N_dim4,L_K_dim4[k]) for k in range(5)]
-
-		# Computing the gluing isogeny in dimension 4
-		self._gluing_isogeny_dim4=GluingIsogenyDim4(self.domain_base_change,L_K_dim4,[(1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1),(1,1,0,0)])
-
-		# Translates for the evaluation of the gluing isogeny in dimension 4
-		self.L_trans=[2*B_K_dim4[k] for k in range(2)]
-		self.L_trans_ind=[1,2] # Corresponds to multi indices (1,0,0,0) and (0,1,0,0)
-
-		self._codomain=self._gluing_isogeny_dim4._codomain
-
-	def evaluate(self,P):
-		if not isinstance(P, TuplePoint):
-			raise TypeError("KaniGluingIsogenyChainDim4 isogeny expects as input a TuplePoint on the domain product E1^2 x E2^2")
-
-		# Translating P
-		L_P_trans=[P+T for T in self.L_trans]
-		
-		# dim4 --> dim2 x dim2
-		eval_P=[TuplePoint(P[0],P[2]),TuplePoint(P[1],P[3])]
-		eval_L_P_trans=[[TuplePoint(Q[0],Q[2]),TuplePoint(Q[1],Q[3])] for Q in L_P_trans]
-
-		# evaluating through the dimension 2 isogenies
-		eval_P=[self._isogenies_dim2(eval_P[0]),self._isogenies_dim2(eval_P[1])]
-		eval_L_P_trans=[[self._isogenies_dim2(Q[0]),self._isogenies_dim2(Q[1])] for Q in eval_L_P_trans]
-
-		# Product Theta structure and base change
-		eval_P=self.domain_product.product_theta_point(eval_P[0],eval_P[1])
-		eval_P=self.domain_base_change.base_change_coords(self.N_dim4,eval_P)
-
-		eval_L_P_trans=[self.domain_product.product_theta_point(Q[0],Q[1]) for Q in eval_L_P_trans]
-		eval_L_P_trans=[self.domain_base_change.base_change_coords(self.N_dim4,Q) for Q in eval_L_P_trans]
-
-		return self._gluing_isogeny_dim4.special_image(eval_P,eval_L_P_trans,self.L_trans_ind)
-
-	def __call__(self,P):
-		return self.evaluate(P)
-
-class KaniGluingIsogenyChainDim4Half:
-	def __init__(self, points_m, a1, a2, q, m, Theta12, M_product_dim2, M_start_dim4, M_gluing_dim4, e4, dual=False,strategy_dim2=None):#points_m,points_4,a1,a2,q,m,precomputed_data=None,dual=False,strategy_dim2=None):
-		r"""
-
-		INPUT: 
-		- points_m: list of 4 points P1_m, Q1_m, R2_m, S2_m of order 2**(m+3)
-		such that (P1_m,Q1_m) generates E1[2**(m+3)] and (R2_m,S2_m) is
-		its image by sigma: E1 --> E2.
-		- points_4: list of 4 points P1_4, Q1_4, R2_4, S2_4 of order 4 obtained by
-		multiplying P1_m, Q1_m, R2_m, S2_m by 2**(m+1).
-		- a1, a2, q: integers such that a1**2+a2**2+q=2**e.
-		- m: 2-adic valuation of a2.
-
-		OUTPUT: Composition of the m+1 first isogenies in the isogeny chained
-		E1^2*E1^2 --> E1^2*E2^2 parametrized by a1, a2 and sigma via Kani's lemma.
-		"""
-
-		P1_m, Q1_m, R2_m, S2_m = points_m
-
-		E1=P1_m.curve()
-		E2=R2_m.curve()
-
-		Fp2=E1.base_field()
-
-		self.M_product_dim2 = M_product_dim2
-
-		self.Theta12=Theta12
-
-		self.e4=e4
-
-		# Gluing base change in dimension 2
-		if not dual:
-			M1=gluing_base_change_matrix_dim2_F1(a1,a2,q)
-		else:
-			M1=gluing_base_change_matrix_dim2_F2(a1,a2,q)
-		
-		M10=M_product_dim2*M1
-
-		self.M_gluing_dim2=M1
-
-		self.e4=e4
-
-		N_dim2=base_change_theta_dim2(M10,e4)
-		#N_dim2=montgomery_to_theta_matrix_dim2(Theta12.zero().coords(),N1)
-
-		# Gluing base change in dimension 4
-
-		self.M_gluing_dim4 = M_gluing_dim4
-
-		self.N_dim4 = base_change_theta_dim4(M_gluing_dim4, e4)
-
-		# Kernel of the 2**m-isogeny chain in dimension 2
-		a1_red=a1%(2**(m+2))
-		a2_red=a2%(2**(m+2))
-		if not dual:
-			B_K_dim2=[TuplePoint(2*a1_red*P1_m-2*a2_red*Q1_m,2*R2_m),TuplePoint(2*a1_red*Q1_m+2*a2_red*P1_m,2*S2_m)]
-		else:
-			B_K_dim2=[TuplePoint(2*a1_red*P1_m+2*a2_red*Q1_m,-2*R2_m),TuplePoint(2*a1_red*Q1_m-2*a2_red*P1_m,-2*S2_m)]
-
-		# Computation of the 2**m-isogeny chain in dimension 2
-		self._isogenies_dim2=IsogenyChainDim2(B_K_dim2,Theta12,N_dim2,m,strategy_dim2)
-
-		# Kernel of the (m+1)-th isogeny in dimension 4 f_{m+1}: A_m^2 --> B (gluing isogeny)
-		lamb=inverse_mod(q,2**(m+3))
-		B_K_dim4=kernel_basis(M_start_dim4,m+1,P1_m,Q1_m,R2_m,lamb*S2_m)
-		L_K_dim4=B_K_dim4+[B_K_dim4[0]+B_K_dim4[1]]
-
-		L_K_dim4=[[self._isogenies_dim2(TuplePoint(L_K_dim4[k][0],L_K_dim4[k][2])),self._isogenies_dim2(TuplePoint(L_K_dim4[k][1],L_K_dim4[k][3]))] for k in range(5)]
-
-		# Product Theta structure on A_m^2
-		self.domain_product=ProductThetaStructureDim2To4(self._isogenies_dim2._codomain,self._isogenies_dim2._codomain)
-		
-		# Theta structure on A_m^2 after base change
-		self.domain_base_change=self.domain_product.base_change_struct(self.N_dim4)
-		
-		# Converting the kernel to the Theta structure domain_base_change
-		L_K_dim4=[self.domain_product.product_theta_point(L_K_dim4[k][0],L_K_dim4[k][1]) for k in range(5)]
-		L_K_dim4=[self.domain_base_change.base_change_coords(self.N_dim4,L_K_dim4[k]) for k in range(5)]
-
-		# Computing the gluing isogeny in dimension 4
-		self._gluing_isogeny_dim4=GluingIsogenyDim4(self.domain_base_change,L_K_dim4,[(1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1),(1,1,0,0)])
-
-		# Translates for the evaluation of the gluing isogeny in dimension 4
-		self.L_trans=[2*B_K_dim4[k] for k in range(2)]
-		self.L_trans_ind=[1,2] # Corresponds to multi indices (1,0,0,0) and (0,1,0,0)
-
-		self._codomain=self._gluing_isogeny_dim4._codomain
-
-	def evaluate(self,P):
-		if not isinstance(P, TuplePoint):
-			raise TypeError("KaniGluingIsogenyChainDim4 isogeny expects as input a TuplePoint on the domain product E1^2 x E2^2")
-
-		# Translating P
-		L_P_trans=[P+T for T in self.L_trans]
-		
-		# dim4 --> dim2 x dim2
-		eval_P=[TuplePoint(P[0],P[2]),TuplePoint(P[1],P[3])]
-		eval_L_P_trans=[[TuplePoint(Q[0],Q[2]),TuplePoint(Q[1],Q[3])] for Q in L_P_trans]
-
-		# evaluating through the dimension 2 isogenies
-		eval_P=[self._isogenies_dim2(eval_P[0]),self._isogenies_dim2(eval_P[1])]
-		eval_L_P_trans=[[self._isogenies_dim2(Q[0]),self._isogenies_dim2(Q[1])] for Q in eval_L_P_trans]
-
-		# Product Theta structure and base change
-		eval_P=self.domain_product.product_theta_point(eval_P[0],eval_P[1])
-		eval_P=self.domain_base_change.base_change_coords(self.N_dim4,eval_P)
-
-		eval_L_P_trans=[self.domain_product.product_theta_point(Q[0],Q[1]) for Q in eval_L_P_trans]
-		eval_L_P_trans=[self.domain_base_change.base_change_coords(self.N_dim4,Q) for Q in eval_L_P_trans]
-
-		return self._gluing_isogeny_dim4.special_image(eval_P,eval_L_P_trans,self.L_trans_ind)
-
-	def __call__(self,P):
-		return self.evaluate(P)
-
-	def dual(self):
-		domain = self._codomain.hadamard()
-		codomain_base_change = self.domain_base_change
-		codomain_product = self.domain_product
-		N_dim4 = self.N_dim4.inverse()
-		isogenies_dim2 = self._isogenies_dim2.dual()
-		splitting_isogeny_dim4 = self._gluing_isogeny_dim4.dual()
-
-		return KaniSplittingIsogenyChainDim4(domain, codomain_base_change, codomain_product, N_dim4, isogenies_dim2, splitting_isogeny_dim4)
-
-class KaniSplittingIsogenyChainDim4:
-	def __init__(self, domain, codomain_base_change, codomain_product, N_dim4, isogenies_dim2, splitting_isogeny_dim4):
-		self._domain = domain
-		self.codomain_base_change = codomain_base_change
-		self.codomain_product = codomain_product
-		self.N_dim4 = N_dim4
-		self._isogenies_dim2 = isogenies_dim2
-		self._splitting_isogeny_dim4 = splitting_isogeny_dim4
-
-	def evaluate(self,P):
-		if not isinstance(P, ThetaPointDim4):
-			raise TypeError("KaniSplittingIsogenyChainDim4 isogeny expects as input a ThetaPointDim4")
-
-		Q = self._splitting_isogeny_dim4(P)
-		Q = self.codomain_product.base_change_coords(self.N_dim4, Q)
-		Q1, Q2 = self.codomain_product.to_theta_points(Q)
-		Q1, Q2 = self._isogenies_dim2._domain(Q1.hadamard()), self._isogenies_dim2._domain(Q2.hadamard()) 
-
-		Q1 = self._isogenies_dim2(Q1)
-		Q2 = self._isogenies_dim2(Q2)
-
-		return TuplePoint(Q1[0],Q2[0],Q1[1],Q2[1])
-
-	def __call__(self,P):
-		return self.evaluate(P)
diff --git a/theta_lib/isogenies/gluing_isogeny_dim4.py b/theta_lib/isogenies/gluing_isogeny_dim4.py
deleted file mode 100644
index cd738c9..0000000
--- a/theta_lib/isogenies/gluing_isogeny_dim4.py
+++ /dev/null
@@ -1,200 +0,0 @@
-from sage.all import *
-
-from ..theta_structures.Theta_dim4 import ThetaStructureDim4
-from ..theta_structures.theta_helpers_dim4 import hadamard, squared, batch_inversion, multindex_to_index
-from .tree import Tree
-from .isogeny_dim4 import IsogenyDim4, DualIsogenyDim4
-
-def proj_equal(P1, P2):
-    if len(P1) != len(P2):
-        return False
-    for i in range(0, len(P1)):
-        if P1[i]==0:
-            if P2[i] != 0:
-                return False
-        else:
-            break
-    r=P1[i]
-    s=P2[i]
-    for i in range(0, len(P1)):
-        if P1[i]*s != P2[i]*r:
-            return False
-    return True
-
-class GluingIsogenyDim4(IsogenyDim4):
-	def __init__(self,domain,L_K_8,L_K_8_ind, coerce=None):
-		r"""
-		Input:
-		- domain: a ThetaStructureDim4.
-		- L_K_8: list of points of 8-torsion in the kernel.
-		- L_K_8_ind: list of corresponding multindices (i0,i1,i2,i3) of points in L_K_8
-		(L_K_8[i]=i0*P0+i1*P1+i2*P2+i3*P3, where (P0,..,P3) is a basis of K_8 (compatible with
-		the canonical basis of K_2) and L_K_8_ind[i]=(i0,i1,i2,i3)).
-		"""
-
-		if not isinstance(domain, ThetaStructureDim4):
-			raise ValueError("Argument domain should be a ThetaStructureDim4 object.")
-		self._domain = domain
-		self._precomputation=None
-		self._coerce=coerce
-		self._special_compute_codomain(L_K_8,L_K_8_ind)
-
-		#a_i2=squared(self._domain.zero())
-		#HB_i2=hadamard(squared(hadamard(self._codomain.zero())))
-		#for i in range(16):
-			#print(HB_i2[i]/a_i2[i])
-
-	def _special_compute_codomain(self,L_K_8,L_K_8_ind):
-		r"""
-		Input:
-		- L_K_8: list of points of 8-torsion in the kernel.
-		- L_K_8_ind: list of corresponding multindices (i0,i1,i2,i3) of points in L_K_8
-		(L_K_8[i]=i0*P0+i1*P1+i2*P2+i3*P3, where (P0,..,P3) is a basis of K_8 (compatible with
-		the canonical basis of K_2) and L_K_8_ind[i]=(i0,i1,i2,i3)).
-		
-		Output:
-		- codomain of the isogeny.
-		Also initializes self._precomputation, containing the inverse of theta-constants.
-		"""
-		HSK_8=[hadamard(squared(P.coords())) for P in L_K_8]
-		
-		# Choice of reference index j_0<->chi_0 corresponding to a non-vanishing theta-constant.
-		found_tree=False
-		j_0=0
-		while not found_tree:
-			found_k0=False
-			for k in range(len(L_K_8)):
-				if HSK_8[k][j_0]!=0:
-					k_0=k
-					found_k0=True
-					break
-			if not found_k0:
-				j_0+=1
-			else:
-				j0pk0=j_0^multindex_to_index(L_K_8_ind[k_0])
-				# List of tuples of indices (index chi of the denominator: HS(f(P_k))_chi, 
-				#index chi.chi_k of the numerator: HS(f(P_k))_chi.chi_k, index k).
-				L_ratios_ind=[(j_0,j0pk0,k_0)]
-				L_covered_ind=[j_0,j0pk0]
-
-				# Tree containing the the theta-null points indices as nodes and the L_ratios_ind reference indices as edges.
-				tree_ratios=Tree(j_0)
-				tree_ratios.add_child(Tree(j0pk0),0)
-
-				# Filling in the tree
-				tree_filled=False
-				while not tree_filled:
-					found_j=False
-					for j in L_covered_ind:
-						for k in range(len(L_K_8)):
-							jpk=j^multindex_to_index(L_K_8_ind[k])
-							if jpk not in L_covered_ind and HSK_8[k][j]!=0:
-								L_covered_ind.append(jpk)
-								L_ratios_ind.append((j,jpk,k))
-								tree_j=tree_ratios.look_node(j)
-								tree_j.add_child(Tree(jpk),len(L_ratios_ind)-1)
-								found_j=True
-								#break
-							#if found_j:
-								#break
-					if not found_j or len(L_covered_ind)==16:
-						tree_filled=True
-				if len(L_covered_ind)!=16:
-					j_0+=1
-				else:
-					found_tree=True
-
-		L_denom=[HSK_8[t[2]][t[0]] for t in L_ratios_ind]
-		L_denom_inv=batch_inversion(L_denom)
-		L_num=[HSK_8[t[2]][t[1]] for t in L_ratios_ind]
-		L_ratios=[L_num[i]*L_denom_inv[i] for i in range(15)]
-
-		L_coords_ind=tree_ratios.edge_product(L_ratios)
-
-		O_coords=[ZZ(0) for i in range(16)]
-		for t in L_coords_ind:
-			if self._coerce:
-				O_coords[t[1]]=self._coerce(t[0])
-			else:
-				O_coords[t[1]]=t[0]
-
-		# Precomputation
-		# TODO: optimize inversions and give precomputation to the codomain _arithmetic_precomputation
-		L_prec=[]
-		L_prec_ind=[]
-		for i in range(16):
-			if O_coords[i]!=0:
-				L_prec.append(O_coords[i])
-				L_prec_ind.append(i)
-		L_prec_inv=batch_inversion(L_prec)
-		precomputation=[None for i in range(16)]
-		for i in range(len(L_prec)):
-			precomputation[L_prec_ind[i]]=L_prec_inv[i]
-
-		self._precomputation=precomputation
-
-		for k in range(len(L_K_8)):
-			for j in range(16):
-				jpk=j^multindex_to_index(L_K_8_ind[k])
-				assert HSK_8[k][j]*O_coords[jpk]==HSK_8[k][jpk]*O_coords[j]
-		
-		assert proj_equal(squared(self._domain._null_point.coords()), hadamard(squared(O_coords)))
-
-		self._codomain=ThetaStructureDim4(hadamard(O_coords),null_point_dual=O_coords)
-
-	def special_image(self,P,L_trans,L_trans_ind):
-		r"""Used when we cannot evaluate the isogeny self because the codomain has zero 
-		dual theta constants.
-
-		Input:
-		- P: ThetaPointDim4 of the domain.
-		- L_trans: list of translates of P+T of P by points of 4-torsion T above the kernel.
-		- L_trans_ind: list of indices of the translation 4-torsion points T.
-		If L_trans[i]=\sum i_j*B_K4[j] then L_trans_ind[j]=\sum 2**j*i_j.
-
-		Output:
-		- the image of P by the isogeny self.
-		"""
-		HS_P=hadamard(squared(P.coords()))
-		HSL_trans=[hadamard(squared(Q.coords())) for Q in L_trans]
-		O_coords=self._codomain.null_point_dual()
-		
-		# L_lambda_inv: List of inverses of lambda_i such that: 
-		# HS(P+Ti)=(lambda_i*U_{chi.chi_i,0}(f(P))*U_{chi,0}(0))_chi.
-		L_lambda_inv_num=[]
-		L_lambda_inv_denom=[]
-
-		for k in range(len(L_trans)):
-			for j in range(16):
-				jpk=j^L_trans_ind[k]
-				if HSL_trans[k][j]!=0 and O_coords[jpk]!=0:
-					L_lambda_inv_num.append(HS_P[jpk]*O_coords[j])
-					L_lambda_inv_denom.append(HSL_trans[k][j]*O_coords[jpk])
-					break
-		L_lambda_inv_denom=batch_inversion(L_lambda_inv_denom)
-		L_lambda_inv=[L_lambda_inv_num[i]*L_lambda_inv_denom[i] for i in range(len(L_trans))]
-
-		for k in range(len(L_trans)):
-			for j in range(16):
-				jpk=j^L_trans_ind[k]
-				assert HS_P[jpk]*O_coords[j]==L_lambda_inv[k]*HSL_trans[k][j]*O_coords[jpk]
-
-		U_fP=[]
-		for i in range(16):
-			if self._precomputation[i]!=None:
-				U_fP.append(self._precomputation[i]*HS_P[i])
-			else:
-				for k in range(len(L_trans)):
-					ipk=i^L_trans_ind[k]
-					if self._precomputation[ipk]!=None:
-						U_fP.append(self._precomputation[ipk]*HSL_trans[k][ipk]*L_lambda_inv[k])
-						break
-
-		fP=hadamard(U_fP)
-		if self._coerce:
-			fP=[self._coerce(x) for x in fP]
-
-		return self._codomain(fP)
-
-	def dual(self):
-		return DualIsogenyDim4(self._codomain,self._domain, hadamard=False)
diff --git a/theta_lib/isogenies/isogeny_chain_dim4.py b/theta_lib/isogenies/isogeny_chain_dim4.py
deleted file mode 100644
index 8c0b78e..0000000
--- a/theta_lib/isogenies/isogeny_chain_dim4.py
+++ /dev/null
@@ -1,114 +0,0 @@
-from sage.all import *
-from ..utilities.strategy import precompute_strategy_with_first_eval, precompute_strategy_with_first_eval_and_splitting
-from .isogeny_dim4 import IsogenyDim4
-
-
-class IsogenyChainDim4:
-	def __init__(self, B_K, first_isogenies, e, m, splitting=True, strategy = None):
-		self.e=e
-		self.m=m
-
-		if strategy == None:
-			strategy = self.get_strategy(splitting)
-		self.strategy = strategy
-
-		self._isogenies=self.isogeny_chain(B_K, first_isogenies)
-
-
-	def get_strategy(self,splitting):
-		if splitting:
-			return precompute_strategy_with_first_eval_and_splitting(self.e,self.m,M=1,S=0.8,I=100)
-		else:
-			return precompute_strategy_with_first_eval(self.e,self.m,M=1,S=0.8,I=100)
-
-	def isogeny_chain(self, B_K, first_isogenies):
-		"""
-		Compute the isogeny chain and store intermediate isogenies for evaluation
-		"""
-		# Store chain of (2,2)-isogenies
-		isogeny_chain = []
-
-		# Bookkeeping for optimal strategy
-		strat_idx = 0
-		level = [0]
-		ker = B_K
-		kernel_elements = [ker]
-
-		# Length of the chain
-		n=self.e-self.m
-		
-		for k in range(n):
-			prev = sum(level)
-			ker = kernel_elements[-1]
-
-			while prev != (n - 1 - k):
-				level.append(self.strategy[strat_idx])
-				prev += self.strategy[strat_idx]
-
-				# Perform the doublings and update kernel elements
-				# Prevent the last unnecessary doublings for first isogeny computation
-				if k>0 or prev!=n-1:
-					ker = [ker[i].double_iter(self.strategy[strat_idx]) for i in range(4)]
-					kernel_elements.append(ker)
-
-				# Update bookkeeping variable
-				strat_idx += 1
-
-			# Compute the codomain from the 8-torsion
-			if k==0:
-				phi = first_isogenies
-			else:
-				phi = IsogenyDim4(Th,ker)
-
-			# Update the chain of isogenies
-			Th = phi._codomain
-			# print(parent(Th.null_point().coords()[0]))
-			isogeny_chain.append(phi)
-
-			# Remove elements from list
-			if k>0:
-				kernel_elements.pop()
-			level.pop()
-
-			# Push through points for the next step
-			kernel_elements = [[phi(T) for T in kernel] for kernel in kernel_elements]
-			# print([[parent(T.coords()[0]) for T in kernel] for kernel in kernel_elements])
-
-		return isogeny_chain
-
-	def evaluate_isogeny(self,P):
-		Q=P
-		for f in self._isogenies:
-			Q=f(Q)
-		return Q
-
-	def __call__(self,P):
-		return self.evaluate_isogeny(P)
-
-	def dual(self):
-		n=len(self._isogenies)
-		isogenies=[]
-		for i in range(n):
-			isogenies.append(self._isogenies[n-1-i].dual())
-		return DualIsogenyChainDim4(isogenies)
-
-
-class DualIsogenyChainDim4:
-	def __init__(self,isogenies):
-		self._isogenies=isogenies
-
-	def evaluate_isogeny(self,P):
-		n=len(self._isogenies)
-		Q=P
-		for j in range(n):
-			Q=self._isogenies[j](Q)
-		return Q
-
-	def __call__(self,P):
-		return self.evaluate_isogeny(P)
-
-
-
-
-		
-
diff --git a/theta_lib/isogenies/isogeny_dim4.py b/theta_lib/isogenies/isogeny_dim4.py
deleted file mode 100644
index 2f483bf..0000000
--- a/theta_lib/isogenies/isogeny_dim4.py
+++ /dev/null
@@ -1,162 +0,0 @@
-from sage.all import *
-
-from ..theta_structures.Theta_dim4 import ThetaStructureDim4
-from ..theta_structures.theta_helpers_dim4 import hadamard, squared, batch_inversion
-from .tree import Tree
-
-class IsogenyDim4:
-	def __init__(self,domain,K_8,codomain=None,precomputation=None):
-		r"""
-		Input:
-		- domain: a ThetaStructureDim4.
-		- K_8: a list of 4 points of 8-torision (such that 4*K_8 is a kernel basis), used to compute the codomain.
-		- codomain: a ThetaStructureDim4 (for the codomain, used only when K_8 is None).
-		- precomputation: list of inverse of dual theta constants of the codomain, used to compute the image.
-		"""
-
-		if not isinstance(domain, ThetaStructureDim4):
-			raise ValueError("Argument domain should be a ThetaStructureDim4 object.")
-		self._domain = domain
-		self._precomputation=None
-		if K_8!=None:
-			self._compute_codomain(K_8)
-		else:
-			self._codomain=codomain
-			self._precomputation=precomputation
-
-	def _compute_codomain(self,K_8):
-		r"""
-		Input:
-		- K_8: a list of 4 points of 8-torision (such that 4*K_8 is a kernel basis).
-
-		Output:
-		- codomain of the isogeny.
-		Also initializes self._precomputation, containing the inverse of theta-constants.
-		"""
-		HSK_8=[hadamard(squared(P.coords())) for P in K_8]
-		
-		# Choice of reference index j_0<->chi_0 corresponding to a non-vanishing theta-constant.
-		found_tree=False
-		j_0=0
-		while not found_tree:
-			found_k0=False
-			for k in range(4):
-				if j_0>15:
-					raise NotImplementedError("The codomain of this 2-isogeny could not be computed.\nWe may have encountered a product of abelian varieties\nsomewhere unexpected along the chain.\nThis is exceptionnal and should not happen in larger characteristic.")
-				if HSK_8[k][j_0]!=0:
-					k_0=k
-					found_k0=True
-					break
-			if not found_k0:
-				j_0+=1
-			else:
-				j0pk0=j_0^(2**k_0)
-				# List of tuples of indices (index chi of the denominator: HS(f(P_k))_chi, 
-				#index chi.chi_k of the numerator: HS(f(P_k))_chi.chi_k, index k).
-				L_ratios_ind=[(j_0,j0pk0,k_0)]
-				L_covered_ind=[j_0,j0pk0]
-
-				# Tree containing the the theta-null points indices as nodes and the L_ratios_ind reference indices as edges.
-				tree_ratios=Tree(j_0)
-				tree_ratios.add_child(Tree(j0pk0),k_0)
-
-				# Filling in the tree
-				tree_filled=False
-				while not tree_filled:
-					found_j=False
-					for j in L_covered_ind:
-						for k in range(4):
-							jpk=j^(2**k)
-							if jpk not in L_covered_ind and HSK_8[k][j]!=0:
-								L_covered_ind.append(jpk)
-								L_ratios_ind.append((j,jpk,k))
-								tree_j=tree_ratios.look_node(j)
-								tree_j.add_child(Tree(jpk),len(L_ratios_ind)-1)
-								found_j=True
-								break
-						if found_j:
-							break
-					if not found_j or len(L_covered_ind)==16:
-						tree_filled=True
-				if len(L_covered_ind)!=16:
-					j_0+=1
-				else:
-					found_tree=True
-
-		L_denom=[HSK_8[t[2]][t[0]] for t in L_ratios_ind]
-		L_denom_inv=batch_inversion(L_denom)
-		L_num=[HSK_8[t[2]][t[1]] for t in L_ratios_ind]
-		L_ratios=[L_num[i]*L_denom_inv[i] for i in range(15)]
-
-		L_coords_ind=tree_ratios.edge_product(L_ratios)
-
-		O_coords=[ZZ(0) for i in range(16)]
-		for t in L_coords_ind:
-			O_coords[t[1]]=t[0]
-
-		# Precomputation
-		# TODO: optimize inversions
-		L_prec=[]
-		L_prec_ind=[]
-		for i in range(16):
-			if O_coords[i]!=0:
-				L_prec.append(O_coords[i])
-				L_prec_ind.append(i)
-		L_prec_inv=batch_inversion(L_prec)
-		precomputation=[None for i in range(16)]
-		for i in range(len(L_prec)):
-			precomputation[L_prec_ind[i]]=L_prec_inv[i]
-
-		self._precomputation=precomputation
-		# Assumes there is no zero theta constant. Otherwise, squared(precomputation) will raise an error (None**2 does not exist)
-		self._codomain=ThetaStructureDim4(hadamard(O_coords),null_point_dual=O_coords)
-
-	def codomain(self):
-		return self._codomain
-
-	def domain(self):
-		return self._domain
-
-	def image(self,P):
-		HS_P=list(hadamard(squared(P.coords())))
-
-		for i in range(16):
-			HS_P[i] *=self._precomputation[i]
-
-		return self._codomain(hadamard(HS_P))
-
-	def dual(self):
-		return DualIsogenyDim4(self._codomain,self._domain, hadamard=True)
-
-	def __call__(self,P):
-		return self.image(P)
-
-
-class DualIsogenyDim4:
-	def __init__(self,domain,codomain,hadamard=True):
-		# domain and codomain are respectively the domain and codomain of \tilde{f}: domain-->codomain,
-		# so respectively  the codomain and domain of f: codomain-->domain.
-		# By convention, domain input is given in usual coordinates (ker(\tilde{f})=K_2).
-		# codomain is in usual coordinates if hadamard, in dual coordinates otherwise.
-		self._domain=domain.hadamard()
-		self._hadamard=hadamard
-		if hadamard:
-			self._codomain=codomain.hadamard()
-			self._precomputation=batch_inversion(codomain.zero().coords())
-		else:
-			self._codomain=codomain
-			self._precomputation=batch_inversion(codomain.zero().coords())
-
-	def image(self,P):
-		# When ker(f)=K_2, ker(\tilde{f})=K_1 so ker(\tilde{f})=K_2 after hadamard transformation of the 
-		# new domain (ex codomain)
-		HS_P=list(hadamard(squared(P.coords())))
-		for i in range(16):
-			HS_P[i] *=self._precomputation[i]
-		if self._hadamard:
-			return self._codomain(hadamard(HS_P))
-		else:
-			return self._codomain(HS_P)
-
-	def __call__(self,P):
-		return self.image(P)
diff --git a/theta_lib/isogenies/tree.py b/theta_lib/isogenies/tree.py
deleted file mode 100644
index a6e3da3..0000000
--- a/theta_lib/isogenies/tree.py
+++ /dev/null
@@ -1,28 +0,0 @@
-from sage.all import *
-
-class Tree:
-	def __init__(self,node):
-		self._node=node
-		self._edges=[]
-		self._children=[]
-
-	def add_child(self,child,edge):
-		self._children.append(child)
-		self._edges.append(edge)
-
-	def look_node(self,node):
-		if self._node==node:
-			return self
-		elif len(self._children)>0:
-			for child in self._children:
-				t_node=child.look_node(node)
-				if t_node!=None:
-					return t_node
-
-	def edge_product(self,L_factors,factor_node=ZZ(1)):
-		n=len(self._children)
-		L_prod=[(factor_node,self._node)]
-		for i in range(n):
-			L_prod+=self._children[i].edge_product(L_factors,factor_node*L_factors[self._edges[i]])
-		return L_prod
-
diff --git a/theta_lib/isogenies_dim2/gluing_isogeny_dim2.py b/theta_lib/isogenies_dim2/gluing_isogeny_dim2.py
deleted file mode 100644
index 2dac387..0000000
--- a/theta_lib/isogenies_dim2/gluing_isogeny_dim2.py
+++ /dev/null
@@ -1,292 +0,0 @@
-"""
-This code is based on a copy of:
-https://github.com/ThetaIsogenies/two-isogenies
-
-MIT License
-
-Copyright (c) 2023 Pierrick Dartois, Luciano Maino, Giacomo Pope and Damien Robert
-
-Permission is hereby granted, free of charge, to any person obtaining a copy
-of this software and associated documentation files (the "Software"), to deal
-in the Software without restriction, including without limitation the rights
-to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-copies of the Software, and to permit persons to whom the Software is
-furnished to do so, subject to the following conditions:
-
-The above copyright notice and this permission notice shall be included in all
-copies or substantial portions of the Software.
-
-THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
-SOFTWARE.
-"""
-
-from sage.all import *
-from ..theta_structures.Tuple_point import TuplePoint
-from ..theta_structures.Theta_dim2 import ThetaStructureDim2, ThetaPointDim2
-from ..theta_structures.theta_helpers_dim2 import batch_inversion
-from ..basis_change.base_change_dim2 import montgomery_to_theta_matrix_dim2, apply_base_change_theta_dim2
-from ..theta_structures.montgomery_theta import lift_kummer_montgomery_point
-
-class GluingThetaIsogenyDim2:
-    """
-    Compute the gluing isogeny from E1 x E2 (Elliptic Product) -> A (Theta Model)
-
-    Expected input:
-
-    - (K1_8, K2_8) The 8-torsion above the kernel generating the isogeny
-    - M (Optional) a base change matrix, if this is not including, it can
-      be derived from [2](K1_8, K2_8)
-    """
-
-    def __init__(self, K1_8, K2_8, Theta12, N):
-        # Double points to get four-torsion, we always need one of these, used
-        # for the image computations but we'll need both if we wish to derived
-        # the base change matrix as well
-        K1_4 = 2*K1_8
-
-        # Initalise self
-        # This is the base change matrix for product Theta coordinates (not used, except in the dual)
-        self._base_change_matrix_theta = N
-        # Here, base change matrix directly applied to the Montgomery coordinates. null_point_bc is the 
-        # theta null point obtained after applying the base change to the product Theta-structure.
-        self._base_change_matrix, null_point_bc = montgomery_to_theta_matrix_dim2(Theta12.zero().coords(),N, return_null_point = True)
-        self._domain_bc = ThetaStructureDim2(null_point_bc)
-        self.T_shift = K1_4
-        self._precomputation = None
-        self._zero_idx = 0
-        self._domain_product = Theta12
-        self._domain=(K1_8[0].curve(), K1_8[1].curve())
-
-        # Map points from elliptic product onto the product theta structure
-        # using the base change matrix
-        T1_8 = self.base_change(K1_8)
-        T2_8 = self.base_change(K2_8)
-
-        # Compute the codomain of the gluing isogeny
-        self._codomain = self._special_compute_codomain(T1_8, T2_8)
-
-    def apply_base_change(self, coords):
-        """
-        Apply the basis change by acting with matrix multiplication, treating
-        the coordinates as a vector
-        """
-        N = self._base_change_matrix
-        x, y, z, t = coords
-        X = N[0, 0] * x + N[0, 1] * y + N[0, 2] * z + N[0, 3] * t
-        Y = N[1, 0] * x + N[1, 1] * y + N[1, 2] * z + N[1, 3] * t
-        Z = N[2, 0] * x + N[2, 1] * y + N[2, 2] * z + N[2, 3] * t
-        T = N[3, 0] * x + N[3, 1] * y + N[3, 2] * z + N[3, 3] * t
-
-        return (X, Y, Z, T)
-
-    def base_change(self, P):
-        """
-        Compute the basis change on a TuplePoint to recover a ThetaPointDim2 of
-        compatible form
-        """
-        if not isinstance(P, TuplePoint):
-            raise TypeError("Function assumes that the input is of type `TuplePoint`")
-
-        # extract X,Z coordinates on pairs of points
-        P1, P2 = P.points()
-        X1, Z1 = P1[0], P1[2]
-        X2, Z2 = P2[0], P2[2]
-
-        # Correct in the case of (0 : 0)
-        if X1 == 0 and Z1 == 0:
-            X1 = 1
-            Z1 = 0
-        if X2 == 0 and Z2 == 0:
-            X2 = 1
-            Z2 = 0
-
-        # Apply the basis transformation on the product
-        coords = self.apply_base_change([X1 * X2, X1 * Z2, Z1 * X2, Z1 * Z2])
-        return coords
-
-    def _special_compute_codomain(self, T1, T2):
-        """
-        Given twzero_matro isotropic points of 8-torsion T1 and T2, compatible with
-        the theta null point, compute the level two theta null point A/K_2
-        """
-        xAxByCyD = ThetaPointDim2.to_squared_theta(*T1)
-        zAtBzYtD = ThetaPointDim2.to_squared_theta(*T2)
-
-        # Find the value of the non-zero index
-        zero_idx = next((i for i, x in enumerate(xAxByCyD) if x == 0), None)
-        self._zero_idx = zero_idx
-
-        # Dumb check to make sure everything is OK
-        assert xAxByCyD[self._zero_idx] == zAtBzYtD[self._zero_idx] == 0
-
-        # Initialize lists
-        # The zero index described the permutation
-        ABCD = [0 for _ in range(4)]
-        precomp = [0 for _ in range(4)]
-
-        # Compute non-trivial numerators (Others are either 1 or 0)
-        num_1 = zAtBzYtD[1 ^ self._zero_idx]
-        num_2 = xAxByCyD[2 ^ self._zero_idx]
-        num_3 = zAtBzYtD[3 ^ self._zero_idx]
-        num_4 = xAxByCyD[3 ^ self._zero_idx]
-
-        # Compute and invert non-trivial denominators
-        den_1, den_2, den_3, den_4 = batch_inversion([num_1, num_2, num_3, num_4])
-
-        # Compute A, B, C, D
-        ABCD[0 ^ self._zero_idx] = 0
-        ABCD[1 ^ self._zero_idx] = num_1 * den_3
-        ABCD[2 ^ self._zero_idx] = num_2 * den_4
-        ABCD[3 ^ self._zero_idx] = 1
-
-        # Compute precomputation for isogeny images
-        precomp[0 ^ self._zero_idx] = 0
-        precomp[1 ^ self._zero_idx] = den_1 * num_3
-        precomp[2 ^ self._zero_idx] = den_2 * num_4
-        precomp[3 ^ self._zero_idx] = 1
-        self._precomputation = precomp
-
-        # Final Hadamard of the above coordinates
-        a, b, c, d = ThetaPointDim2.to_hadamard(*ABCD)
-
-        return ThetaStructureDim2([a, b, c, d])
-
-    def special_image(self, P, translate):
-        """
-        When the domain is a non product theta structure on a product of
-        elliptic curves, we will have one of A,B,C,D=0, so the image is more
-        difficult. We need to give the coordinates of P but also of
-        P+Ti, Ti one of the point of 4-torsion used in the isogeny
-        normalisation
-        """
-        AxByCzDt = ThetaPointDim2.to_squared_theta(*P)
-
-        # We are in the case where at most one of A, B, C, D is
-        # zero, so we need to account for this
-        #
-        # To recover values, we use the translated point to get
-        AyBxCtDz = ThetaPointDim2.to_squared_theta(*translate)
-
-        # Directly compute y,z,t
-        y = AxByCzDt[1 ^ self._zero_idx] * self._precomputation[1 ^ self._zero_idx]
-        z = AxByCzDt[2 ^ self._zero_idx] * self._precomputation[2 ^ self._zero_idx]
-        t = AxByCzDt[3 ^ self._zero_idx]
-
-        # We can compute x from the translation
-        # First we need a normalisation
-        if z != 0:
-            zb = AyBxCtDz[3 ^ self._zero_idx]
-            lam = z / zb
-        else:
-            tb = AyBxCtDz[2 ^ self._zero_idx] * self._precomputation[2 ^ self._zero_idx]
-            lam = t / tb
-
-        # Finally we recover x
-        xb = AyBxCtDz[1 ^ self._zero_idx] * self._precomputation[1 ^ self._zero_idx]
-        x = xb * lam
-
-        xyzt = [0 for _ in range(4)]
-        xyzt[0 ^ self._zero_idx] = x
-        xyzt[1 ^ self._zero_idx] = y
-        xyzt[2 ^ self._zero_idx] = z
-        xyzt[3 ^ self._zero_idx] = t
-
-        image = ThetaPointDim2.to_hadamard(*xyzt)
-        return self._codomain(image)
-
-    def __call__(self, P):
-        """
-        Take into input the theta null point of A/K_2, and return the image
-        of the point by the isogeny
-        """
-        if not isinstance(P, TuplePoint):
-            raise TypeError(
-                "Isogeny image for the gluing isogeny is defined to act on TuplePoints"
-            )
-
-        # Compute sum of points on elliptic curve
-        P_sum_T = P + self.T_shift
-
-        # Push both the point and the translation through the
-        # completion
-        iso_P = self.base_change(P)
-        iso_P_sum_T = self.base_change(P_sum_T)
-
-        return self.special_image(iso_P, iso_P_sum_T)
-
-    def dual(self):
-        domain = self._codomain.hadamard()
-        codomain_bc = self._domain_bc.hadamard()
-        codomain = self._domain
-
-        precomputation = batch_inversion(codomain_bc.null_point_dual())
-
-        N_split = self._base_change_matrix.inverse()
-
-        return DualGluingThetaIsogenyDim2(domain, codomain_bc, codomain, N_split, precomputation)
-
-
-class DualGluingThetaIsogenyDim2:
-    def __init__(self, domain, codomain_bc, codomain, N_split, precomputation):
-        self._domain = domain
-        self._codomain_bc = codomain_bc # Theta structure
-        self._codomain = codomain # Elliptic curves E1 and E2
-        self._precomputation = precomputation
-        self._splitting_matrix = N_split
-
-    def __call__(self,P):
-        # Returns a TuplePoint.
-        if not isinstance(P, ThetaPointDim2):
-            raise TypeError("Isogeny evaluation expects a ThetaPointDim2 as input")
-
-        xx, yy, zz, tt = P.squared_theta()
-
-        Ai, Bi, Ci, Di = self._precomputation
-
-        xx = xx * Ai
-        yy = yy * Bi
-        zz = zz * Ci
-        tt = tt * Di
-
-        image_coords = (xx, yy, zz, tt)
-
-        X1X2, X1Z2, Z1X2, Z1Z2 = apply_base_change_theta_dim2(self._splitting_matrix, image_coords)
-
-        E1, E2 = self._codomain
-
-        if Z1Z2!=0:
-            #Z1=1, Z2=Z1Z2
-
-            Z2_inv=1/Z1Z2
-            X2=Z1X2*Z2_inv# Normalize (X2:Z2)=(X2/Z2:1)
-
-            X1=X1Z2*Z2_inv
-
-            assert X1*Z1X2==X1X2
-            P1 = lift_kummer_montgomery_point(E1, X1)
-            P2 = lift_kummer_montgomery_point(E2, X2)
-            return TuplePoint(P1,P2)
-        elif Z1X2==0 and X1Z2!=0:
-            # Case (X1:Z1)=0, X1!=0 and (X2:Z2)!=0
-
-            X2=X1X2/X1Z2
-            P2 = lift_kummer_montgomery_point(E2, X2)
-            return TuplePoint(E1(0),P2)
-        elif Z1X2!=0 and X1Z2==0:
-            # Case (X1:Z1)!=0 and (X2:Z2)=0, X2!=0
-
-            X1=X1X2/Z1X2
-            P1 = lift_kummer_montgomery_point(E1, X1)
-            return TuplePoint(P1,E2(0))
-        else:
-            return TuplePoint(E1(0),E2(0))
-
-
-
-
-
diff --git a/theta_lib/isogenies_dim2/isogeny_chain_dim2.py b/theta_lib/isogenies_dim2/isogeny_chain_dim2.py
deleted file mode 100644
index 23adb23..0000000
--- a/theta_lib/isogenies_dim2/isogeny_chain_dim2.py
+++ /dev/null
@@ -1,178 +0,0 @@
-"""
-This code is based on a copy of:
-https://github.com/ThetaIsogenies/two-isogenies
-
-MIT License
-
-Copyright (c) 2023 Pierrick Dartois, Luciano Maino, Giacomo Pope and Damien Robert
-
-Permission is hereby granted, free of charge, to any person obtaining a copy
-of this software and associated documentation files (the "Software"), to deal
-in the Software without restriction, including without limitation the rights
-to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-copies of the Software, and to permit persons to whom the Software is
-furnished to do so, subject to the following conditions:
-
-The above copyright notice and this permission notice shall be included in all
-copies or substantial portions of the Software.
-
-THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
-SOFTWARE.
-"""
-
-from sage.all import *
-from ..theta_structures.Tuple_point import TuplePoint
-from ..theta_structures.Theta_dim2 import ThetaPointDim2
-from .gluing_isogeny_dim2 import GluingThetaIsogenyDim2
-from .isogeny_dim2 import ThetaIsogenyDim2
-from ..utilities.strategy import optimised_strategy
-
-
-class IsogenyChainDim2:
-    r"""
-    Given (P1, P2), (Q1, Q2) in (E1 x E2)[2^(n+2)] as the generators of a kernel
-    of a (2^n, 2^n)-isogeny
-
-    ker(Phi) = <(P1, P2), (Q1, Q2)>
-
-    Input:
-
-    - kernel = TuplePoint(P1, P2), TuplePoint(Q1, Q2):
-      where points are on the elliptic curves E1, E2 of order 2^(n+2)
-    - n: the length of the chain
-    - strategy: the optimises strategy to compute a walk through the graph of
-      images and doublings with a quasli-linear number of steps
-    """
-
-    def __init__(self, kernel, Theta12, M, n, strategy=None):
-        self.n = n
-        self.E1, self.E2 = kernel[0].parent_curves()
-        assert kernel[1].parent_curves() == [self.E1, self.E2]
-
-        self._domain = (self.E1, self.E2)
-
-        if strategy is None:
-            strategy = self.get_strategy()
-        self.strategy = strategy
-
-        self._phis = self.isogeny_chain(kernel, Theta12, M)
-
-        self._codomain=self._phis[-1]._codomain
-
-    def get_strategy(self):
-        return optimised_strategy(self.n)
-
-    def isogeny_chain(self, kernel, Theta12, M):
-        """
-        Compute the isogeny chain and store intermediate isogenies for evaluation
-        """
-        # Extract the CouplePoints from the Kernel
-        Tp1, Tp2 = kernel
-
-        # Store chain of (2,2)-isogenies
-        isogeny_chain = []
-
-        # Bookkeeping for optimal strategy
-        strat_idx = 0
-        level = [0]
-        ker = (Tp1, Tp2)
-        kernel_elements = [ker]
-
-        for k in range(self.n):
-            prev = sum(level)
-            ker = kernel_elements[-1]
-
-            while prev != (self.n - 1 - k):
-                level.append(self.strategy[strat_idx])
-
-                # Perform the doublings
-                Tp1 = ker[0].double_iter(self.strategy[strat_idx])
-                Tp2 = ker[1].double_iter(self.strategy[strat_idx])
-
-                ker = (Tp1, Tp2)
-
-                # Update kernel elements and bookkeeping variables
-                kernel_elements.append(ker)
-                prev += self.strategy[strat_idx]
-                strat_idx += 1
-
-            # Compute the codomain from the 8-torsion
-            Tp1, Tp2 = ker
-            if k == 0:
-                phi = GluingThetaIsogenyDim2(Tp1, Tp2, Theta12, M)
-            else:
-                phi = ThetaIsogenyDim2(Th, Tp1, Tp2)
-
-            # Update the chain of isogenies
-            Th = phi._codomain
-            isogeny_chain.append(phi)
-
-            # Remove elements from list
-            kernel_elements.pop()
-            level.pop()
-
-            # Push through points for the next step
-            kernel_elements = [(phi(T1), phi(T2)) for T1, T2 in kernel_elements]
-
-        return isogeny_chain
-
-    def evaluate_isogeny(self, P):
-        """
-        Given a point P, of type TuplePoint on the domain E1 x E2, computes the
-        ThetaPointDim2 on the codomain ThetaStructureDim2.
-        """
-        if not isinstance(P, TuplePoint):
-            raise TypeError(
-                "IsogenyChainDim2 isogeny expects as input a TuplePoint on the domain product E1 x E2"
-            )
-        n=len(self._phis)
-        for i in range(n):
-            P = self._phis[i](P)
-        return P
-
-    def __call__(self, P):
-        """
-        Evaluate a TuplePoint under the action of this isogeny.
-        """
-        return self.evaluate_isogeny(P)
-
-    def dual(self):
-        domain = self._codomain
-        codomain = self._domain
-        n=len(self._phis)
-        isogenies=[]
-        for i in range(n):
-            isogenies.append(self._phis[n-1-i].dual())
-        return DualIsogenyChainDim2(domain, codomain, isogenies)
-
-
-class DualIsogenyChainDim2:
-    def __init__(self, domain, codomain, isogenies):
-        self._domain = domain
-        self._codomain = codomain
-        self._phis = isogenies
-
-    def evaluate_isogeny(self, P):
-        """
-        Given a ThetaPointDim2 point P on the codomain ThetaStructureDim2, 
-        computes the image TuplePoint on the codomain E1 x E2.
-        """
-        if not isinstance(P, ThetaPointDim2):
-            raise TypeError(
-                "DualIsogenyChainDim2 isogeny expects as input a ThetaPointDim2."
-            )
-        n=len(self._phis)
-        for i in range(n):
-            P = self._phis[i](P)
-        return P
-
-    def __call__(self, P):
-        """
-        Evaluate a ThetaPointDim2 under the action of this isogeny.
-        """
-        return self.evaluate_isogeny(P)
diff --git a/theta_lib/isogenies_dim2/isogeny_dim2.py b/theta_lib/isogenies_dim2/isogeny_dim2.py
deleted file mode 100644
index b740ab1..0000000
--- a/theta_lib/isogenies_dim2/isogeny_dim2.py
+++ /dev/null
@@ -1,229 +0,0 @@
-"""
-This code is based on a copy of:
-https://github.com/ThetaIsogenies/two-isogenies
-
-MIT License
-
-Copyright (c) 2023 Pierrick Dartois, Luciano Maino, Giacomo Pope and Damien Robert
-
-Permission is hereby granted, free of charge, to any person obtaining a copy
-of this software and associated documentation files (the "Software"), to deal
-in the Software without restriction, including without limitation the rights
-to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
-copies of the Software, and to permit persons to whom the Software is
-furnished to do so, subject to the following conditions:
-
-The above copyright notice and this permission notice shall be included in all
-copies or substantial portions of the Software.
-
-THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
-IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
-FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
-AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
-LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
-OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
-SOFTWARE.
-"""
-
-from sage.all import ZZ
-from ..theta_structures.Theta_dim2 import ThetaStructureDim2, ThetaPointDim2
-from ..theta_structures.theta_helpers_dim2 import batch_inversion
-
-
-class ThetaIsogenyDim2:
-    def __init__(self, domain, T1_8, T2_8, hadamard=(False, True)):
-        """
-        Compute a (2,2)-isogeny in the theta model. Expects as input:
-
-        - domain: the ThetaStructureDim2 from which we compute the isogeny
-        - (T1_8, T2_8): points of 8-torsion above the kernel generating the isogeny
-
-        When the 8-torsion is not available (for example at the end of a long
-        (2,2)-isogeny chain), the the helper functions in isogeny_sqrt.py
-        must be used.
-
-        NOTE: on the hadamard bools:
-
-        The optional parameter 'hadamard' controls if we are in standard or dual
-        coordinates, and if the codomain is in standard or dual coordinates. By
-        default this is (False, True), meaning we use standard coordinates on
-        the domain A and the codomain B.
-
-        The kernel is then the kernel K_2 where the action is by sign. Other
-        possibilities: - (False, False): standard coordinates on A, dual
-        coordinates on B - (True, True): start in dual coordinates on A
-        (alternatively: standard coordinates on A but quotient by K_1 whose
-        action is by permutation), and standard coordinates on B. - (True,
-        False): dual coordinates on A and B
-
-        These can be composed as follows for A -> B -> C:
-
-        - (False, True) -> (False, True) (False, False) -> (True, True):
-          - standard coordinates on A and C,
-          - standard/resp dual coordinates on B
-        - (False, True) -> (False, False) (False, False) -> (True, False):
-          - standard coordinates on A,
-          - dual coordinates on C,
-          - standard/resp dual coordinates on B
-        - (True, True) -> (False, True) (True, False) -> (True, True):
-          - dual coordinates on A,
-          - standard coordinates on C,
-          - standard/resp dual coordiantes on B
-        - (True, True) -> (False, False) (True, False) -> (True, False):
-          - dual coordinates on A and C
-          - standard/resp dual coordinates on B
-
-        On the other hand, these gives the multiplication by [2] on A:
-
-        - (False, False) -> (False, True) (False, True) -> (True, True):
-          - doubling in standard coordinates on A
-          - going through dual/standard coordinates on B=A/K_2
-        - (True, False) -> (False, False) (True, True) -> (True, False):
-          - doubling in dual coordinates on A
-          - going through dual/standard coordinates on B=A/K_2
-            (alternatively: doubling in standard coordinates on A going
-            through B'=A/K_1)
-        - (False, False) -> (False, False) (False, True) -> (True, False):
-          - doubling from standard to dual coordinates on A
-        - (True, False) -> (False, True) (True, True) -> (True, True):
-          - doubling from dual to standard coordinates on A
-        """
-        if not isinstance(domain, ThetaStructureDim2):
-            raise ValueError
-        self._domain = domain
-
-        self._hadamard = hadamard
-        self._precomputation = None
-        self._codomain = self._compute_codomain(T1_8, T2_8)
-
-    def _compute_codomain(self, T1, T2):
-        """
-        Given two isotropic points of 8-torsion T1 and T2, compatible with
-        the theta null point, compute the level two theta null point A/K_2
-        """
-        if self._hadamard[0]:
-            xA, xB, _, _ = ThetaPointDim2.to_squared_theta(
-                *ThetaPointDim2.to_hadamard(*T1.coords())
-            )
-            zA, tB, zC, tD = ThetaPointDim2.to_squared_theta(
-                *ThetaPointDim2.to_hadamard(*T2.coords())
-            )
-        else:
-            xA, xB, _, _ = T1.squared_theta()
-            zA, tB, zC, tD = T2.squared_theta()
-
-        if not self._hadamard[0] and self._domain._precomputation:
-            # Batch invert denominators
-            xA_inv, zA_inv, tB_inv = batch_inversion([xA, zA, tB])
-
-            # Compute A, B, C, D
-            A = ZZ(1)
-            B = xB * xA_inv
-            C = zC * zA_inv
-            D = tD * tB_inv * B
-
-            _, _, _, BBinv, CCinv, DDinv = self._domain._arithmetic_precomputation()
-            B_inv = BBinv * B
-            C_inv = CCinv * C
-            D_inv = DDinv * D
-        else:
-            # Batch invert denominators
-            xA_inv, zA_inv, tB_inv, xB_inv, zC_inv, tD_inv = batch_inversion([xA, zA, tB, xB, zC, tD])
-
-            # Compute A, B, C, D
-            A = ZZ(1)
-            B = xB * xA_inv
-            C = zC * zA_inv
-            D = tD * tB_inv * B
-            B_inv = xB_inv * xA
-            C_inv = zC_inv * zA
-            D_inv = tD_inv * tB * B_inv
-
-            # NOTE: some of the computations we did here could be reused for the
-            # arithmetic precomputations of the codomain However, we are always
-            # in the mode (False, True) except the very last isogeny, so we do
-            # not lose much by not doing this optimisation Just in case we need
-            # it later:
-            # - for hadamard=(False, True): we can reuse the arithmetic
-            #   precomputation; we do this already above
-            # - for hadamard=(False, False): we can reuse the arithmetic
-            #   precomputation as above, and furthermore we could reuse B_inv,
-            #   C_inv, D_inv for the precomputation of the codomain
-            # - for hadamard=(True, False): we could reuse B_inv, C_inv, D_inv
-            #   for the precomputation of the codomain
-            # - for hadamard=(True, True): nothing to reuse!
-
-        self._precomputation = (B_inv, C_inv, D_inv)
-        if self._hadamard[1]:
-            a, b, c, d = ThetaPointDim2.to_hadamard(A, B, C, D)
-            return ThetaStructureDim2([a, b, c, d], null_point_dual=[A, B, C, D])
-        else:
-            return ThetaStructureDim2([A, B, C, D])
-
-    def __call__(self, P):
-        """
-        Take into input the theta null point of A/K_2, and return the image
-        of the point by the isogeny
-        """
-        if not isinstance(P, ThetaPointDim2):
-            raise TypeError("Isogeny evaluation expects a ThetaPointDim2 as input")
-
-        if self._hadamard[0]:
-            xx, yy, zz, tt = ThetaPointDim2.to_squared_theta(
-                *ThetaPointDim2.to_hadamard(*P.coords())
-            )
-        else:
-            xx, yy, zz, tt = P.squared_theta()
-
-        Bi, Ci, Di = self._precomputation
-
-        yy = yy * Bi
-        zz = zz * Ci
-        tt = tt * Di
-
-        image_coords = (xx, yy, zz, tt)
-        if self._hadamard[1]:
-            image_coords = ThetaPointDim2.to_hadamard(*image_coords)
-        return self._codomain(image_coords)
-
-    def dual(self):
-        # Returns the dual isogeny (domain and codomain are inverted).
-        # By convention, the new domain and codomain are in standard coordinates.
-        if self._hadamard[1]:
-            domain=self._codomain.hadamard()
-        else:
-            domain=self._codomain
-        if self._hadamard[0]:
-            codomain=self._domain
-        else:
-            codomain=self._domain.hadamard()
-
-        precomputation = batch_inversion(self._domain.null_point().coords())
-
-        return DualThetaIsogenyDim2(domain,codomain,precomputation)
-
-class DualThetaIsogenyDim2:
-    def __init__(self,domain,codomain,precomputation):
-        self._domain=domain
-        self._codomain=codomain
-        self._precomputation=precomputation
-
-    def __call__(self,P):
-        if not isinstance(P, ThetaPointDim2):
-            raise TypeError("Isogeny evaluation expects a ThetaPointDim2 as input")
-
-        xx, yy, zz, tt = P.squared_theta()
-
-        Ai, Bi, Ci, Di = self._precomputation
-
-        xx = xx * Ai
-        yy = yy * Bi
-        zz = zz * Ci
-        tt = tt * Di
-
-        image_coords = (xx, yy, zz, tt)
-        image_coords = ThetaPointDim2.to_hadamard(*image_coords)
-
-        return self._codomain(image_coords)
-        
-
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
-
-
-
-
diff --git a/theta_lib/utilities/discrete_log.py b/theta_lib/utilities/discrete_log.py
deleted file mode 100644
index dff04dc..0000000
--- a/theta_lib/utilities/discrete_log.py
+++ /dev/null
@@ -1,273 +0,0 @@
-# ===================================== #
-#  Fast DLP solving using Weil pairing  #
-# ===================================== #
-
-"""
-This code has been taken from:
-https://github.com/FESTA-PKE/FESTA-SageMath
-
-Copyright (c) 2023 Andrea Basso, Luciano Maino and Giacomo Pope.
-"""
-
-
-
-# Sage imports
-from sage.all import (
-    ZZ
-)
-
-# import pari for fast dlog
-import cypari2
-
-# Make instance of Pari
-pari = cypari2.Pari()
-
-def discrete_log_pari(a, base, order):
-    """
-    Wrapper around pari discrete log. Works like a.log(b), 
-    but allows us to use the optional argument order. This
-    is important as we skip Pari attempting to factor the
-    full order of Fp^k, which is slow.
-    """
-    x = pari.fflog(a, base, order)
-    return ZZ(x)
-
-def ell_discrete_log_pari(E,P,Base,order):
-    """
-    Wrapper around pari elllog function to compute discrete
-    logarithms in elliptic curves over finite fields. 
-    Works like P.log(Base).
-    """
-    x = pari.elllog(E,P,Base,order)
-    return ZZ(x)
-
-def weil_pairing_pari(P, Q, D, check=False):
-    """
-    Wrapper around Pari's implementation of the Weil pairing
-    Allows the check of whether P,Q are in E[D] to be optional
-    """
-    if check:
-        nP, nQ = D*P, D*Q
-        if nP.is_zero() or nQ.is_zero():
-            raise ValueError("points must both be n-torsion")
-
-    return pari.ellweilpairing(P.curve(), P, Q, D)
-
-def tate_pairing_pari(P, Q, D):
-    """
-    Wrapper around Pari's implementation of the Tate pairing
-    NOTE: this is not the reduced Tate pairing, so to make this
-    match with SageMath you need
-
-    P.tate_pairing(Q, D, k) == pari.elltatepairing(E, P, Q, D)**((p^k - 1) / D)
-    """
-    E = P.curve()
-    return pari.elltatepairing(E, P, Q, D)
-
-def _precompute_baby_steps(base, step, e):
-    """
-    Helper function to compute the baby steps for 
-    pohlig_hellman_base and windowed_pohlig_hellman.
-    """
-    baby_steps = [base]
-    for _ in range(e):
-        base = base**step
-        baby_steps.append(base)
-    return baby_steps
-
-def pohlig_hellman_base(a, base, e):
-    """
-    Solve the discrete log for a = base^x for 
-    elements base,a of order 2^e using the 
-    Pohlig-Hellman algorithm.
-    """
-    baby_steps = _precompute_baby_steps(base, 2, e)
-
-    dlog = 0
-    exp  = 2**(e-1)
-
-    # Solve the discrete log mod 2, only 2 choices
-    # for each digit!
-    for i in range(e):
-        if a**exp != 1:
-            a /= baby_steps[i]
-            dlog += (2**i)
-
-        if a ==  1:
-            break
-
-        exp //= 2
-
-    return dlog
-
-def windowed_pohlig_hellman(a, base, e, window):
-    """
-    Solve the discrete log for a = base^x for 
-    elements base,a of order 2^e using the 
-    windowed Pohlig-Hellman algorithm following
-    https://ia.cr/2016/963.
-
-    Algorithm runs recursively, computing in windows
-    l^wi for window=[w1, w2, w3, ...]. 
-
-    Runs the base case when window = []
-    """  
-    # Base case when all windows have been used
-    if not window:
-        return pohlig_hellman_base(a, base, e)
-
-    # Collect the next window
-    w, window = window[0], window[1:]
-    step = 2**w
-
-    # When the window is not a divisor of e, we compute
-    # e mod w and solve for both the lower e - e mod w
-    # bits and then the e mod w bits at the end.
-    e_div_w, e_rem = divmod(e, w)
-    e_prime = e - e_rem
-
-    # First force elements to have order e - e_rem
-    a_prime = a**(2**e_rem)
-    base_prime = base**(2**e_rem)
-
-    # Compute base^(2^w*i) for i in (0, ..., e/w-1)
-    baby_steps = _precompute_baby_steps(base_prime, step, e_div_w)
-    
-    # Initialise some pieces
-    dlog = 0
-    if e_prime:
-        exp = 2**(e_prime - w)
-
-        # Work in subgroup of size 2^w
-        s = base_prime**(exp)
-
-        # Windowed Pohlig-Hellman to recover dlog as
-        # alpha = sum l^(i*w) * alpha_i
-        for i in range(e_div_w):
-            # Solve the dlog in 2^w
-            ri = a_prime**exp
-            alpha_i = windowed_pohlig_hellman(ri, s, w, window)
-
-            # Update a value and dlog computation
-            a_prime /= baby_steps[i]**(alpha_i)
-            dlog  += alpha_i * step**i
-
-            if a_prime == 1:
-                break
-
-            exp //= step
-
-    # TODO: 
-    # I don't know if this is a nice way to do
-    # this last step... Works well enough but could
-    # be improved I imagine...
-    exp = 2**e_prime
-    if e_rem:
-        base_last = base**exp
-        a_last = a / base**dlog
-        dlog_last = pohlig_hellman_base(a_last, base_last, e_rem)
-        dlog += exp * dlog_last
-
-    return dlog
-
-def BiDLP(R, P, Q, D, ePQ=None):
-    """
-    Given a basis P,Q of E[D] finds
-    a,b such that R = [a]P + [b]Q.
-
-    Uses the fact that 
-        e([a]P + [b]Q, [c]P + [d]Q) = e(P,Q)^(ad-bc)
-
-    Optional: include the pairing e(P,Q) which can be precomputed
-    which is helpful when running multiple BiDLP problems with P,Q
-    as input. This happens, for example, during compression.
-    """
-    # e(P,Q)
-    if ePQ:
-        pair_PQ = ePQ
-    else:
-        pair_PQ = weil_pairing_pari(P, Q, D)
-
-    # Write R = aP + bQ for unknown a,b
-    # e(R, Q) = e(P, Q)^a
-    pair_a = weil_pairing_pari(R, Q, D)
-
-    # e(R,-P) = e(P, Q)^b
-    pair_b = weil_pairing_pari(R, -P, D)
-
-    # Now solve the dlog in Fq
-    a = discrete_log_pari(pair_a, pair_PQ, D)
-    b = discrete_log_pari(pair_b, pair_PQ, D)
-
-    return a, b
-
-def BiDLP_power_two(R, P, Q, e, window, ePQ=None):
-    r"""
-    Same as the above, but uses optimisations using that
-    D = 2^e.
-
-    First, rather than use the Weil pairing, we can use
-    the Tate pairing which is approx 2x faster. 
-
-    Secondly, rather than solve the discrete log naively,
-    we use an optimised windowed Pohlig-Hellman.
-
-    NOTE: this could be optimised further, following the
-    SIKE key-compression algorithms and for the Tate pairing
-    we could reuse Miller loops to save even more time, but
-    that seemed a little overkill for a SageMath PoC
-
-    Finally, as the Tate pairing produces elements in \mu_n
-    we also have fast inversion from conjugation, but SageMath
-    has slow conjugation, so this doesn't help for now.
-    """
-    p = R.curve().base_ring().characteristic()
-    D = 2**e
-    exp = (p**2 - 1) // D
-
-    # e(P,Q)
-    if ePQ:
-        pair_PQ = ePQ
-    else:
-        pair_PQ = tate_pairing_pari(P, Q, D)**exp
-
-    # Write R = aP + bQ for unknown a,b
-    # e(R, Q) = e(P, Q)^a
-    pair_a = tate_pairing_pari(Q, -R, D)**exp
-
-    # e(R,-P) = e(P, Q)^b
-    pair_b = tate_pairing_pari(P, R, D)**exp
-
-    # Now solve the dlog in Fq
-    a = windowed_pohlig_hellman(pair_a, pair_PQ, e, window)
-    b = windowed_pohlig_hellman(pair_b, pair_PQ, e, window)
-
-    return a, b
-
-def DLP_power_two(R, P, Q, e, window, ePQ=None, first=True):
-    r"""
-    This is the same as BiDLP but it only returns either a or b 
-    depending on whether first is true or false.
-    This is used in compression, where we only send 3 of the 4 
-    scalars from BiDLP
-    """
-    p = R.curve().base_ring().characteristic()
-    D = 2**e
-    exp = (p**2 - 1) // D
-
-    # e(P,Q)
-    if ePQ:
-        pair_PQ = ePQ
-    else:
-        pair_PQ = tate_pairing_pari(P, Q, D)**exp
-
-    if first:
-        pair_a = tate_pairing_pari(Q, -R, D)**exp
-        x = windowed_pohlig_hellman(pair_a, pair_PQ, e, window)
-
-    else:
-        pair_b = tate_pairing_pari(P, R, D)**exp
-        x = windowed_pohlig_hellman(pair_b, pair_PQ, e, window)
-
-    return x
-
diff --git a/theta_lib/utilities/fast_sqrt.py b/theta_lib/utilities/fast_sqrt.py
deleted file mode 100644
index 0609fe5..0000000
--- a/theta_lib/utilities/fast_sqrt.py
+++ /dev/null
@@ -1,55 +0,0 @@
-
-# ============================================ #
-#     Fast square root and quadratic roots     #
-# ============================================ #
-
-"""
-Most of this code has been taken from:
-https://github.com/FESTA-PKE/FESTA-SageMath
-
-Copyright (c) 2023 Andrea Basso, Luciano Maino and Giacomo Pope.
-
-Functions with another Copyright mention are not from the above authors.
-"""
-
-def sqrt_Fp2(a):
-    """
-    Efficiently computes the sqrt
-    of an element in Fp2 using that
-    we always have a prime p such that
-    p ≡ 3 mod 4.
-    """
-    Fp2 = a.parent()
-    p = Fp2.characteristic()
-    i = Fp2.gen() # i = √-1
-
-    a1 = a ** ((p - 3) // 4)
-    x0 = a1 * a
-    alpha = a1 * x0
-
-    if alpha == -1:
-        x = i * x0
-    else:
-        b = (1 + alpha) ** ((p - 1) // 2)
-        x = b * x0
-
-    return x
-
-def n_sqrt(a, n):
-    for _ in range(n):
-        a = sqrt_Fp2(a)
-    return a
-
-def sqrt_Fp(a):
-    """
-    Efficiently computes the sqrt
-    of an element in Fp using that
-    we always have a prime p such that
-    p ≡ 3 mod 4.
-
-    Copyright (c) Pierrick Dartois 2025.
-    """
-    Fp = a.parent()
-    p = Fp.characteristic()
-
-    return a**((p+1)//4)
diff --git a/theta_lib/utilities/order.py b/theta_lib/utilities/order.py
deleted file mode 100644
index 223f568..0000000
--- a/theta_lib/utilities/order.py
+++ /dev/null
@@ -1,117 +0,0 @@
-# ================================================== #
-#  Code to check whether a group element has order D #
-# ================================================== #
-
-"""
-This code has been taken from:
-https://github.com/FESTA-PKE/FESTA-SageMath
-
-Copyright (c) 2023 Andrea Basso, Luciano Maino and Giacomo Pope.
-"""
-
-from sage.all import(
-    ZZ,
-    prod,
-    cached_function
-)
-
-def batch_cofactor_mul_generic(G_list, pis, group_action, lower, upper):
-    """
-    Input:  A list of elements `G_list`, such that
-                G is the first entry and the rest is empty
-                in the sublist G_list[lower:upper]
-            A list `pis` of primes p such that
-                their product is D
-            The `group_action` of the group
-            Indices lower and upper
-    Output: None
-
-    NOTE: G_list is created in place
-    """
-
-    # check that indices are valid
-    if lower > upper:
-        raise ValueError(f"Wrong input to cofactor_multiples()")
-
-    # last recursion step does not need any further splitting
-    if upper - lower == 1:
-        return
-
-    # Split list in two parts,
-    # multiply to get new start points for the two sublists,
-    # and call the function recursively for both sublists.
-    mid = lower + (upper - lower + 1) // 2
-    cl, cu = 1, 1
-    for i in range(lower, mid):
-        cu = cu * pis[i]
-    for i in range(mid, upper):
-        cl = cl * pis[i]
-    # cl = prod(pis[lower:mid])
-    # cu = prod(pis[mid:upper])
-
-    G_list[mid] = group_action(G_list[lower], cu)
-    G_list[lower] = group_action(G_list[lower], cl)
-
-    batch_cofactor_mul_generic(G_list, pis, group_action, lower, mid)
-    batch_cofactor_mul_generic(G_list, pis, group_action, mid, upper)
-
-
-@cached_function
-def has_order_constants(D):
-    """
-    Helper function, finds constants to
-    help with has_order_D
-    """
-    D = ZZ(D)
-    pis = [p for p, _ in D.factor()]
-    D_radical = prod(pis)
-    Dtop = D // D_radical
-    return Dtop, pis
-
-
-def has_order_D(G, D, multiplicative=False):
-    """
-    Given an element G in a group, checks if the
-    element has order exactly D. This is much faster
-    than determining its order, and is enough for 
-    many checks we need when computing the torsion
-    basis.
-
-    We allow both additive and multiplicative groups
-    which means we can use this when computing the order
-    of points and elements in Fp^k when checking the 
-    multiplicative order of the Weil pairing output
-    """
-    # For the case when we work with elements of Fp^k
-    if multiplicative:
-        group_action = lambda a, k: a**k
-        is_identity = lambda a: a == 1
-        identity = 1
-    # For the case when we work with elements of E / Fp^k
-    else:
-        group_action = lambda a, k: k * a
-        is_identity = lambda a: a.is_zero()
-        identity = G.curve()(0)
-
-    if is_identity(G):
-        return False
-
-    D_top, pis = has_order_constants(D)
-
-    # If G is the identity after clearing the top
-    # factors, we can abort early
-    Gtop = group_action(G, D_top)
-    if is_identity(Gtop):
-        return False
-
-    G_list = [identity for _ in range(len(pis))]
-    G_list[0] = Gtop
-
-    # Lastly we have to determine whether removing any prime 
-    # factors of the order gives the identity of the group
-    if len(pis) > 1:
-        batch_cofactor_mul_generic(G_list, pis, group_action, 0, len(pis))
-        if not all([not is_identity(G) for G in G_list]):
-            return False
-
-    return True
\ No newline at end of file
diff --git a/theta_lib/utilities/strategy.py b/theta_lib/utilities/strategy.py
deleted file mode 100644
index 550ef09..0000000
--- a/theta_lib/utilities/strategy.py
+++ /dev/null
@@ -1,229 +0,0 @@
-# ============================================================================ #
-#     Compute optimised strategy for 2-isogeny chains (in dimensions 2 and 4)  #
-# ============================================================================ #
-
-"""
-The function optimised_strategy has been taken from:
-https://github.com/FESTA-PKE/FESTA-SageMath
-
-Copyright (c) 2023 Andrea Basso, Luciano Maino and Giacomo Pope.
-
-Other functions are original work.
-"""
-
-from sage.all import log, cached_function
-import logging
-logger = logging.getLogger(__name__)
-#logger.setLevel("DEBUG")
-
-@cached_function
-def optimised_strategy(n, mul_c=1):
-    """
-    Algorithm 60: https://sike.org/files/SIDH-spec.pdf
-    Shown to be appropriate for (l,l)-chains in 
-    https://ia.cr/2023/508
-    
-    Note: the costs we consider are:
-       eval_c: the cost of one isogeny evaluation
-       mul_c:  the cost of one element doubling
-    """
-
-    eval_c = 1.000
-    mul_c  = mul_c
-
-    S = {1:[]}
-    C = {1:0 }
-    for i in range(2, n+1):
-        b, cost = min(((b, C[i-b] + C[b] + b*mul_c + (i-b)*eval_c) for b in range(1,i)), key=lambda t: t[1])
-        S[i] = [b] + S[i-b] + S[b]
-        C[i] = cost
-
-    return S[n]
-
-@cached_function
-def optimised_strategy_with_first_eval(n,mul_c=1,first_eval_c=1):
-    r"""
-    Adapted from optimised_strategy when the fist isogeny evaluation is more costly.
-    This is well suited to gluing comptations. Computes optimal strategies with constraint
-    at the beginning. This takes into account the fact that doublings on the codomain of 
-    the first isogeny are impossible (because of zero dual theta constants).
-
-    CAUTION: When splittings are involved, do not use this function. Use 
-    optimised_strategy_with_first_eval_and_splitting instead.
-
-    INPUT:
-    - n: number of leaves of the strategy (length of the isogeny).
-    - mul_c: relative cost of one doubling compared to one generic 2-isogeny evaluation.
-    - first_eval_c: relative cost of an evaluation of the first 2-isogeny (gluing) 
-    compared to one generic 2-isogeny evaluation.
-
-    OUTPUT:
-    - S_left[n]: an optimal strategy of depth n with constraint at the beginning
-    represented as a sequence [s_0,...,s_{t-2}], where there is an index i for every 
-    internal node of the strategy, where indices are ordered depth-first left-first 
-    (as the way we move on the strategy) and s_i is the number of leaves to the right 
-    of internal node i (see https://sike.org/files/SIDH-spec.pdf, pp. 16-17).
-    """
-
-    S_left, _, _, _ = optimised_strategies_with_first_eval(n,mul_c,first_eval_c)
-
-    return S_left[n]
-
-@cached_function
-def optimised_strategies_with_first_eval(n,mul_c=1,first_eval_c=1):
-    r"""
-    Adapted from optimised_strategy when the fist isogeny evaluation is more costly.
-    This is well suited to gluing comptations. Computes optimal strategies with constraint
-    at the beginning. This takes into account the fact that doublings on the codomain of 
-    the first isogeny are impossible (because of zero dual theta constants).
-
-    CAUTION: When splittings are involved, do not use this function. Use 
-    optimised_strategy_with_first_eval_and_splitting instead.
-
-    INPUT:
-    - n: number of leaves of the strategy (length of the isogeny).
-    - mul_c: relative cost of one doubling compared to one generic 2-isogeny evaluation.
-    - first_eval_c: relative cost of an evaluation of the first 2-isogeny (gluing) 
-    compared to one generic 2-isogeny evaluation.
-
-    OUTPUT:
-    - S_left: Optimal strategies "on the left"/with constraint at the beginning i.e. meeting the 
-    first left edge that do not contain any left edge on the line y=sqrt(3)*(x-1).
-    - S_right: Optimal strategies "on the right" i.e. not meeting the fisrt left edge (no constraint).
-    """
-
-    # print(f"Strategy eval: n={n}, mul_c={mul_c}, first_eval_c={first_eval_c}")
-
-    eval_c = 1.000
-    first_eval_c = first_eval_c
-    mul_c  = mul_c
-
-    S_left = {1:[], 2:[1]} # Optimal strategies "on the left" i.e. meeting the first left edge 
-    S_right = {1:[]} # Optimal strategies "on the right" i.e. not meeting the first left edge
-    C_left = {1:0, 2:mul_c+first_eval_c } # Cost of strategies on the left
-    C_right = {1:0 } # Cost of strategies on the right
-    for i in range(2, n+1):
-        # Optimisation on the right
-        b, cost = min(((b, C_right[i-b] + C_right[b] + b*mul_c + (i-b)*eval_c) for b in range(1,i)), key=lambda t: t[1])
-        S_right[i] = [b] + S_right[i-b] + S_right[b]
-        C_right[i] = cost
-
-    for i in range(3,n+1):
-        # Optimisation on the left (b<i-1 to avoid doublings on the codomain of the first isogeny)
-        b, cost = min(((b, C_left[i-b] + C_right[b] + b*mul_c + (i-1-b)*eval_c+first_eval_c) for b in range(1,i-1)), key=lambda t: t[1])
-        S_left[i] = [b] + S_left[i-b] + S_right[b]
-        C_left[i] = cost
-
-    return S_left, S_right, C_left, C_right
-
-@cached_function
-def optimised_strategy_with_first_eval_and_splitting(n,m,mul_c=1,first_eval_c=1):
-    eval_c = 1.000
-    first_eval_c = first_eval_c
-    mul_c  = mul_c
-
-    S_left, S_middle, C_left, C_middle = optimised_strategies_with_first_eval(n-m,mul_c,first_eval_c)
-
-    ## Optimization of the right part (with constraint at the end only)
-
-    # Dictionnary of dictionnaries of translated strategies "on the right".
-    # trans_S_right[d][i] is an optimal strategy of depth i 
-    # without left edge on the line y=sqrt(3)*(x-(i-1-d))
-    trans_S_right={}
-    trans_C_right={}
-
-    for d in range(1,m+1):
-        trans_S_right[d]={1:[]}
-        trans_C_right[d]={1:0}
-        if d==1:
-            for i in range(3,n-m+d):
-                b, cost = min(((b, C_middle[i-b] + trans_C_right[1][b] + b*mul_c + (i-b)*eval_c) for b in [1]+list(range(3,i))), key=lambda t: t[1])
-                trans_S_right[1][i] = [b] + S_middle[i-b] + trans_S_right[1][b]
-                trans_C_right[1][i] = cost
-        else:
-            for i in range(2,n-m+d):
-                if i!=d+1:
-                    b = 1
-                    cost = trans_C_right[d-b][i-b] + C_middle[b] + b*mul_c + (i-b)*eval_c
-                    for k in range(2,min(i,d)):
-                        cost_k = trans_C_right[d-k][i-k] + C_middle[k] + k*mul_c + (i-k)*eval_c
-                        if cost_k<cost:
-                            b = k
-                            cost = cost_k
-                    # k=d
-                    if i>d:
-                        cost_k = C_middle[i-d] + C_middle[d] + d*mul_c + (i-d)*eval_c
-                        if cost_k<cost:
-                            b = d
-                            cost = cost_k
-                    for k in range(d+2,i):
-                        #print(d,i,k)
-                        cost_k = C_middle[i-k] + trans_C_right[d][k] + k*mul_c + (i-k)*eval_c
-                        if cost_k<cost:
-                            b = k
-                            cost = cost_k
-                    if b<d:
-                        trans_S_right[d][i] = [b] + trans_S_right[d-b][i-b] + S_middle[b]
-                        trans_C_right[d][i] = cost
-                    else:
-                        trans_S_right[d][i] = [b] + S_middle[i-b] + trans_S_right[d][b]
-                        trans_C_right[d][i] = cost
-
-    ## Optimization on the left (last part) taking into account the constraints at the beginning and at the end
-    for i in range(n-m+1,n+1):
-        d = i-(n-m)
-        b = 1
-        cost = C_left[i-b] + trans_C_right[d][b] + b*mul_c + (i-1-b)*eval_c + first_eval_c
-        for k in range(2,i):
-            if k!=d+1 and k!=n-1:
-                cost_k = C_left[i-k] + trans_C_right[d][k] + k*mul_c + (i-1-k)*eval_c + first_eval_c
-                if cost_k<cost:
-                    b = k
-                    cost = cost_k
-
-        S_left[i] = [b] + S_left[i-b] + trans_S_right[d][b]
-        C_left[i] = cost
-
-    return S_left[n]
-
-@cached_function
-def precompute_strategy_with_first_eval(e,m,M=1,S=0.8,I=100):
-    r"""
-    INPUT:
-    - e: isogeny chain length.
-    - m: length of the chain in dimension 2 before gluing in dimension 4.
-    - M: multiplication cost.
-    - S: squaring cost.
-    - I: inversion cost.
-
-    OUTPUT: Optimal strategy to compute an isogeny chain without splitting of
-    length e with m steps in dimension 2 before gluing in dimension 4.
-    """
-    n = e - m
-    eval_c = 4*(16*M+16*S)
-    mul_c = 4*(32*M+32*S)+(48*M+I)/log(n*1.0)
-    first_eval_c = 4*(2*I+(244+6*m)*M+(56+8*m)*S)
-
-    return optimised_strategy_with_first_eval(n, mul_c = mul_c/eval_c, first_eval_c = first_eval_c/eval_c)
-
-@cached_function
-def precompute_strategy_with_first_eval_and_splitting(e,m,M=1,S=0.8,I=100):
-    r"""
-    INPUT:
-    - e: isogeny chain length.
-    - m: length of the chain in dimension 2 before gluing in dimension 4.
-    - M: multiplication cost.
-    - S: squaring cost.
-    - I: inversion cost.
-
-    OUTPUT: Optimal strategy to compute an isogeny chain of length e 
-    with m steps in dimension 2 before gluing in dimension 4 and
-    with splitting m steps before the end.
-    """
-    logger.debug(f"Strategy eval split: e={e}, m={m}")
-    n = e - m
-    eval_c = 4*(16*M+16*S)
-    mul_c = 4*(32*M+32*S)+(48*M+I)/log(n*1.0)
-    first_eval_c = 4*(2*I+(244+6*m)*M+(56+8*m)*S)
-
-    return optimised_strategy_with_first_eval_and_splitting(n, m, mul_c = mul_c/eval_c, first_eval_c = first_eval_c/eval_c)
diff --git a/theta_lib/utilities/supersingular.py b/theta_lib/utilities/supersingular.py
deleted file mode 100644
index e760c1b..0000000
--- a/theta_lib/utilities/supersingular.py
+++ /dev/null
@@ -1,413 +0,0 @@
-"""
-Helper functions for the supersingular elliptic curve computations in FESTA
-"""
-
-# =========================================== #
-# Compute points of order D and Torsion Bases #
-# =========================================== #
-
-"""
-Most of this code has been taken from:
-https://github.com/FESTA-PKE/FESTA-SageMath
-
-Copyright (c) 2023 Andrea Basso, Luciano Maino and Giacomo Pope.
-
-Functions with another Copyright mention are not from the above authors.
-"""
-
-# Sage Imports
-from sage.all import ZZ, GF, Integers
-
-# Local imports
-from .order import has_order_D
-from .discrete_log import weil_pairing_pari, ell_discrete_log_pari, discrete_log_pari
-from .fast_sqrt import sqrt_Fp2, sqrt_Fp
-
-
-
-
-def random_point(E):
-    """
-    Returns a random point on the elliptic curve E
-    assumed to be in Montgomery form with a base 
-    field which characteristic p = 3 mod 4
-    """
-    A = E.a_invariants()[1]
-    if E.a_invariants() != (0,A,0,1,0):
-        raise ValueError("Function `generate_point` assumes the curve E is in the Montgomery model")
-
-    # Try 10000 times then give up, just protection
-    # for infinite loops
-    F = E.base_ring()
-    for _ in range(10000):
-        x = F.random_element()
-        y2 = x*(x**2 + A*x + 1)
-        if y2.is_square():
-            y = sqrt_Fp2(y2)
-            return E(x, y)
-
-    raise ValueError("Generated 10000 points, something is probably going wrong somewhere.")
-
-def generate_point(E, x_start=0):
-    """
-    Generate points on a curve E with x-coordinate 
-    i + x for x in Fp and i is the generator of Fp^2
-    such that i^2 = -1.
-    """
-    F = E.base_ring()
-    one = F.one()
-
-    if x_start:
-        x = x_start + one
-    else:
-        x = F.gen() + one
-
-    A = E.a_invariants()[1]
-    if E.a_invariants() != (0,A,0,1,0):
-        raise ValueError("Function `generate_point` assumes the curve E is in the Montgomery model")
-
-    # Try 10000 times then give up, just protection
-    # for infinite loops
-    for _ in range(10000):
-        y2 = x*(x**2 + A*x + 1)
-        if y2.is_square():
-            y = sqrt_Fp2(y2)
-            yield E(x, y)
-        x += one
-
-    raise ValueError("Generated 10000 points, something is probably going wrong somewhere.")
-
-def generate_Fp_point(E):
-    """
-    Returns a random Fp-rational point on the 
-    elliptic curve E assumed to be in Montgomery 
-    form with a base field which characteristic 
-    p = 3 mod 4 and coefficients defined over Fp.
-
-    Copyright (c) Pierrick Dartois 2025.
-    """
-
-    A = E.a_invariants()[1]
-    if E.a_invariants() != (0,A,0,1,0):
-        raise ValueError("Function `generate_Fp_point` assumes the curve E is in the Montgomery model")
-    if A[1] != 0:
-        raise ValueError("The curve E should be Fp-rational")
-
-    p=E.base_field().characteristic()
-    Fp=GF(p,proof=False)
-    one=Fp(1)
-    x=one
-
-    # Try 10000 times then give up, just protection
-    # for infinite loops
-    for _ in range(10000):
-        y2 = x*(x**2 + A*x + 1)
-        if y2.is_square():
-            y = sqrt_Fp(y2)
-            yield E(x, y)
-        x += one
-
-    raise ValueError("Generated 10000 points, something is probably going wrong somewhere.")
-
-def generate_point_order_D(E, D, x_start=0):
-    """
-    Input:  An elliptic curve E / Fp2
-            An integer D dividing (p +1)
-    Output: A point P of order D.
-    """
-    p = E.base().characteristic()
-    n = (p + 1) // D
-
-    Ps = generate_point(E, x_start=x_start)
-    for G in Ps:
-        P = n * G
-
-        # Case when we randomly picked
-        # a point in the n-torsion
-        if P.is_zero():
-            continue
-
-        # Check that P has order exactly D
-        if has_order_D(P, D):
-            P._order = ZZ(D)
-            yield P
-
-    raise ValueError(f"Never found a point P of order D.")
-
-def generate_Fp_point_order_D(E, D):
-    """
-    Input:  An elliptic curve E / Fp2
-            An integer D dividing (p +1)
-    Output: A point P of order D.
-
-    Copyright (c) Pierrick Dartois 2025.
-    """
-    p = E.base().characteristic()
-    n = (p + 1) // D
-    if D%2==0:
-        n=n//2
-
-    Ps = generate_Fp_point(E)
-    for G in Ps:
-        P = n * G
-
-        # Case when we randomly picked
-        # a point in the n-torsion
-        if P.is_zero():
-            continue
-
-        # Check that P has order exactly D
-        if has_order_D(P, D):
-            P._order = ZZ(D)
-            yield P
-
-    raise ValueError(f"Never found a point P of order D.")
-
-def compute_point_order_D(E, D, x_start=0):
-    """
-    Wrapper function around a generator which returns the first
-    point of order D
-    """
-    return generate_point_order_D(E, D, x_start=x_start).__next__()
-
-def compute_Fp_point_order_D(E, D):
-    """
-    Wrapper function around a generator which returns the first
-    point of order D
-
-    Copyright (c) Pierrick Dartois 2025.
-    """
-    return generate_Fp_point_order_D(E, D).__next__()
-
-def compute_linearly_independent_point_with_pairing(E, P, D, x_start=0):
-    """
-    Input:  An elliptic curve E / Fp2
-            A point P ∈ E[D]
-            An integer D dividing (p +1)
-    Output: A point Q such that E[D] = <P, Q>
-            The Weil pairing e(P,Q)
-    """
-    Qs = generate_point_order_D(E, D, x_start=x_start)
-    for Q in Qs:
-        # Make sure the point is linearly independent
-        pair = weil_pairing_pari(P, Q, D)
-        if has_order_D(pair, D, multiplicative=True):
-            Q._order = ZZ(D)
-            return Q, pair
-    raise ValueError("Never found a point Q linearly independent to P")
-
-def compute_linearly_independent_point(E, P, D, x_start=0):
-    """
-    Wrapper function around `compute_linearly_independent_point_with_pairing`
-    which only returns a linearly independent point
-    """
-    Q, _ = compute_linearly_independent_point_with_pairing(E, P, D, x_start=x_start)
-    return Q
-
-def torsion_basis_with_pairing(E, D):
-    """
-    Generate basis of E(Fp^2)[D] of supersingular curve
-
-    While computing E[D] = <P, Q> we naturally compute the 
-    Weil pairing e(P,Q), which we also return as in some cases
-    the Weil pairing is then used when solving the BiDLP
-    """
-    p = E.base().characteristic()
-
-    # Ensure D divides the curve's order
-    if (p + 1) % D != 0:
-        print(f"{ZZ(D).factor() = }")
-        print(f"{ZZ(p+1).factor() = }")
-        raise ValueError(f"D must divide the point's order")
-
-    P = compute_point_order_D(E, D)
-    Q, ePQ = compute_linearly_independent_point_with_pairing(E, P, D, x_start=P[0])
-
-    return P, Q, ePQ
-
-def torsion_basis(E, D):
-    """
-    Wrapper function around torsion_basis_with_pairing which only
-    returns the torsion basis <P,Q> = E[D]
-    """
-    P, Q, _ = torsion_basis_with_pairing(E, D)
-    return P, Q
-
-def torsion_basis_2f_frobenius(E,f):
-    """
-    Returns a basis (P, Q) of E[2^f] with pi_p(P)=P and pi_p(Q)=-Q,
-    pi_p being the p-th Frobenius endomorphism of E.
-
-    Copyright (c) Pierrick Dartois 2025.
-    """
-
-    P = compute_Fp_point_order_D(E, 2**f)
-
-    i = E.base_field().gen()
-    p = E.base_field().characteristic() 
-
-    R = compute_linearly_independent_point(E, P, 2**f)
-    piR = E(R[0]**p,R[1]**p)
-    piR_p_R = piR + R
-
-    two_lamb = ell_discrete_log_pari(E,piR_p_R,P,2**f)
-    lamb = two_lamb//2
-    Q = R-lamb*P
-
-    piQ = E(Q[0]**p,Q[1]**p)
-
-    assert piQ == -Q
-
-    return P, Q 
-
-def torsion_basis_to_Fp_rational_point(E,P,Q,D):
-    """
-    Returns an Fp-rational point R=[lamb]P+[mu]Q
-    of D-torsion given a D-torsion basis (P,Q) of E.
-
-    Copyright (c) Pierrick Dartois 2025.
-    """
-
-    p=E.base_field().characteristic()
-    piP = E(P[0]**p,P[1]**p)
-    piQ = E(Q[0]**p,Q[1]**p)
-
-    e_PQ = weil_pairing_pari(P,Q,D)
-    e_PpiP = weil_pairing_pari(P, piP, D)
-    e_QpiP = weil_pairing_pari(Q, piP, D)
-    e_PpiQ = weil_pairing_pari(P, piQ, D)
-    e_QpiQ = weil_pairing_pari(Q, piQ, D)
-
-    # pi(P)=[a]P+[b]Q, pi(Q)=[c]P+[d]Q
-    a = -discrete_log_pari(e_QpiP,e_PQ,D)
-    b = discrete_log_pari(e_PpiP,e_PQ,D)
-    c = -discrete_log_pari(e_QpiQ,e_PQ,D)
-    d = discrete_log_pari(e_PpiQ,e_PQ,D)
-
-    ZD = Integers(D)
-
-    if a%D==1:
-        mu = ZD(b/d)
-        R=P+mu*Q
-    elif c==0:
-        # c=0 and a!=1 => d=1 so pi(Q)=Q
-        R=Q
-    else:
-        # c!=0 and a!=1
-        mu = ZD((1-a)/c)
-        R=P+mu*Q
-
-    piR = E(R[0]**p,R[1]**p)
-
-    assert piR==R
-
-    return R
-
-# =========================================== #
-#   Entangled torsion basis for fast E[2^k]   #
-# =========================================== #
-
-def precompute_elligator_tables(F):
-    """
-    Precomputes tables of quadratic non-residue or 
-    quadratic residue in Fp2. Used to compute entangled
-    torsion bases following https://ia.cr/2017/1143
-    """
-    u = 2*F.gen()
-
-    T1 = dict()
-    T2 = dict()
-    # TODO: estimate how large r should be
-    for r in range(1, 30):
-        v = 1 / (1 + u*r**2)
-        if v.is_square():
-            T2[r] = v
-        else:
-            T1[r] = v
-    return T1, T2
-
-def entangled_torsion_basis(E, elligator_tables, cofactor):
-    """
-    Optimised algorithm following https://ia.cr/2017/1143
-    which modifies the elligator method of hashing to points
-    to find points P,Q of order k*2^b. Clearing the cofactor
-    gives the torsion basis without checking the order or
-    computing a Weil pairing. 
-
-    To do this, we need tables TQNR, TQR of pairs of values
-    (r, v) where r is an integer and v = 1/(1 + ur^2) where
-    v is either a quadratic non-residue or quadratic residue 
-    in Fp2 and u = 2i = 2*sqrt(-1).
-    """
-    F = E.base_ring()
-    p = F.characteristic()
-    p_sqrt = (p+1)//4
-
-    i = F.gen()
-    u =  2 * i
-    u0 = 1 + i
-
-    TQNR, TQR = elligator_tables
-    
-    # Pick the look up table depending on whether
-    # A = a + ib is a QR or NQR
-    A = E.a_invariants()[1]
-    if (0,A,0,1,0) != E.a_invariants():
-        raise ValueError("The elliptic curve E must be in Montgomery form")
-    if A.is_square():
-        T = TQNR
-    else:
-        T = TQR
-
-    # Look through the table to find point with 
-    # rational (x,y)
-    y = None
-    for r, v in T.items():
-        x = -A * v
-
-        t = x * (x**2 + A*x + 1)
-        
-        # Break when we find rational y: t = y^2 
-        c, d = t.list()
-        z = c**2 + d**2
-        s = z**p_sqrt
-        if s**2 == z:
-            y = sqrt_Fp2(t)
-            break
-    
-    if y is None:
-        raise ValueError("Never found a y-coordinate, increase the lookup table size")
-
-    z = (c + s) // 2
-    alpha = z**p_sqrt
-    beta  = d / (2*alpha)
-
-    if alpha**2 == z:
-        y =  F([alpha, beta])
-    else:
-        y = -F([beta, alpha])
-
-    S1 = E([x, y])
-    S2 = E([u*r**2*x, u0*r*y])
-
-    return cofactor*S1, cofactor*S2
-
-# =============================================== #
-#  Ensure Basis <P,Q> of E[2^k] has (0,0) under Q #
-# =============================================== #
-
-def fix_torsion_basis_renes(P, Q, k):
-    """
-    Set the torsion basis P,Q such that
-    2^(k-1)Q = (0,0) to ensure that (0,0)
-    is never a kernel of a two isogeny
-    """
-    cofactor = 2**(k-1)
-
-    R = cofactor*P
-    if R[0] == 0:
-        return Q, P
-    R = cofactor*Q
-    if R[0] == 0:
-        return P, Q
-    return P, P + Q
diff --git a/time_pegasis.py b/time_pegasis.py
index db0ab62..c81c745 100644
--- a/time_pegasis.py
+++ b/time_pegasis.py
@@ -5,7 +5,8 @@ import itertools
 from time import time
 
 from pegasis import PEGASIS
-from theta_lib.utilities.strategy import precompute_strategy_with_first_eval_and_splitting as precompute
+#from theta_lib.utilities.strategy import precompute_strategy_with_first_eval_and_splitting as precompute
+from Theta_dim4.Theta_dim4_sage.pkg.utilities.strategy import precompute_strategy_with_first_eval as precompute
 
 
 def precompute_strategies(emax):