Initial Commit

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

6 files changed, 1329 insertions, 0 deletions

diff --git a/theta_lib/isogenies/Kani_clapoti.py b/theta_lib/isogenies/Kani_clapoti.py
new file mode 100644
index 0000000..66030e2
--- /dev/null
+++ b/theta_lib/isogenies/Kani_clapoti.py
@@ -0,0 +1,258 @@
+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
new file mode 100644
index 0000000..282219c
--- /dev/null
+++ b/theta_lib/isogenies/Kani_gluing_isogeny_chain_dim4.py
@@ -0,0 +1,567 @@
+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
new file mode 100644
index 0000000..cd738c9
--- /dev/null
+++ b/theta_lib/isogenies/gluing_isogeny_dim4.py
@@ -0,0 +1,200 @@
+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
new file mode 100644
index 0000000..8c0b78e
--- /dev/null
+++ b/theta_lib/isogenies/isogeny_chain_dim4.py
@@ -0,0 +1,114 @@
+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
new file mode 100644
index 0000000..2f483bf
--- /dev/null
+++ b/theta_lib/isogenies/isogeny_dim4.py
@@ -0,0 +1,162 @@
+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
new file mode 100644
index 0000000..a6e3da3
--- /dev/null
+++ b/theta_lib/isogenies/tree.py
@@ -0,0 +1,28 @@
+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
+