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