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