diff options
Diffstat (limited to 'theta_lib/theta_structures')
| -rw-r--r-- | theta_lib/theta_structures/Theta_dim1.py | 98 | ||||
| -rw-r--r-- | theta_lib/theta_structures/Theta_dim2.py | 440 | ||||
| -rw-r--r-- | theta_lib/theta_structures/Theta_dim4.py | 351 | ||||
| -rw-r--r-- | theta_lib/theta_structures/Tuple_point.py | 106 | ||||
| -rw-r--r-- | theta_lib/theta_structures/montgomery_theta.py | 68 | ||||
| -rw-r--r-- | theta_lib/theta_structures/theta_helpers_dim2.py | 32 | ||||
| -rw-r--r-- | theta_lib/theta_structures/theta_helpers_dim4.py | 259 |
7 files changed, 0 insertions, 1354 deletions
diff --git a/theta_lib/theta_structures/Theta_dim1.py b/theta_lib/theta_structures/Theta_dim1.py deleted file mode 100644 index e8e94dd..0000000 --- a/theta_lib/theta_structures/Theta_dim1.py +++ /dev/null @@ -1,98 +0,0 @@ -from sage.all import * -from ..utilities.supersingular import compute_linearly_independent_point -from .montgomery_theta import ( - montgomery_point_to_theta_point, - theta_point_to_montgomery_point, - torsion_to_theta_null_point, -) - - -class ThetaStructureDim1: - def __init__(self,E,P=None,Q=None): - self.E=E - - a_inv=E.a_invariants() - - A =a_inv[1] - if a_inv != (0,A,0,1,0): - raise ValueError("The elliptic curve E is not in the Montgomery model.") - - if Q==None: - y2=A-2 - y=y2.sqrt() - Q=E([-1,y,1]) - P=compute_linearly_independent_point(E,Q,4) - else: - if Q[0]!=-1: - raise ValueError("You should enter a canonical 4-torsion basis. Q[0] should be -1.") - - self.P=P - self.Q=Q - - self._base_ring=E.base_ring() - - self._point=ThetaPointDim1 - self._null_point=self._point(self,torsion_to_theta_null_point(P)) - - def null_point(self): - """ - """ - return self._null_point - - def base_ring(self): - """ - """ - return self._base_ring - - def zero(self): - """ - """ - return self.null_point() - - def elliptic_curve(self): - return self.E - - def torsion_basis(self): - return (self.P,self.Q) - - def __call__(self,coords): - r""" - Input: either a tuple or list of 2 coordinates or an elliptic curve point. - - Output: the corresponding theta point for the self theta structure. - """ - if isinstance(coords,tuple): - return self._point(self,coords) - elif isinstance(coords,list): - return self._point(self,coords) - else: - return self._point(self,montgomery_point_to_theta_point(self.null_point().coords(),coords)) - - def __repr__(self): - return f"Theta structure on {self.elliptic_curve()} with null point: {self.null_point()} induced by the 4-torsion basis {self.torsion_basis()}" - - def to_montgomery_point(self,P): - return theta_point_to_montgomery_point(self.E,self.null_point().coords(),P.coords()) - - -class ThetaPointDim1: - def __init__(self, parent, coords): - """ - """ - if not isinstance(parent, ThetaStructureDim1): - raise ValueError("Entry parent should be a ThetaStructureDim1 object.") - - self._parent = parent - self._coords = tuple(coords) - - - def coords(self): - return self._coords - - def __repr__(self): - return f"Theta point with coordinates: {self.coords()}" - - - - - diff --git a/theta_lib/theta_structures/Theta_dim2.py b/theta_lib/theta_structures/Theta_dim2.py deleted file mode 100644 index 0d9955a..0000000 --- a/theta_lib/theta_structures/Theta_dim2.py +++ /dev/null @@ -1,440 +0,0 @@ -# Sage Imports -from sage.all import Integer -from sage.structure.element import get_coercion_model, RingElement - -cm = get_coercion_model() - -from .theta_helpers_dim2 import batch_inversion, product_theta_point -from .Theta_dim1 import ThetaStructureDim1 -from .Tuple_point import TuplePoint -from ..basis_change.base_change_dim2 import apply_base_change_theta_dim2 - -# =========================================== # -# Class for Theta Structure (level-2) # -# =========================================== # - - -class ThetaStructureDim2: - """ - Class for the Theta Structure in dimension 2, defined by its theta null point. This type - represents the generic domain/codomain of the (2,2)-isogeny in the theta model. - """ - - def __init__(self, null_point, null_point_dual=None): - if not len(null_point) == 4: - raise ValueError - - self._base_ring = cm.common_parent(*(c.parent() for c in null_point)) - self._point = ThetaPointDim2 - self._precomputation = None - - self._null_point = self._point(self, null_point) - self._null_point_dual = null_point_dual - - def null_point(self): - """ - Return the null point of the given theta structure - """ - return self._null_point - - def null_point_dual(self): - if self._null_point_dual==None: - self._null_point_dual = self._point.to_hadamard(*self.coords()) - return self._null_point_dual - - def base_ring(self): - """ - Return the base ring of the common parent of the coordinates of the null point - """ - return self._base_ring - - def zero(self): - """ - The additive identity is the theta null point - """ - return self.null_point() - - def zero_dual(self): - return self.null_point_dual() - - def __repr__(self): - return f"Theta structure over {self.base_ring()} with null point: {self.null_point()}" - - def coords(self): - """ - Return the coordinates of the theta null point of the theta structure - """ - return self.null_point().coords() - - def hadamard(self): - """ - Compute the Hadamard transformation of the theta structure - """ - return ThetaStructureDim2(self.null_point_dual(),null_point_dual=self.coords()) - - def squared_theta(self): - """ - Square the coefficients and then compute the Hadamard transformation of - the theta null point of the theta structure - """ - return self.null_point().squared_theta() - - def _arithmetic_precomputation(self): - """ - Precompute 6 field elements used in arithmetic and isogeny computations - """ - if self._precomputation is None: - a, b, c, d = self.null_point().coords() - - # Technically this computes 4A^2, 4B^2, ... - # but as we take quotients this doesnt matter - # Cost: 4S - AA, BB, CC, DD = self.squared_theta() - - # Precomputed constants for addition and doubling - b_inv, c_inv, d_inv, BB_inv, CC_inv, DD_inv = batch_inversion([ - b, c, d, BB, CC, DD] - ) - - y0 = a * b_inv - z0 = a * c_inv - t0 = a * d_inv - - Y0 = AA * BB_inv - Z0 = AA * CC_inv - T0 = AA * DD_inv - - self._precomputation = (y0, z0, t0, Y0, Z0, T0) - return self._precomputation - - def __call__(self, coords): - return self._point(self, coords) - - def base_change_struct(self,N): - null_coords=self.null_point().coords() - new_null_coords=apply_base_change_theta_dim2(N,null_coords) - return ThetaStructure(new_null_coords) - - def base_change_coords(self,N,P): - coords=P.coords() - new_coords=apply_base_change_theta_dim2(N,coords) - return self.__call__(new_coords) - - -# =================================================== # -# Class for Product Theta Structure (level-2) # -# =================================================== # - - -class ProductThetaStructureDim2(ThetaStructureDim2): - def __init__(self,*args): - r"""Defines the product theta structure at level 2 of 2 elliptic curves. - - Input: Either - - 2 theta structures of dimension 1: T0, T1; - - 2 elliptic curves: E0, E1. - - 2 elliptic curves E0, E1 and their respective canonical 4-torsion basis B0, B1. - """ - if len(args)==2: - theta_structures=list(args) - for k in range(2): - if not isinstance(theta_structures[k],ThetaStructureDim1): - theta_structures[k]=ThetaStructureDim1(theta_structures[k]) - elif len(args)==4: - theta_structures=[ThetaStructureDim1(args[k],args[2+k][0],args[2+k][1]) for k in range(2)] - else: - raise ValueError("2 or 4 arguments expected but {} were given.\nYou should enter a list of 2 elliptic curves or ThetaStructureDim1\nor a list of 2 elliptic curves with a 4-torsion basis for each of them.".format(len(args))) - - self._theta_structures=theta_structures - - null_point=product_theta_point(theta_structures[0].zero().coords(),theta_structures[1].zero().coords()) - - ThetaStructureDim2.__init__(self,null_point) - - def product_theta_point(self,theta_points): - t0,t1=theta_points[0].coords() - u0,u1=theta_points[1].coords() - return self._point(self,[t0*u0,t1*u0,t0*u1,t1*u1]) - - def __call__(self,point): - if isinstance(point,TuplePoint): - theta_points=[] - theta_structures=self._theta_structures - for i in range(2): - theta_points.append(theta_structures[i](point[i])) - return self.product_theta_point(theta_points) - else: - return self._point(self,point) - - def to_theta_points(self,P): - coords=P.coords() - theta_coords=[(coords[0],coords[1]),(coords[1],coords[3])] - theta_points=[self._theta_structures[i](theta_coords[i]) for i in range(2)] - return theta_points - - def to_tuple_point(self,P): - theta_points=self.to_theta_points(P) - montgomery_points=[self._theta_structures[i].to_montgomery_point(theta_points[i]) for i in range(2)] - return TuplePoint(montgomery_points) - - -# ======================================= # -# Class for Theta Point (level-2) # -# ======================================= # - - -class ThetaPointDim2: - """ - A Theta Point in the level-2 Theta Structure is defined with four projective - coordinates - - We cannot perform arbitrary arithmetic, but we can compute doubles and - differential addition, which like x-only points on the Kummer line, allows - for scalar multiplication - """ - - def __init__(self, parent, coords): - if not isinstance(parent, ThetaStructureDim2) and not isinstance(parent, ProductThetaStructureDim2): - raise ValueError - - self._parent = parent - self._coords = tuple(coords) - - self._hadamard = None - self._squared_theta = None - - def parent(self): - """ - Return the parent of the element, of type ThetaStructureDim2 - """ - return self._parent - - def theta(self): - """ - Return the parent theta structure of this ThetaPointDim2""" - return self.parent() - - def coords(self): - """ - Return the projective coordinates of the ThetaPointDim2 - """ - return self._coords - - def is_zero(self): - """ - An element is zero if it is equivalent to the null point of the parent - ThetaStrcuture - """ - return self == self.parent().zero() - - @staticmethod - def to_hadamard(x_00, x_10, x_01, x_11): - """ - Compute the Hadamard transformation of four coordinates, using recursive - formula. - """ - x_00, x_10 = (x_00 + x_10, x_00 - x_10) - x_01, x_11 = (x_01 + x_11, x_01 - x_11) - return x_00 + x_01, x_10 + x_11, x_00 - x_01, x_10 - x_11 - - def hadamard(self): - """ - Compute the Hadamard transformation of this element - """ - if self._hadamard is None: - self._hadamard = self.to_hadamard(*self.coords()) - return self._hadamard - - @staticmethod - def to_squared_theta(x, y, z, t): - """ - Square the coordinates and then compute the Hadamard transform of the - input - """ - return ThetaPointDim2.to_hadamard(x * x, y * y, z * z, t * t) - - def squared_theta(self): - """ - Compute the Squared Theta transformation of this element - which is the square operator followed by Hadamard. - """ - if self._squared_theta is None: - self._squared_theta = self.to_squared_theta(*self.coords()) - return self._squared_theta - - def double(self): - """ - Computes [2]*self - - NOTE: Assumes that no coordinate is zero - - Cost: 8S 6M - """ - # If a,b,c,d = 0, then the codomain of A->A/K_2 is a product of - # elliptic curves with a non product theta structure. - # Unless we are very unlucky, A/K_1 will not be in this case, so we - # just need to Hadamard, double, and Hadamard inverse - # If A,B,C,D=0 then the domain itself is a product of elliptic - # curves with a non product theta structure. The Hadamard transform - # will not change this, we need a symplectic change of variable - # that puts us back in a product theta structure - y0, z0, t0, Y0, Z0, T0 = self.parent()._arithmetic_precomputation() - - # Temp coordinates - # Cost 8S 3M - xp, yp, zp, tp = self.squared_theta() - xp = xp**2 - yp = Y0 * yp**2 - zp = Z0 * zp**2 - tp = T0 * tp**2 - - # Final coordinates - # Cost 3M - X, Y, Z, T = self.to_hadamard(xp, yp, zp, tp) - X = X - Y = y0 * Y - Z = z0 * Z - T = t0 * T - - coords = (X, Y, Z, T) - return self._parent(coords) - - def diff_addition(P, Q, PQ): - """ - Given the theta points of P, Q and P-Q computes the theta point of - P + Q. - - NOTE: Assumes that no coordinate is zero - - Cost: 8S 17M - """ - # Extract out the precomputations - Y0, Z0, T0 = P.parent()._arithmetic_precomputation()[-3:] - - # Transform with the Hadamard matrix and multiply - # Cost: 8S 7M - p1, p2, p3, p4 = P.squared_theta() - q1, q2, q3, q4 = Q.squared_theta() - - xp = p1 * q1 - yp = Y0 * p2 * q2 - zp = Z0 * p3 * q3 - tp = T0 * p4 * q4 - - # Final coordinates - PQx, PQy, PQz, PQt = PQ.coords() - - # Note: - # We replace the four divisions by - # PQx, PQy, PQz, PQt by 10 multiplications - # Cost: 10M - PQxy = PQx * PQy - PQzt = PQz * PQt - - X, Y, Z, T = P.to_hadamard(xp, yp, zp, tp) - X = X * PQzt * PQy - Y = Y * PQzt * PQx - Z = Z * PQxy * PQt - T = T * PQxy * PQz - - coords = (X, Y, Z, T) - return P.parent()(coords) - - def scale(self, n): - """ - Scale all coordinates of the ThetaPointDim2 by `n` - """ - x, y, z, t = self.coords() - if not isinstance(n, RingElement): - raise ValueError(f"Cannot scale by element {n} of type {type(n)}") - scaled_coords = (n * x, n * y, n * z, n * t) - return self._parent(scaled_coords) - - def double_iter(self, m): - """ - Compute [2^n] Self - - NOTE: Assumes that no coordinate is zero at any point during the doubling - """ - if not isinstance(m, Integer): - try: - m = Integer(m) - except: - raise TypeError(f"Cannot coerce input scalar {m = } to an integer") - - if m.is_zero(): - return self.parent().zero() - - P1 = self - for _ in range(m): - P1 = P1.double() - return P1 - - def __mul__(self, m): - """ - Uses Montgomery ladder to compute [m] Self - - NOTE: Assumes that no coordinate is zero at any point during the doubling - """ - # Make sure we're multiplying by something value - if not isinstance(m, (int, Integer)): - try: - m = Integer(m) - except: - raise TypeError(f"Cannot coerce input scalar {m = } to an integer") - - # If m is zero, return the null point - if not m: - return self.parent().zero() - - # We are with ±1 identified, so we take the absolute value of m - m = abs(m) - - P0, P1 = self, self - P2 = P1.double() - # If we are multiplying by two, the chain stops here - if m == 2: - return P2 - - # Montgomery double and add. - for bit in bin(m)[3:]: - Q = P2.diff_addition(P1, P0) - if bit == "1": - P2 = P2.double() - P1 = Q - else: - P1 = P1.double() - P2 = Q - - return P1 - - def __rmul__(self, m): - return self * m - - def __imul__(self, m): - self = self * m - return self - - def __eq__(self, other): - """ - Check the quality of two ThetaPoints. Note that as this is a - projective equality, we must be careful for when certain coefficients may - be zero. - """ - if not isinstance(other, ThetaPointDim2): - return False - - a1, b1, c1, d1 = self.coords() - a2, b2, c2, d2 = other.coords() - - if d1 != 0 or d2 != 0: - return all([a1 * d2 == a2 * d1, b1 * d2 == b2 * d1, c1 * d2 == c2 * d1]) - elif c1 != 0 or c2 != 0: - return all([a1 * c2 == a2 * c1, b1 * c2 == b2 * c1]) - elif b1 != 0 or b2 != 0: - return a1 * b2 == a2 * b1 - else: - return True - - def __repr__(self): - return f"Theta point with coordinates: {self.coords()}" diff --git a/theta_lib/theta_structures/Theta_dim4.py b/theta_lib/theta_structures/Theta_dim4.py deleted file mode 100644 index 7851c56..0000000 --- a/theta_lib/theta_structures/Theta_dim4.py +++ /dev/null @@ -1,351 +0,0 @@ -from sage.all import * -from sage.structure.element import get_coercion_model - -from .theta_helpers_dim4 import ( - hadamard, - batch_inversion, - product_theta_point_dim4, - product_to_theta_points_dim4, - product_theta_point_dim2_dim4, - product_to_theta_points_dim4_dim2, - act_point, - squared, -) -from .Theta_dim1 import ThetaStructureDim1 -from .Tuple_point import TuplePoint -from ..basis_change.base_change_dim4 import ( - apply_base_change_theta_dim4, - random_symplectic_matrix, - base_change_theta_dim4, -) - -cm = get_coercion_model() - - -class ThetaStructureDim4: - def __init__(self,null_point,null_point_dual=None,inv_null_point_dual_sq=None): - r""" - INPUT: - - null_point: theta-constants. - - inv_null_point_dual_sq: inverse of the squares of dual theta-constants, if provided - (meant to prevent duplicate computation, since this data is already computed when the - codomain of an isogeny is computed). - """ - if not len(null_point) == 16: - raise ValueError("Entry null_point should have 16 coordinates.") - - self._base_ring = cm.common_parent(*(c.parent() for c in null_point)) - self._point = ThetaPointDim4 - self._null_point = self._point(self, null_point) - self._null_point_dual=null_point_dual - self._inv_null_point=None - self._inv_null_point_dual_sq=inv_null_point_dual_sq - - def null_point(self): - """ - """ - return self._null_point - - def null_point_dual(self): - if self._null_point_dual==None: - self._null_point_dual=hadamard(self._null_point.coords()) - return self._null_point_dual - - def base_ring(self): - """ - """ - return self._base_ring - - def zero(self): - """ - """ - return self.null_point() - - def zero_dual(self): - return self.null_point_dual() - - def __repr__(self): - return f"Theta structure over {self.base_ring()} with null point: {self.null_point()}" - - def __call__(self,coords): - return self._point(self,coords) - - def act_null(self,I,J): - r""" - Point of 2-torsion. - - INPUT: - - I, J: two 4-tuples of indices in {0,1}. - - OUTPUT: the action of (I,\chi_J) on the theta null point given by: - (I,\chi_J)*P=(\chi_J(I+K)^{-1}P[I+K])_K - """ - return self.null_point().act_point(I,J) - - def is_K2(self,B): - r""" - Given a symplectic decomposition A[2]=K_1\oplus K_2 canonically - induced by the theta-null point, determines if B is the canonical - basis of K_2 given by act_nul(0,\delta_i)_{1\leq i\leq 4}. - - INPUT: - - B: Basis of 4 points of 2-torsion. - - OUTPUT: Boolean True if and only if B is the canonical basis of K_2. - """ - I0=(0,0,0,0) - if B[0]!=self.act_null(I0,(1,0,0,0)): - return False - if B[1]!=self.act_null(I0,(0,1,0,0)): - return False - if B[2]!=self.act_null(I0,(0,0,1,0)): - return False - if B[3]!=self.act_null(I0,(0,0,0,1)): - return False - return True - - def base_change_struct(self,N): - null_coords=self.null_point().coords() - new_null_coords=apply_base_change_theta_dim4(N,null_coords) - return ThetaStructureDim4(new_null_coords) - - def base_change_coords(self,N,P): - coords=P.coords() - new_coords=apply_base_change_theta_dim4(N,coords) - return self.__call__(new_coords) - - #@cached_method - def _arithmetic_precomputation(self): - r""" - Initializes the precomputation containing the inverse of the theta-constants in standard - and dual (Hadamard transformed) theta-coordinates. Assumes no theta-constant is zero. - """ - O=self.null_point() - if all([O[k]!=0 for k in range(16)]): - self._inv_null_point=batch_inversion(O.coords()) - if self._inv_null_point_dual_sq==None: - U_chi_0_sq=hadamard(squared(O.coords())) - if all([U_chi_0_sq[k]!=0 for k in range(16)]): - self._inv_null_point_dual_sq=batch_inversion(U_chi_0_sq) - self._arith_base_change=False - else: - self._arith_base_change=True - self._arithmetic_base_change() - #print("Zero dual theta constants.\nRandom symplectic base change is being used for duplication.\nDoublings are more costly than expected.") - else: - self._arith_base_change=True - self._arithmetic_base_change() - #print("Zero theta constants.\nRandom symplectic base change is being used for duplication.\nDoublings are more costly than expected.") - - def _arithmetic_base_change(self,max_iter=50): - F=self._base_ring - if F.degree() == 2: - i=self._base_ring.gen() - else: - assert(F.degree() == 1) - Fp2 = GF((F.characteristic(), 2), name='i', modulus=var('x')**2 + 1) - i=Fp2.gen() - - count=0 - O=self.null_point() - while count<max_iter: - count+=1 - M=random_symplectic_matrix(4) - N=base_change_theta_dim4(M,i) - - NO=apply_base_change_theta_dim4(N,O) - NU_chi_0_sq=hadamard(squared(NO)) - - if all([NU_chi_0_sq[k]!=0 for k in range(16)]) and all([NO[k]!=0 for k in range(16)]): - self._arith_base_change_matrix=N - self._arith_base_change_matrix_inv=N.inverse() - self._inv_null_point=batch_inversion(NO) - self._inv_null_point_dual_sq=batch_inversion(NU_chi_0_sq) - break - - def has_suitable_doubling(self): - O=self.null_point() - UO=hadamard(O.coords()) - if all([O[k]!=0 for k in range(16)]) and all([UO[k]!=0 for k in range(16)]): - return True - else: - return False - - def hadamard(self): - return ThetaStructureDim4(self.null_point_dual()) - - def hadamard_change_coords(self,P): - new_coords=hadamard(P) - return self.__call__(new_coords) - - -class ProductThetaStructureDim1To4(ThetaStructureDim4): - def __init__(self,*args): - r"""Defines the product theta structure at level 2 of 4 elliptic curves. - - Input: Either - - 4 theta structures of dimension 1: T0, T1, T2, T3; - - 4 elliptic curves: E0, E1, E2, E3. - - 4 elliptic curves E0, E1, E2, E3 and their respective canonical 4-torsion basis B0, B1, B2, B3. - """ - if len(args)==4: - theta_structures=list(args) - for k in range(4): - if not isinstance(theta_structures[k],ThetaStructureDim1): - try: - theta_structures[k]=ThetaStructureDim1(theta_structures[k]) - except: - pass - elif len(args)==8: - theta_structures=[ThetaStructureDim1(args[k],args[4+k][0],args[4+k][1]) for k in range(4)] - else: - raise ValueError("4 or 8 arguments expected but {} were given.\nYou should enter a list of 4 elliptic curves or ThetaStructureDim1\nor a list of 4 elliptic curves with a 4-torsion basis for each of them.".format(len(args))) - - self._theta_structures=theta_structures - - null_point=product_theta_point_dim4(theta_structures[0].zero().coords(),theta_structures[1].zero().coords(), - theta_structures[2].zero().coords(),theta_structures[3].zero().coords()) - - ThetaStructureDim4.__init__(self,null_point) - - def product_theta_point(self,theta_points): - return self._point(self,product_theta_point_dim4(theta_points[0].coords(),theta_points[1].coords(), - theta_points[2].coords(),theta_points[3].coords())) - - def __call__(self,point): - if isinstance(point,TuplePoint): - theta_points=[] - theta_structures=self._theta_structures - for i in range(4): - theta_points.append(theta_structures[i](point[i])) - return self.product_theta_point(theta_points) - else: - return self._point(self,point) - - def to_theta_points(self,P): - theta_coords=product_to_theta_points_dim4(P) - theta_points=[self._theta_structures[i](theta_coords[i]) for i in range(4)] - return theta_points - - def to_tuple_point(self,P): - theta_points=self.to_theta_points(P) - montgomery_points=[self._theta_structures[i].to_montgomery_point(theta_points[i]) for i in range(4)] - return TuplePoint(montgomery_points) - -class ProductThetaStructureDim2To4(ThetaStructureDim4): - def __init__(self,theta1,theta2): - self._theta_structures=(theta1,theta2) - - null_point=product_theta_point_dim2_dim4(theta1.zero().coords(),theta2.zero().coords()) - - ThetaStructureDim4.__init__(self,null_point) - - def product_theta_point(self,P1,P2): - return self._point(self,product_theta_point_dim2_dim4(P1.coords(),P2.coords())) - - def __call__(self,point): - return self._point(self,point) - - def to_theta_points(self,P): - theta_coords=product_to_theta_points_dim4_dim2(P) - theta_points=[self._theta_structures[i](theta_coords[i]) for i in range(2)] - return theta_points - -class ThetaPointDim4: - def __init__(self, parent, coords): - """ - """ - if not isinstance(parent, ThetaStructureDim4): - raise ValueError("Entry parent should be a ThetaStructureDim4 object.") - - self._parent = parent - self._coords = tuple(coords) - - def parent(self): - """ - """ - return self._parent - - def theta(self): - """ - """ - return self.parent() - - def coords(self): - """ - """ - return self._coords - - def is_zero(self): - """ - """ - return self == self.parent().zero() - - def __eq__(self, other): - P=self.coords() - Q=other.coords() - - k0=0 - while k0<15 and P[k0]==0: - k0+=1 - - for l in range(16): - if P[l]*Q[k0]!=Q[l]*P[k0]: - return False - return True - - def __repr__(self): - return f"Theta point with coordinates: {self.coords()}" - - def __getitem__(self,i): - return self._coords[i] - - def scale(self,lamb): - if lamb==0: - raise ValueError("Entry lamb should be non-zero.") - - P=self.coords() - return self._parent([lamb*x for x in P]) - - def act_point(self,I,J): - r""" - Translation by a point of 2-torsion. - - Input: - - I, J: two 4-tuples of indices in {0,1}. - - Output: the action of (I,\chi_J) on P given by: - (I,\chi_J)*P=(\chi_J(I+K)^{-1}P[I+K])_K - """ - return self._parent(act_point(self._coords,I,J)) - - def double(self): - ## This formula is projective. - ## Works only when theta constants are non-zero. - P=self.coords() - if self.parent()._inv_null_point==None or self.parent()._inv_null_point_dual_sq==None: - self.parent()._arithmetic_precomputation() - inv_O,inv_U_chi_0_sq=self.parent()._inv_null_point,self.parent()._inv_null_point_dual_sq - - if self.parent()._arith_base_change: - P=apply_base_change_theta_dim4(self.parent()._arith_base_change_matrix,P) - - U_chi_P=squared(hadamard(squared(P))) - for chi in range(16): - U_chi_P[chi]*=inv_U_chi_0_sq[chi] - - theta_2P = list(hadamard(U_chi_P)) - for i in range(16): - theta_2P[i] *= inv_O[i] - - if self.parent()._arith_base_change: - theta_2P=apply_base_change_theta_dim4(self.parent()._arith_base_change_matrix_inv,theta_2P) - - return self._parent(theta_2P) - - def double_iter(self,n): - ## Computes 2**n*self - Q=self - for i in range(n): - Q=Q.double() - return Q diff --git a/theta_lib/theta_structures/Tuple_point.py b/theta_lib/theta_structures/Tuple_point.py deleted file mode 100644 index aac248b..0000000 --- a/theta_lib/theta_structures/Tuple_point.py +++ /dev/null @@ -1,106 +0,0 @@ -from sage.all import * -from ..utilities.discrete_log import weil_pairing_pari - -class TuplePoint: - def __init__(self,*args): - if len(args)==1: - self._points=list(args[0]) - else: - self._points=list(args) - - def points(self): - return self._points - - def parent_curves(self): - return [x.curve() for x in self._points] - - def parent_curve(self,i): - return self._points[i].curve() - - def n_points(self): - return len(self._points) - - def is_zero(self): - return all([self._points[i]==0 for i in range(self.n_points())]) - - def __repr__(self): - return str(self._points) - - def __getitem__(self,i): - return self._points[i] - - def __setitem__(self,i,P): - self._points[i]=P - - def __eq__(self,other): - n_self=self.n_points() - n_other=self.n_points() - return n_self==n_other and all([self._points[i]==other._points[i] for i in range(n_self)]) - - def __add__(self,other): - n_self=self.n_points() - n_other=self.n_points() - - if n_self!=n_other: - raise ValueError("Cannot add TuplePoint of distinct lengths {} and {}.".format(n_self,n_other)) - - points=[] - for i in range(n_self): - points.append(self._points[i]+other._points[i]) - return self.__class__(points) - - def __sub__(self,other): - n_self=self.n_points() - n_other=self.n_points() - - if n_self!=n_other: - raise ValueError("Cannot substract TuplePoint of distinct lengths {} and {}.".format(n_self,n_other)) - - points=[] - for i in range(n_self): - points.append(self._points[i]-other._points[i]) - return self.__class__(points) - - def __neg__(self): - n_self=self.n_points() - points=[] - for i in range(n_self): - points.append(-self._points[i]) - return self.__class__(points) - - def __mul__(self,m): - n_self=self.n_points() - points=[] - for i in range(n_self): - points.append(m*self._points[i]) - return self.__class__(points) - - def __rmul__(self,m): - return self*m - - def double_iter(self,n): - result=self - for i in range(n): - result=2*result - return result - - def weil_pairing(self,other,n): - n_self=self.n_points() - n_other=self.n_points() - - if n_self!=n_other: - raise ValueError("Cannot compute the Weil pairing of TuplePoint of distinct lengths {} and {}.".format(n_self,n_other)) - - zeta=1 - for i in range(n_self): - zeta*=weil_pairing_pari(self._points[i],other._points[i],n) - - return zeta - - - - - - - - diff --git a/theta_lib/theta_structures/montgomery_theta.py b/theta_lib/theta_structures/montgomery_theta.py deleted file mode 100644 index 6dddb2b..0000000 --- a/theta_lib/theta_structures/montgomery_theta.py +++ /dev/null @@ -1,68 +0,0 @@ -from sage.all import * - -def torsion_to_theta_null_point(P): - r=P[0] - s=P[2] - return (r+s,r-s) - -def montgomery_point_to_theta_point(O,P): - if P[0]==P[2]==0: - return O - else: - a,b=O - return (a*(P[0]-P[2]),b*(P[0]+P[2])) - -def theta_point_to_montgomery_point(E,O,P,twist=False): - a,b=O - - x=a*P[1]+b*P[0] - z=a*P[1]-b*P[0] - - if twist: - x=-x - - if z==0: - return E(0) - else: - x=x/z - - a_inv=E.a_invariants() - - A =a_inv[1] - if a_inv != (0,A,0,1,0): - raise ValueError("The elliptic curve E is not in the Montgomery model.") - - y2=x*(x**2+A*x+1) - if not is_square(y2): - raise ValueError("The Montgomery point is not on the base field.") - else: - y=y2.sqrt() - return E([x,y,1]) - -def lift_kummer_montgomery_point(E,x,z=1): - if z==0: - return E(0) - elif z!=1: - x=x/z - - a_inv=E.a_invariants() - - A =a_inv[1] - if a_inv != (0,A,0,1,0): - raise ValueError("The elliptic curve E is not in the Montgomery model.") - - y2=x*(x**2+A*x+1) - if not is_square(y2): - raise ValueError("The Montgomery point is not on the base field.") - else: - y=y2.sqrt() - return E([x,y,1]) - -def null_point_to_montgomery_coeff(a,b): - return -2*(a**4+b**4)/(a**4-b**4) - - - - - - diff --git a/theta_lib/theta_structures/theta_helpers_dim2.py b/theta_lib/theta_structures/theta_helpers_dim2.py deleted file mode 100644 index a57789d..0000000 --- a/theta_lib/theta_structures/theta_helpers_dim2.py +++ /dev/null @@ -1,32 +0,0 @@ -from sage.all import * - -def batch_inversion(L): - r"""Does n inversions in 3(n-1)M+1I. - - Input: - - L: list of elements to invert. - - Output: - - [1/x for x in L] - """ - # Given L=[a0,...,an] - # Computes multiples=[a0, a0.a1, ..., a0...an] - multiples=[L[0]] - for ai in L[1:]: - multiples.append(multiples[-1]*ai) - - # Computes inverses=[1/(a0...an),...,1/a0] - inverses=[1/multiples[-1]] - for i in range(1,len(L)): - inverses.append(inverses[-1]*L[-i]) - - # Finally computes [1/a0,...,1/an] - result=[inverses[-1]] - for i in range(2,len(L)+1): - result.append(inverses[-i]*multiples[i-2]) - return result - -def product_theta_point(theta1,theta2): - a,b=theta1 - c,d=theta2 - return [a*c,b*c,a*d,b*d] diff --git a/theta_lib/theta_structures/theta_helpers_dim4.py b/theta_lib/theta_structures/theta_helpers_dim4.py deleted file mode 100644 index 29599b0..0000000 --- a/theta_lib/theta_structures/theta_helpers_dim4.py +++ /dev/null @@ -1,259 +0,0 @@ -from sage.all import * -import itertools - -## Index management -@cached_function -def index_to_multindex(k): - r""" - Input: - - k: integer between 0 and 15. - - Output: binary decomposition of k. - """ - L_ind=[] - l=k - for i in range(4): - L_ind.append(l%2) - l=l//2 - return tuple(L_ind) - -@cached_function -def multindex_to_index(*args): - r""" - Input: 4 elements i0,i1,i2,i3 in {0,1}. - - Output: k=i0+2*i1+4*i2+8*i3. - """ - if len(args)==4: - i0,i1,i2,i3=args - else: - i0,i1,i2,i3=args[0] - return i0+2*i1+4*i2+8*i3 - -@cached_function -def scal_prod(i,j): - r""" - Input: Two integers i and j in {0,...,15}. - - Output: Scalar product of the bits of i and j mod 2. - """ - return (int(i)&int(j)).bit_count()%2 - -def act_point(P,I,J): - r""" - Input: - - P: a point with 16 coordinates. - - I, J: two 4-tuples of indices in {0,1}. - - Output: the action of (I,\chi_J) on P given by: - (I,\chi_J)*P=(\chi_J(I+K)^{-1}P[I+K])_K - """ - Q=[] - i=multindex_to_index(I) - j=multindex_to_index(J) - for k in range(16): - ipk=i^k - Q.append((-1)**scal_prod(ipk,j)*P[ipk]) - return Q - -## Product of theta points -def product_theta_point_dim4(P0,P1,P2,P3): - # Computes the product theta coordinates of a product of 4 elliptic curves. - P=[0 for k in range(16)] - for i0, i1, i2, i3 in itertools.product([0,1],repeat=4): - P[multindex_to_index(i0,i1,i2,i3)]=P0[i0]*P1[i1]*P2[i2]*P3[i3] - return P - -def product_theta_point_dim2_dim4(P0,P1): - # Computes the product theta coordinates of a product of 2 abelian surfaces. - P=[0 for k in range(16)] - for i0, i1, i2, i3 in itertools.product([0,1],repeat=4): - P[multindex_to_index(i0,i1,i2,i3)]=P0[i0+2*i1]*P1[i2+2*i3] - return P - -## 4-dimensional product Theta point to 1-dimensional Theta points -def product_to_theta_points_dim4(P): - Fp2=P[0].parent() - d_Pi={0:None,1:None,2:None,3:None} - d_index_ratios={0:None,1:None,2:None,3:None}# Index of numertors/denominators to compute the theta points. - for k in range(4): - is_zero=True# theta_1(Pk)=0 - for I in itertools.product([0,1],repeat=3): - J=list(I) - J.insert(k,1) - j=multindex_to_index(*J) - if P[j]!=0: - is_zero=False - d_index_ratios[k]=[j^(2**k),j] - break - if is_zero: - d_Pi[k]=(Fp2(1),Fp2(0)) - L_num=[] - L_denom=[] - d_index_num_denom={} - for k in range(4): - if d_Pi[k]==None:# Point has non-zero coordinate theta_1 - d_index_num_denom[k]=len(L_num) - L_num.append(P[d_index_ratios[k][0]]) - L_denom.append(P[d_index_ratios[k][1]]) - L_denom=batch_inversion(L_denom) - for k in range(4): - if d_Pi[k]==None: - d_Pi[k]=(L_num[d_index_num_denom[k]]*L_denom[d_index_num_denom[k]],Fp2(1)) - return d_Pi[0],d_Pi[1],d_Pi[2],d_Pi[3] - -## 4-dimensional product Theta point to 2-dimensional Theta points -def product_to_theta_points_dim4_dim2(P): - Fp2=P[0].parent() - k0=0# Index of the denominator k0=multindex_to_index(I0,J0) and - # we divide by theta_{I0,J0}=theta_{I0}*theta_{J0}!=0 - while k0<=15 and P[k0]==0: - k0+=1 - i0, j0 = k0%4, k0//4 - inv_theta_k0=1/P[k0] - theta_P1=[] - theta_P2=[] - for i in range(4): - if i==i0: - theta_P1.append(1) - else: - theta_P1.append(P[i+4*j0]*inv_theta_k0) - for j in range(4): - if j==j0: - theta_P2.append(1) - else: - theta_P2.append(P[i0+4*j]*inv_theta_k0) - return theta_P1, theta_P2 - - - - - -## Usual theta transformations -@cached_function -def inv_16(F): - return 1/F(16) - -def hadamard2(x,y): - return (x+y, x-y) - -def hadamard4(x,y,z,t): - x,y=hadamard2(x,y) - z,t=hadamard2(z,t) - return (x+z, y+t, x-z, y-t) - -def hadamard8(a,b,c,d,e,f,g,h): - a,b,c,d=hadamard4(a,b,c,d) - e,f,g,h=hadamard4(e,f,g,h) - return (a+e, b+f, c+g, d+h, a-e, b-f, c-g, d-h) - -def hadamard16(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p): - a,b,c,d,e,f,g,h=hadamard8(a,b,c,d,e,f,g,h) - i,j,k,l,m,n,o,p=hadamard8(i,j,k,l,m,n,o,p) - return (a+i, b+j, c+k, d+l, e+m, f+n, g+o, h+p, a-i, b-j, c-k, d-l, e-m, f-n, g-o, h-p) - -def hadamard_inline(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p): - a,b=a+b,a-b - c,d=c+d,c-d - e,f=e+f,e-f - g,h=g+h,g-h - i,j=i+j,i-j - k,l=k+l,k-l - m,n=m+n,m-n - o,p=o+p,o-p - a,b,c,d=a+c,b+d,a-c,b-d - e,f,g,h=e+g,f+h,e-g,f-h - i,j,k,l=i+k,j+l,i-k,j-l - m,n,o,p=m+o,n+p,m-o,n-p - a,b,c,d,e,f,g,h=a+e, b+f, c+g, d+h, a-e, b-f, c-g, d-h - i,j,k,l,m,n,o,p=i+m,j+n,k+o,l+p,i-m,j-n,k-o,l-p - return (a+i, b+j, c+k, d+l, e+m, f+n, g+o, h+p, a-i, b-j, c-k, d-l, e-m, f-n, g-o, h-p) - -def hadamard(P): - return hadamard16(*P) - #return hadamard_inline(*P) - -def hadamard_inverse(P): - H_inv_P=[] - C=inv_16(P[0].parent()) - for j in range(16): - HP.append(0) - for k in range(16): - if scal_prod(k,j)==0: - H_inv_P[j]+=P[k] - else: - H_inv_P[j]-=P[k] - H_inv_P[j]=H_inv_P[j]*C - return H_inv_P - -def squared(P): - return [x**2 for x in P] - -def batch_inversion(L): - r"""Does n inversions in 3(n-1)M+1I. - - Input: - - L: list of elements to invert. - - Output: - - [1/x for x in L] - """ - # Given L=[a0,...,an] - # Computes multiples=[a0, a0.a1, ..., a0...an] - multiples=[L[0]] - for ai in L[1:]: - multiples.append(multiples[-1]*ai) - - # Computes inverses=[1/(a0...an),...,1/a0] - inverses=[1/multiples[-1]] - for i in range(1,len(L)): - inverses.append(inverses[-1]*L[-i]) - - # Finally computes [1/a0,...,1/an] - result=[inverses[-1]] - for i in range(2,len(L)+1): - result.append(inverses[-i]*multiples[i-2]) - return result - - -## Functions to handle zero theta coordinates -def find_zeros(P): - L_ind_zeros=[] - for i in range(16): - if P[i]==0: - L_ind_zeros.append(i) - return L_ind_zeros - -def find_translates(L_ind_zeros): - L_ind_non_zero=[] - L_ind_origin=L_ind_zeros.copy() - - for i in range(16): - if i not in L_ind_zeros: - L_ind_non_zero.append(i) - - L_trans=[] - while L_ind_origin!=[]: - n_target_max=0 - ind_trans_max=0 - for k in range(16): - trans=[x^k for x in L_ind_origin] - n_target=0 - for y in trans: - if y in L_ind_non_zero: - n_target+=1 - if n_target>n_target_max: - n_target_max=n_target - ind_trans_max=k - L_trans.append(ind_trans_max) - L_ind_remove=[] - for x in L_ind_origin: - if x^ind_trans_max in L_ind_non_zero: - L_ind_remove.append(x) - for x in L_ind_remove: - L_ind_origin.remove(x) - return L_trans - - - - |