diff options
Diffstat (limited to 'theta_lib')
25 files changed, 5822 insertions, 0 deletions
diff --git a/theta_lib/basis_change/base_change_dim2.py b/theta_lib/basis_change/base_change_dim2.py new file mode 100644 index 0000000..4994529 --- /dev/null +++ b/theta_lib/basis_change/base_change_dim2.py @@ -0,0 +1,150 @@ +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 new file mode 100644 index 0000000..aaf685f --- /dev/null +++ b/theta_lib/basis_change/base_change_dim4.py @@ -0,0 +1,152 @@ +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 new file mode 100644 index 0000000..e1c3d1f --- /dev/null +++ b/theta_lib/basis_change/canonical_basis_dim1.py @@ -0,0 +1,76 @@ +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 new file mode 100644 index 0000000..e6de2e2 --- /dev/null +++ b/theta_lib/basis_change/kani_base_change.py @@ -0,0 +1,975 @@ +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 new file mode 100644 index 0000000..66030e2 --- /dev/null +++ b/theta_lib/isogenies/Kani_clapoti.py @@ -0,0 +1,258 @@ +from sage.all import * +import itertools + +from ..basis_change.kani_base_change import clapoti_cob_splitting_matrix +from ..basis_change.base_change_dim4 import base_change_theta_dim4 +from ..theta_structures.Tuple_point import TuplePoint +from ..theta_structures.montgomery_theta import null_point_to_montgomery_coeff, theta_point_to_montgomery_point +from ..theta_structures.theta_helpers_dim4 import product_to_theta_points_dim4 +from ..utilities.supersingular import torsion_basis_to_Fp_rational_point +from .Kani_gluing_isogeny_chain_dim4 import KaniClapotiGluing +from .isogeny_chain_dim4 import IsogenyChainDim4 + +class KaniClapotiIsog(IsogenyChainDim4): + r"""Class representing the 4-dimensional isogeny obtained via Kani's lemma F: Eu^2*Ev^2 --> Ea^2*A + where Ea=[\mf{a}]*E is the result of the ideal class group action by \mf{a} when given relevant + constants and torsion point information. + + INPUT: + - Pu, Qu = phi_u(P, Q)\in Eu; + - Pv, Qv = phi_v*phi_{ck}*\hat{\phi}_{bk}(P, Q)\in Ev; + - gu, xu, yu, gv, xv, yv, Nbk, Nck, e: positive integers; + where: + * gu(xu^2+yu^2)Nbk+gv(xv^2+yv^2)Nck=2^e; + * gcd(u*Nbk,v*Nck)=1 with u:=gu(xu^2+yu^2) and v:=gv(xv^2+yv^2); + * xu and xv are odd and yu and yv are even; + * \mf{b}=\mf{be}*\mf{bk} is a product of ideals of norms Nbe and Nbk respectively, + where Nbe is a product of small Elkies primes; + * \mf{c}=\mf{ce}*\mf{ck} is a product of ideals of norms Nce and Nck respectively, + where Nbe is a product of small Elkies primes; + * phi_{bk}: E --> E1 and phi_{ck}: E --> E2 are induced by the action of + ideals \mf{bk} and \mf{ck} respectively; + * <P,Q>=E_1[2^{e+2}]; + * phi_u: E1 --> Eu and phi_v: E2 --> Ev are gu and gv-isogenies respectively. + + OUTPUT: F: Eu^2*Ev^2 --> Ea^2*A is the isogeny: + + F := [[Phi_{bp}*\tilde{Phi}_u, Phi_{cp}*\tilde{Phi}_v], + [-Psi, \tilde{Phi}]] + + obtained from the Kani isogeny diamond: + + A --------------------Phi------------------> Ev^2 + ^ ^ + | | + | Phi_v + | | + Psi E2^2 + | ^ + | | + | \tilde{Phi}_{ck} + | | + Eu^2 --\tilde{Phi}_{u}--> E1^2 --Phi_{bk}--> Ea^2 + + where Phi_{bk}:=Diag(phi_{bk},phi_{bk}), Phi_{ck}:=Diag(phi_{ck},phi_{ck}), + + Phi_u := [[xu, -yu], + [yu, xu]] * Diag(phi_u,phi_u) + + Phi_v := [[xv, -yv], + [yv, xv]] * Diag(phi_v,phi_v) + """ + + def __init__(self,points,integers,strategy=None): + gu,xu,yu,gv,xv,yv,Nbk,Nck,e = integers + Pu,Qu,Pv,Qv = points + if gu*(xu**2+yu**2)*Nbk+gv*(xv**2+yv**2)*Nck!=2**e: + raise ValueError("Wrong parameters: gu(xu^2+yu^2)Nbk + gv(xv^2+yv^2)Nck != 2^e") + if gcd(ZZ(gu*(xu**2+yu**2)*Nbk),ZZ(gv*(xv**2+yv**2)*Nck))!=1: + raise ValueError("Non coprime parameters: gcd(gu(xu^2+yu^2)Nbk, gv(xv^2+yv^2)Nck) != 1") + if xu%2==0: + xu,yu=yu,xu + if xv%2==0: + xv,yv=yv,xv + + self.Eu = Pu.curve() + self.Ev = Pv.curve() + Fp2 = parent(Pu[0]) + Fp = Fp2.base_ring() + + # Number of dimension 2 steps before dimension 4 gluing Am^2-->B + m=valuation(xv*yu-xu*yv,2) + integers=[gu,xu,yu,gv,xv,yv,Nbk,Nck,e,m] + + points_mp3=[(2**(e-m-1))*P for P in points] + points_mp2=[2*P for P in points_mp3] + points_4=[(2**m)*P for P in points_mp2] + + self.Ru_Fp = torsion_basis_to_Fp_rational_point(self.Eu,points_4[0],points_4[1],4) + + self.gluing_isogeny_chain = KaniClapotiGluing(points_mp3,points_mp2,points_4,integers, coerce=Fp) + + xuNbk = xu*Nbk + yuNbk = yu*Nbk + two_ep2 = 2**(e+2) + inv_Nbk = inverse_mod(Nbk,two_ep2) + u = gu*(xu**2+yu**2) + inv_u = inverse_mod(u,4) + lambxu = ((1-2**e*inv_u*inv_Nbk)*xu)%two_ep2 + lambyu = ((1-2**e*inv_u*inv_Nbk)*yu)%two_ep2 + xv_Nbk = (xv*inv_Nbk)%two_ep2 + yv_Nbk = (yv*inv_Nbk)%two_ep2 + + + B_Kpp = [TuplePoint(xuNbk*Pu,yuNbk*Pu,xv*Pv,yv*Pv), + TuplePoint(-yuNbk*Pu,xuNbk*Pu,-yv*Pv,xv*Pv), + TuplePoint(lambxu*Qu,lambyu*Qu,xv_Nbk*Qv,yv_Nbk*Qv), + TuplePoint(-lambyu*Qu,lambxu*Qu,-yv_Nbk*Qv,xv_Nbk*Qv)] + + IsogenyChainDim4.__init__(self, B_Kpp, self.gluing_isogeny_chain, e, m, splitting=True, strategy=strategy) + + # Splitting + M_split = clapoti_cob_splitting_matrix(integers) + + self.N_split = base_change_theta_dim4(M_split,self.gluing_isogeny_chain.e4) + + self.codomain_product = self._isogenies[-1]._codomain.base_change_struct(self.N_split) + + # Extracting the group action image Ea=[\mathfrak{a}]*E from the codomain Ea^2*E'^2 + self.theta_null_Ea, self.theta_null_Ep, self.Ea, self.Ep = self.extract_montgomery_curve() + + + + def extract_montgomery_curve(self): + + # Computing the theta null point of Ea + null_point=self.codomain_product.zero() + Fp2=parent(null_point[0]) + Fp = Fp2.base_ring() + for i3, i4 in itertools.product([0,1],repeat=2): + if null_point[4*i3+8*i4]!=0: + i30=i3 + i40=i4 + theta_Ea_0=Fp(null_point[4*i3+8*i4]) + theta_Ea_1=Fp(null_point[1+4*i3+8*i4]) + break + for i1, i2 in itertools.product([0,1],repeat=2): + if null_point[i1+2*i2]!=0: + i10=i1 + i20=i2 + theta_Ep_0=Fp(null_point[i1+2*i2]) + theta_Ep_1=Fp(null_point[i1+2*i2+4]) + break + + # Sanity check: is the codomain of F a product of the form Ea^2*E'^2 ? + theta_Ea=[Fp(theta_Ea_0),Fp(theta_Ea_1)] + theta_Ep=[Fp(theta_Ep_0),Fp(theta_Ep_1)] + + theta_Ea2Ep2=[0 for i in range(16)] + for i1,i2,i3,i4 in itertools.product([0,1],repeat=4): + theta_Ea2Ep2[i1+2*i2+4*i3+8*i4]=theta_Ea[i1]*theta_Ea[i2]*theta_Ep[i3]*theta_Ep[i4] + theta_Ea2Ep2=self.codomain_product(theta_Ea2Ep2) + + assert theta_Ea2Ep2.is_zero() + + A_Ep = null_point_to_montgomery_coeff(theta_Ep_0,theta_Ep_1) + Ep = EllipticCurve([0,A_Ep,0,1,0]) + + ## ## Recovering Ea over Fp and not Fp2 + ## self.find_Fp_rational_theta_struct_Ea() + + ## theta_Ea = self.iso_Ea(theta_Ea) + ## A_Ea = null_point_to_montgomery_coeff(theta_Ea[0],theta_Ea[1]) + + ## # Sanity check : the curve Ea should be defined over Fp + ## # assert A_Ea[1] == 0 + + ## # Twisting Ea if necessary: if A_Ea+2 is not a square in Fp, then we take the twist (A_Ea --> -A_Ea) + ## p=self.Eu.base_field().characteristic() + ## self.twist = False + ## if (A_Ea+2)**((p-1)//2)==-1: + ## A_Ea = -A_Ea + ## self.twist = True + + ## Ea = EllipticCurve([0,A_Ea,0,1,0]) + + A = null_point_to_montgomery_coeff(theta_Ea_0, theta_Ea_1) + Ab = null_point_to_montgomery_coeff(theta_Ea_0+theta_Ea_1, theta_Ea_0-theta_Ea_1) + Acan = min([A, -A, Ab, -Ab]) + Acan = A + if (Acan == A or Acan == -A): + # 'Id' corresponds to the point on the twist + self.iso_type = 'Id' + else: + # 'Hadamard' corresponds to the point on the curve + self.iso_type = 'Hadamard' + if ((self.iso_type == 'Hadamard' and not (Acan+2).is_square()) or (self.iso_type == 'Id' and (Acan+2).is_square())): + Acan=-Acan + if (Acan == A or Acan == Ab): + self.twist = False + else: + self.twist = True + Ea = EllipticCurve([0,Acan,0,1,0]) + + # Find the dual null point + return theta_Ea, theta_Ep, Ea, Ep + + def eval_rational_point_4_torsion(self): + T = TuplePoint(self.Ru_Fp,self.Eu(0),self.Ev(0),self.Ev(0)) + + FPu_4 = self.evaluate_isogeny(T) + FPu_4=self.codomain_product.base_change_coords(self.N_split,FPu_4) + FPu_4=product_to_theta_points_dim4(FPu_4) + + return FPu_4[0] + + def find_Fp_rational_theta_struct_Ea(self): + Pa = self.eval_rational_point_4_torsion() + + HPa = (Pa[0]+Pa[1],Pa[0]-Pa[1]) + i = self.Eu.base_field().gen() + self.i = i + iHPa = (Pa[0]+i*Pa[1],Pa[0]-i*Pa[1]) + + if Pa[0]==0 or Pa[1]==0: + self.iso_type='Id' + elif HPa[0]==0 or HPa[1]==0: + self.iso_type='Hadamard' + elif iHPa[0]==0 or iHPa[1]==0: + self.iso_type='iHadamard' + else: + raise ValueError("A rational theta point should be mapped to (0:1) or (1:0) after change of theta coordinates on Ea.") + + def iso_Ea(self,P): + # Change of theta coordinates to obtain Fp-rational theta coordinates on Ea + + if self.iso_type == 'Id': + return P + elif self.iso_type == 'Hadamard': + return (P[0]+P[1],P[0]-P[1]) + else: + return (P[0]+self.i*P[1],P[0]-self.i*P[1]) + + + def evaluate(self,P): + FP=self.evaluate_isogeny(P) + FP=self.codomain_product.base_change_coords(self.N_split,FP) + + FP=product_to_theta_points_dim4(FP) + FP=TuplePoint([theta_point_to_montgomery_point(self.Ea,self.theta_null_Ea,self.iso_Ea(FP[0]),self.twist), + theta_point_to_montgomery_point(self.Ea,self.theta_null_Ea,self.iso_Ea(FP[1]),self.twist), + theta_point_to_montgomery_point(self.Ep,self.theta_null_Ep,FP[2]), + theta_point_to_montgomery_point(self.Ep,self.theta_null_Ep,FP[3])]) + + return FP + + def __call__(self,P): + return self.evaluate(P) + + + + + + + + + + + diff --git a/theta_lib/isogenies/Kani_gluing_isogeny_chain_dim4.py b/theta_lib/isogenies/Kani_gluing_isogeny_chain_dim4.py new file mode 100644 index 0000000..282219c --- /dev/null +++ b/theta_lib/isogenies/Kani_gluing_isogeny_chain_dim4.py @@ -0,0 +1,567 @@ +from sage.all import * +from ..utilities.discrete_log import weil_pairing_pari +from ..basis_change.canonical_basis_dim1 import make_canonical +from ..basis_change.kani_base_change import ( + fixed_deg_gluing_matrix_Phi1, + fixed_deg_gluing_matrix_Phi2, + fixed_deg_gluing_matrix_dim4, + clapoti_cob_matrix_dim2, + clapoti_cob_matrix_dim2_dim4, + gluing_base_change_matrix_dim2, + gluing_base_change_matrix_dim2_dim4, + gluing_base_change_matrix_dim2_F1, + gluing_base_change_matrix_dim2_F2, + kernel_basis, +) +from ..basis_change.base_change_dim2 import base_change_theta_dim2 +from ..basis_change.base_change_dim4 import base_change_theta_dim4 +from ..theta_structures.Theta_dim1 import ThetaStructureDim1 +from ..theta_structures.Theta_dim2 import ProductThetaStructureDim2 +from ..theta_structures.Tuple_point import TuplePoint +from ..theta_structures.Theta_dim4 import ProductThetaStructureDim2To4, ThetaPointDim4 +from ..isogenies_dim2.isogeny_chain_dim2 import IsogenyChainDim2 +from .gluing_isogeny_dim4 import GluingIsogenyDim4 + +class KaniFixedDegDim2Gluing: + def __init__(self,P_mp3,Q_mp3,a,b,c,d,u,f,m,strategy_dim2=None): + r""" + INPUT: + - P_mp3, Q_mp3: basis of E[2^(m+3)] such that pi(P_mp3)=P_mp3 and pi(Q_mp3)=-Q_mp3. + - a,b,c,d,u,f: integers such that a**2+c**2+p*(b**2+d**2)=u*(2**f-u), where p is + ths characteristic of the base field. + - m: integer such that m=min(v_2(a-b),v_2(a+b)). + + OUTPUT: Gluing isogeny chain F_{m+1}\circ...\circ F_1 containing the first m+1 steps of + the isogeny F: E^4 --> A*A' representing a u-isogeny in dimension 2. + """ + + P_mp2 = 2*P_mp3 + Q_mp2 = 2*Q_mp3 + P_4 = 2**m*P_mp2 + Q_4 = 2**m*Q_mp2 + + E = P_mp3.curve() + + # Canonical basis with S_4=(1,0) + _, _, R_4, S_4, M_dim1 = make_canonical(P_4,Q_4,4,preserve_pairing=True) + + Z4 = Integers(4) + M0 = matrix(Z4,[[M_dim1[0,0],0,M_dim1[0,1],0], + [0,M_dim1[0,0],0,M_dim1[0,1]], + [M_dim1[1,0],0,M_dim1[1,1],0], + [0,M_dim1[1,0],0,M_dim1[1,1]]]) + + # Theta structures + Theta_E = ThetaStructureDim1(E,R_4,S_4) + Theta_EE = ProductThetaStructureDim2(Theta_E,Theta_E) + + # Gluing change of basis in dimension 2 + M1 = fixed_deg_gluing_matrix_Phi1(u,a,b,c,d) + M2 = fixed_deg_gluing_matrix_Phi2(u,a,b,c,d) + + M10 = M0*M1 + M20 = M0*M2 + + Fp2 = E.base_field() + e4 = Fp2(weil_pairing_pari(R_4,S_4,4)) + + N_Phi1 = base_change_theta_dim2(M10,e4) + N_Phi2 = base_change_theta_dim2(M20,e4) + + # Gluing change of basis dimension 2 * dimension 2 --> dimension 4 + M_dim4 = fixed_deg_gluing_matrix_dim4(u,a,b,c,d,m) + + self.N_dim4 = base_change_theta_dim4(M_dim4,e4) + + # Kernel of Phi1 : E^2 --> A_m1 and Phi2 : E^2 --> A_m2 + two_mp2 = 2**(m+2) + two_mp3 = 2*two_mp2 + mu = inverse_mod(u,two_mp3) + + B_K_Phi1 = [TuplePoint((u%two_mp2)*P_mp2,((c+d)%two_mp2)*P_mp2), + TuplePoint((((d**2-c**2)*mu)%two_mp2)*Q_mp2,((c-d)%two_mp2)*Q_mp2)] + + B_K_Phi2 = [TuplePoint((u%two_mp2)*P_mp2,((d-c)%two_mp2)*P_mp2), + TuplePoint((((d**2-c**2)*mu)%two_mp2)*Q_mp2,(-(c+d)%two_mp2)*Q_mp2)] + + # Computation of the 2**m-isogenies Phi1 and Phi2 + self._Phi1=IsogenyChainDim2(B_K_Phi1,Theta_EE,N_Phi1,m,strategy_dim2) + self._Phi2=IsogenyChainDim2(B_K_Phi2,Theta_EE,N_Phi2,m,strategy_dim2) + + # Kernel of the (m+1)-th isogeny in dimension 4 F_{m+1}: A_m1*A_m2 --> B (gluing isogeny) + + B_K_dim4 =[TuplePoint((u%two_mp3)*P_mp3,E(0),((a+b)%two_mp3)*P_mp3,((c+d)%two_mp3)*P_mp3), + TuplePoint(E(0),(u%two_mp3)*P_mp3,((d-c)%two_mp3)*P_mp3,((a-b)%two_mp3)*P_mp3), + TuplePoint(((u-2**f)%two_mp3)*Q_mp3,E(0),((a-b)%two_mp3)*Q_mp3,((c-d)%two_mp3)*Q_mp3), + TuplePoint(E(0),((u-2**f)%two_mp3)*Q_mp3,((-c-d)%two_mp3)*Q_mp3,((a+b)%two_mp3)*Q_mp3)] + + L_K_dim4=B_K_dim4+[B_K_dim4[0]+B_K_dim4[1]] + + L_K_dim4=[[self._Phi1(TuplePoint(T[0],T[3])),self._Phi2(TuplePoint(T[1],T[2]))] for T in L_K_dim4] + + # Product Theta structure on A_m1*A_m2 + self.domain_product=ProductThetaStructureDim2To4(self._Phi1._codomain,self._Phi2._codomain) + + # Theta structure on A_m1*A_m2 after change of theta coordinates + self.domain_base_change=self.domain_product.base_change_struct(self.N_dim4) + + # Converting the kernel to the Theta structure domain_base_change + L_K_dim4=[self.domain_product.product_theta_point(T[0],T[1]) for T in L_K_dim4] + L_K_dim4=[self.domain_base_change.base_change_coords(self.N_dim4,T) for T in L_K_dim4] + + # Computing the gluing isogeny in dimension 4 + self._gluing_isogeny_dim4=GluingIsogenyDim4(self.domain_base_change,L_K_dim4,[(1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1),(1,1,0,0)]) + + # Translates for the evaluation of the gluing isogeny in dimension 4 + self.L_trans=[2*B_K_dim4[k] for k in range(2)] + self.L_trans_ind=[1,2] # Corresponds to multi indices (1,0,0,0) and (0,1,0,0) + + self._codomain=self._gluing_isogeny_dim4._codomain + + def evaluate(self,P): + if not isinstance(P, TuplePoint): + raise TypeError("KaniGluingIsogenyChainDim4 isogeny expects as input a TuplePoint on the domain product E^4") + + # Translating P + L_P_trans=[P+T for T in self.L_trans] + + # dim4 --> dim2 x dim2 + eval_P=[TuplePoint(P[0],P[3]),TuplePoint(P[1],P[2])] + eval_L_P_trans=[[TuplePoint(Q[0],Q[3]),TuplePoint(Q[1],Q[2])] for Q in L_P_trans] + + # evaluating through the dimension 2 isogenies + eval_P=[self._Phi1(eval_P[0]),self._Phi2(eval_P[1])] + eval_L_P_trans=[[self._Phi1(Q[0]),self._Phi2(Q[1])] for Q in eval_L_P_trans] + + # Product Theta structure and base change + eval_P=self.domain_product.product_theta_point(eval_P[0],eval_P[1]) + eval_P=self.domain_base_change.base_change_coords(self.N_dim4,eval_P) + + eval_L_P_trans=[self.domain_product.product_theta_point(Q[0],Q[1]) for Q in eval_L_P_trans] + eval_L_P_trans=[self.domain_base_change.base_change_coords(self.N_dim4,Q) for Q in eval_L_P_trans] + + return self._gluing_isogeny_dim4.special_image(eval_P,eval_L_P_trans,self.L_trans_ind) + + def __call__(self,P): + return self.evaluate(P) + + +class KaniClapotiGluing: + def __init__(self,points_mp3,points_mp2,points_4,integers,strategy_dim2=None,coerce=None): + self._coerce=coerce + Pu_mp3,Qu_mp3,Pv_mp3,Qv_mp3 = points_mp3 + Pu_mp2,Qu_mp2,Pv_mp2,Qv_mp2 = points_mp2 + Pu_4,Qu_4,Pv_4,Qv_4 = points_4 + gu,xu,yu,gv,xv,yv,Nbk,Nck,e,m = integers + + Eu=Pu_4.curve() + Ev=Pv_4.curve() + + lamb_u = inverse_mod(ZZ(gu),4) + lamb_v = inverse_mod(ZZ(gv*Nbk*Nck),4) + + + # 4-torsion canonical change of basis in Eu and Ev (Su=(1,0), Sv=(1,0)) + _,_,Ru,Su,Mu=make_canonical(Pu_4,lamb_u*Qu_4,4,preserve_pairing=True) + _,_,Rv,Sv,Mv=make_canonical(Pv_4,lamb_v*Qv_4,4,preserve_pairing=True) + + Z4 = Integers(4) + M0=matrix(Z4,[[Mu[0,0],0,Mu[1,0],0], + [0,Mv[0,0],0,Mv[1,0]], + [Mu[0,1],0,Mu[1,1],0], + [0,Mv[0,1],0,Mv[1,1]]]) + + self.M_product_dim2=M0 + + # Theta structures in dimension 1 and 2 + Theta_u=ThetaStructureDim1(Eu,Ru,Su) + Theta_v=ThetaStructureDim1(Ev,Rv,Sv) + + Theta_uv=ProductThetaStructureDim2(Theta_u,Theta_v) + + # Gluing change of basis in dimension 2 + M1 = clapoti_cob_matrix_dim2(integers) + M10 = M0*M1 + + Fp2 = Eu.base_field() + e4 = Fp2(weil_pairing_pari(Ru,Su,4)) + self.e4 = e4 + + N_dim2 = base_change_theta_dim2(M10,e4) + + # Gluing change of basis dimension 2 * dimension 2 --> dimension 4 + M2 = clapoti_cob_matrix_dim2_dim4(integers) + + self.N_dim4 = base_change_theta_dim4(M2,e4) + + # Kernel of the 2**m-isogeny chain in dimension 2 + two_mp2=2**(m+2) + two_mp3=2*two_mp2 + u=ZZ(gu*(xu**2+yu**2)) + mu=inverse_mod(u,two_mp2) + suv=ZZ(xu*xv+yu*yv) + duv=ZZ(xv*yu-xu*yv) + uNbk=(u*Nbk)%two_mp2 + gusuv=(gu*suv)%two_mp2 + xK2=(uNbk+gu*gv*mu*Nck*duv**2)%two_mp2 + B_K_dim2 = [TuplePoint(uNbk*Pu_mp2,gusuv*Pv_mp2),TuplePoint(xK2*Qu_mp2,gusuv*Qv_mp2)] + + # Computation of the 2**m-isogeny chain in dimension 2 + self._isogenies_dim2=IsogenyChainDim2(B_K_dim2,Theta_uv,N_dim2,m,strategy_dim2) + + # Kernel of the (m+1)-th isogeny in dimension 4 f_{m+1}: A_m^2 --> B (gluing isogeny) + xuNbk = (xu*Nbk)%two_mp3 + yuNbk = (yu*Nbk)%two_mp3 + inv_Nbk = inverse_mod(Nbk,two_mp3) + lambxu = ((1-2**e)*xu)%two_mp3 # extreme case m=e-2, 2^e = 2^(m+2) so 2^e/(u*Nbk) = 2^e mod 2^(m+3). + lambyu = ((1-2**e)*yu)%two_mp3 + xv_Nbk = (xv*inv_Nbk)%two_mp3 + yv_Nbk = (yv*inv_Nbk)%two_mp3 + + B_K_dim4 = [TuplePoint(xuNbk*Pu_mp3,yuNbk*Pu_mp3,xv*Pv_mp3,yv*Pv_mp3), + TuplePoint(-yuNbk*Pu_mp3,xuNbk*Pu_mp3,-yv*Pv_mp3,xv*Pv_mp3), + TuplePoint(lambxu*Qu_mp3,lambyu*Qu_mp3,xv_Nbk*Qv_mp3,yv_Nbk*Qv_mp3), + TuplePoint(-lambyu*Qu_mp3,lambxu*Qu_mp3,-yv_Nbk*Qv_mp3,xv_Nbk*Qv_mp3)] + + L_K_dim4=B_K_dim4+[B_K_dim4[0]+B_K_dim4[1]] + + L_K_dim4=[[self._isogenies_dim2(TuplePoint(L_K_dim4[k][0],L_K_dim4[k][2])),self._isogenies_dim2(TuplePoint(L_K_dim4[k][1],L_K_dim4[k][3]))] for k in range(5)] + + # Product Theta structure on A_m^2 + self.domain_product=ProductThetaStructureDim2To4(self._isogenies_dim2._codomain,self._isogenies_dim2._codomain) + + # Theta structure on A_m^2 after change of theta coordinates + self.domain_base_change=self.domain_product.base_change_struct(self.N_dim4) + + # Converting the kernel to the Theta structure domain_base_change + L_K_dim4=[self.domain_product.product_theta_point(L_K_dim4[k][0],L_K_dim4[k][1]) for k in range(5)] + L_K_dim4=[self.domain_base_change.base_change_coords(self.N_dim4,L_K_dim4[k]) for k in range(5)] + + # Computing the gluing isogeny in dimension 4 + self._gluing_isogeny_dim4=GluingIsogenyDim4(self.domain_base_change,L_K_dim4,[(1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1),(1,1,0,0)], coerce=self._coerce) + + # Translates for the evaluation of the gluing isogeny in dimension 4 + self.L_trans=[2*B_K_dim4[k] for k in range(2)] + self.L_trans_ind=[1,2] # Corresponds to multi indices (1,0,0,0) and (0,1,0,0) + + self._codomain=self._gluing_isogeny_dim4._codomain + + def evaluate(self,P): + if not isinstance(P, TuplePoint): + raise TypeError("KaniGluingIsogenyChainDim4 isogeny expects as input a TuplePoint on the domain product Eu^2 x Ev^2") + + # Translating P + L_P_trans=[P+T for T in self.L_trans] + + # dim4 --> dim2 x dim2 + eval_P=[TuplePoint(P[0],P[2]),TuplePoint(P[1],P[3])] + eval_L_P_trans=[[TuplePoint(Q[0],Q[2]),TuplePoint(Q[1],Q[3])] for Q in L_P_trans] + + # evaluating through the dimension 2 isogenies + eval_P=[self._isogenies_dim2(eval_P[0]),self._isogenies_dim2(eval_P[1])] + eval_L_P_trans=[[self._isogenies_dim2(Q[0]),self._isogenies_dim2(Q[1])] for Q in eval_L_P_trans] + + # Product Theta structure and base change + eval_P=self.domain_product.product_theta_point(eval_P[0],eval_P[1]) + eval_P=self.domain_base_change.base_change_coords(self.N_dim4,eval_P) + + eval_L_P_trans=[self.domain_product.product_theta_point(Q[0],Q[1]) for Q in eval_L_P_trans] + eval_L_P_trans=[self.domain_base_change.base_change_coords(self.N_dim4,Q) for Q in eval_L_P_trans] + + return self._gluing_isogeny_dim4.special_image(eval_P,eval_L_P_trans,self.L_trans_ind) + + def __call__(self,P): + return self.evaluate(P) + + + +class KaniGluingIsogenyChainDim4: + def __init__(self,points_m,points_4,a1,a2,q,m,strategy_dim2=None): + r""" + + INPUT: + - points_m: list of 4 points P1_m, Q1_m, R2_m, S2_m of order 2**(m+3) + such that (P1_m,Q1_m) generates E1[2**(m+3)] and (R2_m,S2_m) is + its image by sigma: E1 --> E2. + - points_4: list of 4 points P1_4, Q1_4, R2_4, S2_4 of order 4 obtained by + multiplying P1_m, Q1_m, R2_m, S2_m by 2**(m+1). + - a1, a2, q: integers such that a1**2+a2**2+q=2**e. + - m: 2-adic valuation of a2. + + OUTPUT: Composition of the m+1 first isogenies in the isogeny chained + E1^2*E1^2 --> E1^2*E2^2 parametrized by a1, a2 and sigma via Kani's lemma. + """ + + P1_m, Q1_m, R2_m, S2_m = points_m + P1_4, Q1_4, R2_4, S2_4 = points_4 + + E1=P1_m.curve() + E2=R2_m.curve() + + Fp2=E1.base_field() + + lamb=inverse_mod(q,4) + + _,_,T1,T2,MT=make_canonical(P1_4,Q1_4,4,preserve_pairing=True) + _,_,U1,U2,MU=make_canonical(R2_4,lamb*S2_4,4,preserve_pairing=True) + + Z4=Integers(4) + M0=matrix(Z4,[[MT[0,0],0,MT[1,0],0], + [0,MU[0,0],0,MU[1,0]], + [MT[0,1],0,MT[1,1],0], + [0,MU[0,1],0,MU[1,1]]]) + + self.M_product_dim2=M0 + + # Theta structures in dimension 1 and 2 + Theta1=ThetaStructureDim1(E1,T1,T2) + Theta2=ThetaStructureDim1(E2,U1,U2) + + Theta12=ProductThetaStructureDim2(Theta1,Theta2) + + self.Theta1=Theta1 + self.Theta2=Theta2 + self.Theta12=Theta12 + + # Gluing base change in dimension 2 + M1=gluing_base_change_matrix_dim2(a1,a2,q) + M10=M0*M1 + + self.M_gluing_dim2=M1 + + e4=Fp2(weil_pairing_pari(T1,T2,4)) + + self.e4=e4 + + N_dim2=base_change_theta_dim2(M10,e4) + #N_dim2=montgomery_to_theta_matrix_dim2(Theta12.zero().coords(),N1) + + # Gluing base change in dimension 4 + mua2=-M1[3,1] + M2=gluing_base_change_matrix_dim2_dim4(a1,a2,m,mua2) + + self.M_gluing_dim4=M2 + + self.N_dim4=base_change_theta_dim4(M2,e4) + + # Kernel of the 2**m-isogeny chain in dimension 2 + a1_red=a1%(2**(m+2)) + a2_red=a2%(2**(m+2)) + B_K_dim2=[TuplePoint(2*a1_red*P1_m-2*a2_red*Q1_m,2*R2_m),TuplePoint(2*a1_red*Q1_m+2*a2_red*P1_m,2*S2_m)] + + # Computation of the 2**m-isogeny chain in dimension 2 + self._isogenies_dim2=IsogenyChainDim2(B_K_dim2,Theta12,N_dim2,m,strategy_dim2) + + # Kernel of the (m+1)-th isogeny in dimension 4 f_{m+1}: A_m^2 --> B (gluing isogeny) + a1_red=a1%(2**(m+3)) + a2_red=a2%(2**(m+3)) + + a1P1_m=(a1_red)*P1_m + a2P1_m=(a2_red)*P1_m + a1Q1_m=(a1_red)*Q1_m + a2Q1_m=(a2_red)*Q1_m + + OE2=E2(0) + + B_K_dim4=[TuplePoint(a1P1_m,a2P1_m,R2_m,OE2),TuplePoint(a1Q1_m,a2Q1_m,S2_m,OE2), + TuplePoint(-a2P1_m,a1P1_m,OE2,R2_m),TuplePoint(-a2Q1_m,a1Q1_m,OE2,S2_m)] + L_K_dim4=B_K_dim4+[B_K_dim4[0]+B_K_dim4[1]] + + L_K_dim4=[[self._isogenies_dim2(TuplePoint(L_K_dim4[k][0],L_K_dim4[k][2])),self._isogenies_dim2(TuplePoint(L_K_dim4[k][1],L_K_dim4[k][3]))] for k in range(5)] + + # Product Theta structure on A_m^2 + self.domain_product=ProductThetaStructureDim2To4(self._isogenies_dim2._codomain,self._isogenies_dim2._codomain) + + # Theta structure on A_m^2 after base change + self.domain_base_change=self.domain_product.base_change_struct(self.N_dim4) + + # Converting the kernel to the Theta structure domain_base_change + L_K_dim4=[self.domain_product.product_theta_point(L_K_dim4[k][0],L_K_dim4[k][1]) for k in range(5)] + L_K_dim4=[self.domain_base_change.base_change_coords(self.N_dim4,L_K_dim4[k]) for k in range(5)] + + # Computing the gluing isogeny in dimension 4 + self._gluing_isogeny_dim4=GluingIsogenyDim4(self.domain_base_change,L_K_dim4,[(1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1),(1,1,0,0)]) + + # Translates for the evaluation of the gluing isogeny in dimension 4 + self.L_trans=[2*B_K_dim4[k] for k in range(2)] + self.L_trans_ind=[1,2] # Corresponds to multi indices (1,0,0,0) and (0,1,0,0) + + self._codomain=self._gluing_isogeny_dim4._codomain + + def evaluate(self,P): + if not isinstance(P, TuplePoint): + raise TypeError("KaniGluingIsogenyChainDim4 isogeny expects as input a TuplePoint on the domain product E1^2 x E2^2") + + # Translating P + L_P_trans=[P+T for T in self.L_trans] + + # dim4 --> dim2 x dim2 + eval_P=[TuplePoint(P[0],P[2]),TuplePoint(P[1],P[3])] + eval_L_P_trans=[[TuplePoint(Q[0],Q[2]),TuplePoint(Q[1],Q[3])] for Q in L_P_trans] + + # evaluating through the dimension 2 isogenies + eval_P=[self._isogenies_dim2(eval_P[0]),self._isogenies_dim2(eval_P[1])] + eval_L_P_trans=[[self._isogenies_dim2(Q[0]),self._isogenies_dim2(Q[1])] for Q in eval_L_P_trans] + + # Product Theta structure and base change + eval_P=self.domain_product.product_theta_point(eval_P[0],eval_P[1]) + eval_P=self.domain_base_change.base_change_coords(self.N_dim4,eval_P) + + eval_L_P_trans=[self.domain_product.product_theta_point(Q[0],Q[1]) for Q in eval_L_P_trans] + eval_L_P_trans=[self.domain_base_change.base_change_coords(self.N_dim4,Q) for Q in eval_L_P_trans] + + return self._gluing_isogeny_dim4.special_image(eval_P,eval_L_P_trans,self.L_trans_ind) + + def __call__(self,P): + return self.evaluate(P) + +class KaniGluingIsogenyChainDim4Half: + def __init__(self, points_m, a1, a2, q, m, Theta12, M_product_dim2, M_start_dim4, M_gluing_dim4, e4, dual=False,strategy_dim2=None):#points_m,points_4,a1,a2,q,m,precomputed_data=None,dual=False,strategy_dim2=None): + r""" + + INPUT: + - points_m: list of 4 points P1_m, Q1_m, R2_m, S2_m of order 2**(m+3) + such that (P1_m,Q1_m) generates E1[2**(m+3)] and (R2_m,S2_m) is + its image by sigma: E1 --> E2. + - points_4: list of 4 points P1_4, Q1_4, R2_4, S2_4 of order 4 obtained by + multiplying P1_m, Q1_m, R2_m, S2_m by 2**(m+1). + - a1, a2, q: integers such that a1**2+a2**2+q=2**e. + - m: 2-adic valuation of a2. + + OUTPUT: Composition of the m+1 first isogenies in the isogeny chained + E1^2*E1^2 --> E1^2*E2^2 parametrized by a1, a2 and sigma via Kani's lemma. + """ + + P1_m, Q1_m, R2_m, S2_m = points_m + + E1=P1_m.curve() + E2=R2_m.curve() + + Fp2=E1.base_field() + + self.M_product_dim2 = M_product_dim2 + + self.Theta12=Theta12 + + self.e4=e4 + + # Gluing base change in dimension 2 + if not dual: + M1=gluing_base_change_matrix_dim2_F1(a1,a2,q) + else: + M1=gluing_base_change_matrix_dim2_F2(a1,a2,q) + + M10=M_product_dim2*M1 + + self.M_gluing_dim2=M1 + + self.e4=e4 + + N_dim2=base_change_theta_dim2(M10,e4) + #N_dim2=montgomery_to_theta_matrix_dim2(Theta12.zero().coords(),N1) + + # Gluing base change in dimension 4 + + self.M_gluing_dim4 = M_gluing_dim4 + + self.N_dim4 = base_change_theta_dim4(M_gluing_dim4, e4) + + # Kernel of the 2**m-isogeny chain in dimension 2 + a1_red=a1%(2**(m+2)) + a2_red=a2%(2**(m+2)) + if not dual: + B_K_dim2=[TuplePoint(2*a1_red*P1_m-2*a2_red*Q1_m,2*R2_m),TuplePoint(2*a1_red*Q1_m+2*a2_red*P1_m,2*S2_m)] + else: + B_K_dim2=[TuplePoint(2*a1_red*P1_m+2*a2_red*Q1_m,-2*R2_m),TuplePoint(2*a1_red*Q1_m-2*a2_red*P1_m,-2*S2_m)] + + # Computation of the 2**m-isogeny chain in dimension 2 + self._isogenies_dim2=IsogenyChainDim2(B_K_dim2,Theta12,N_dim2,m,strategy_dim2) + + # Kernel of the (m+1)-th isogeny in dimension 4 f_{m+1}: A_m^2 --> B (gluing isogeny) + lamb=inverse_mod(q,2**(m+3)) + B_K_dim4=kernel_basis(M_start_dim4,m+1,P1_m,Q1_m,R2_m,lamb*S2_m) + L_K_dim4=B_K_dim4+[B_K_dim4[0]+B_K_dim4[1]] + + L_K_dim4=[[self._isogenies_dim2(TuplePoint(L_K_dim4[k][0],L_K_dim4[k][2])),self._isogenies_dim2(TuplePoint(L_K_dim4[k][1],L_K_dim4[k][3]))] for k in range(5)] + + # Product Theta structure on A_m^2 + self.domain_product=ProductThetaStructureDim2To4(self._isogenies_dim2._codomain,self._isogenies_dim2._codomain) + + # Theta structure on A_m^2 after base change + self.domain_base_change=self.domain_product.base_change_struct(self.N_dim4) + + # Converting the kernel to the Theta structure domain_base_change + L_K_dim4=[self.domain_product.product_theta_point(L_K_dim4[k][0],L_K_dim4[k][1]) for k in range(5)] + L_K_dim4=[self.domain_base_change.base_change_coords(self.N_dim4,L_K_dim4[k]) for k in range(5)] + + # Computing the gluing isogeny in dimension 4 + self._gluing_isogeny_dim4=GluingIsogenyDim4(self.domain_base_change,L_K_dim4,[(1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1),(1,1,0,0)]) + + # Translates for the evaluation of the gluing isogeny in dimension 4 + self.L_trans=[2*B_K_dim4[k] for k in range(2)] + self.L_trans_ind=[1,2] # Corresponds to multi indices (1,0,0,0) and (0,1,0,0) + + self._codomain=self._gluing_isogeny_dim4._codomain + + def evaluate(self,P): + if not isinstance(P, TuplePoint): + raise TypeError("KaniGluingIsogenyChainDim4 isogeny expects as input a TuplePoint on the domain product E1^2 x E2^2") + + # Translating P + L_P_trans=[P+T for T in self.L_trans] + + # dim4 --> dim2 x dim2 + eval_P=[TuplePoint(P[0],P[2]),TuplePoint(P[1],P[3])] + eval_L_P_trans=[[TuplePoint(Q[0],Q[2]),TuplePoint(Q[1],Q[3])] for Q in L_P_trans] + + # evaluating through the dimension 2 isogenies + eval_P=[self._isogenies_dim2(eval_P[0]),self._isogenies_dim2(eval_P[1])] + eval_L_P_trans=[[self._isogenies_dim2(Q[0]),self._isogenies_dim2(Q[1])] for Q in eval_L_P_trans] + + # Product Theta structure and base change + eval_P=self.domain_product.product_theta_point(eval_P[0],eval_P[1]) + eval_P=self.domain_base_change.base_change_coords(self.N_dim4,eval_P) + + eval_L_P_trans=[self.domain_product.product_theta_point(Q[0],Q[1]) for Q in eval_L_P_trans] + eval_L_P_trans=[self.domain_base_change.base_change_coords(self.N_dim4,Q) for Q in eval_L_P_trans] + + return self._gluing_isogeny_dim4.special_image(eval_P,eval_L_P_trans,self.L_trans_ind) + + def __call__(self,P): + return self.evaluate(P) + + def dual(self): + domain = self._codomain.hadamard() + codomain_base_change = self.domain_base_change + codomain_product = self.domain_product + N_dim4 = self.N_dim4.inverse() + isogenies_dim2 = self._isogenies_dim2.dual() + splitting_isogeny_dim4 = self._gluing_isogeny_dim4.dual() + + return KaniSplittingIsogenyChainDim4(domain, codomain_base_change, codomain_product, N_dim4, isogenies_dim2, splitting_isogeny_dim4) + +class KaniSplittingIsogenyChainDim4: + def __init__(self, domain, codomain_base_change, codomain_product, N_dim4, isogenies_dim2, splitting_isogeny_dim4): + self._domain = domain + self.codomain_base_change = codomain_base_change + self.codomain_product = codomain_product + self.N_dim4 = N_dim4 + self._isogenies_dim2 = isogenies_dim2 + self._splitting_isogeny_dim4 = splitting_isogeny_dim4 + + def evaluate(self,P): + if not isinstance(P, ThetaPointDim4): + raise TypeError("KaniSplittingIsogenyChainDim4 isogeny expects as input a ThetaPointDim4") + + Q = self._splitting_isogeny_dim4(P) + Q = self.codomain_product.base_change_coords(self.N_dim4, Q) + Q1, Q2 = self.codomain_product.to_theta_points(Q) + Q1, Q2 = self._isogenies_dim2._domain(Q1.hadamard()), self._isogenies_dim2._domain(Q2.hadamard()) + + Q1 = self._isogenies_dim2(Q1) + Q2 = self._isogenies_dim2(Q2) + + return TuplePoint(Q1[0],Q2[0],Q1[1],Q2[1]) + + def __call__(self,P): + return self.evaluate(P) diff --git a/theta_lib/isogenies/gluing_isogeny_dim4.py b/theta_lib/isogenies/gluing_isogeny_dim4.py new file mode 100644 index 0000000..cd738c9 --- /dev/null +++ b/theta_lib/isogenies/gluing_isogeny_dim4.py @@ -0,0 +1,200 @@ +from sage.all import * + +from ..theta_structures.Theta_dim4 import ThetaStructureDim4 +from ..theta_structures.theta_helpers_dim4 import hadamard, squared, batch_inversion, multindex_to_index +from .tree import Tree +from .isogeny_dim4 import IsogenyDim4, DualIsogenyDim4 + +def proj_equal(P1, P2): + if len(P1) != len(P2): + return False + for i in range(0, len(P1)): + if P1[i]==0: + if P2[i] != 0: + return False + else: + break + r=P1[i] + s=P2[i] + for i in range(0, len(P1)): + if P1[i]*s != P2[i]*r: + return False + return True + +class GluingIsogenyDim4(IsogenyDim4): + def __init__(self,domain,L_K_8,L_K_8_ind, coerce=None): + r""" + Input: + - domain: a ThetaStructureDim4. + - L_K_8: list of points of 8-torsion in the kernel. + - L_K_8_ind: list of corresponding multindices (i0,i1,i2,i3) of points in L_K_8 + (L_K_8[i]=i0*P0+i1*P1+i2*P2+i3*P3, where (P0,..,P3) is a basis of K_8 (compatible with + the canonical basis of K_2) and L_K_8_ind[i]=(i0,i1,i2,i3)). + """ + + if not isinstance(domain, ThetaStructureDim4): + raise ValueError("Argument domain should be a ThetaStructureDim4 object.") + self._domain = domain + self._precomputation=None + self._coerce=coerce + self._special_compute_codomain(L_K_8,L_K_8_ind) + + #a_i2=squared(self._domain.zero()) + #HB_i2=hadamard(squared(hadamard(self._codomain.zero()))) + #for i in range(16): + #print(HB_i2[i]/a_i2[i]) + + def _special_compute_codomain(self,L_K_8,L_K_8_ind): + r""" + Input: + - L_K_8: list of points of 8-torsion in the kernel. + - L_K_8_ind: list of corresponding multindices (i0,i1,i2,i3) of points in L_K_8 + (L_K_8[i]=i0*P0+i1*P1+i2*P2+i3*P3, where (P0,..,P3) is a basis of K_8 (compatible with + the canonical basis of K_2) and L_K_8_ind[i]=(i0,i1,i2,i3)). + + Output: + - codomain of the isogeny. + Also initializes self._precomputation, containing the inverse of theta-constants. + """ + HSK_8=[hadamard(squared(P.coords())) for P in L_K_8] + + # Choice of reference index j_0<->chi_0 corresponding to a non-vanishing theta-constant. + found_tree=False + j_0=0 + while not found_tree: + found_k0=False + for k in range(len(L_K_8)): + if HSK_8[k][j_0]!=0: + k_0=k + found_k0=True + break + if not found_k0: + j_0+=1 + else: + j0pk0=j_0^multindex_to_index(L_K_8_ind[k_0]) + # List of tuples of indices (index chi of the denominator: HS(f(P_k))_chi, + #index chi.chi_k of the numerator: HS(f(P_k))_chi.chi_k, index k). + L_ratios_ind=[(j_0,j0pk0,k_0)] + L_covered_ind=[j_0,j0pk0] + + # Tree containing the the theta-null points indices as nodes and the L_ratios_ind reference indices as edges. + tree_ratios=Tree(j_0) + tree_ratios.add_child(Tree(j0pk0),0) + + # Filling in the tree + tree_filled=False + while not tree_filled: + found_j=False + for j in L_covered_ind: + for k in range(len(L_K_8)): + jpk=j^multindex_to_index(L_K_8_ind[k]) + if jpk not in L_covered_ind and HSK_8[k][j]!=0: + L_covered_ind.append(jpk) + L_ratios_ind.append((j,jpk,k)) + tree_j=tree_ratios.look_node(j) + tree_j.add_child(Tree(jpk),len(L_ratios_ind)-1) + found_j=True + #break + #if found_j: + #break + if not found_j or len(L_covered_ind)==16: + tree_filled=True + if len(L_covered_ind)!=16: + j_0+=1 + else: + found_tree=True + + L_denom=[HSK_8[t[2]][t[0]] for t in L_ratios_ind] + L_denom_inv=batch_inversion(L_denom) + L_num=[HSK_8[t[2]][t[1]] for t in L_ratios_ind] + L_ratios=[L_num[i]*L_denom_inv[i] for i in range(15)] + + L_coords_ind=tree_ratios.edge_product(L_ratios) + + O_coords=[ZZ(0) for i in range(16)] + for t in L_coords_ind: + if self._coerce: + O_coords[t[1]]=self._coerce(t[0]) + else: + O_coords[t[1]]=t[0] + + # Precomputation + # TODO: optimize inversions and give precomputation to the codomain _arithmetic_precomputation + L_prec=[] + L_prec_ind=[] + for i in range(16): + if O_coords[i]!=0: + L_prec.append(O_coords[i]) + L_prec_ind.append(i) + L_prec_inv=batch_inversion(L_prec) + precomputation=[None for i in range(16)] + for i in range(len(L_prec)): + precomputation[L_prec_ind[i]]=L_prec_inv[i] + + self._precomputation=precomputation + + for k in range(len(L_K_8)): + for j in range(16): + jpk=j^multindex_to_index(L_K_8_ind[k]) + assert HSK_8[k][j]*O_coords[jpk]==HSK_8[k][jpk]*O_coords[j] + + assert proj_equal(squared(self._domain._null_point.coords()), hadamard(squared(O_coords))) + + self._codomain=ThetaStructureDim4(hadamard(O_coords),null_point_dual=O_coords) + + def special_image(self,P,L_trans,L_trans_ind): + r"""Used when we cannot evaluate the isogeny self because the codomain has zero + dual theta constants. + + Input: + - P: ThetaPointDim4 of the domain. + - L_trans: list of translates of P+T of P by points of 4-torsion T above the kernel. + - L_trans_ind: list of indices of the translation 4-torsion points T. + If L_trans[i]=\sum i_j*B_K4[j] then L_trans_ind[j]=\sum 2**j*i_j. + + Output: + - the image of P by the isogeny self. + """ + HS_P=hadamard(squared(P.coords())) + HSL_trans=[hadamard(squared(Q.coords())) for Q in L_trans] + O_coords=self._codomain.null_point_dual() + + # L_lambda_inv: List of inverses of lambda_i such that: + # HS(P+Ti)=(lambda_i*U_{chi.chi_i,0}(f(P))*U_{chi,0}(0))_chi. + L_lambda_inv_num=[] + L_lambda_inv_denom=[] + + for k in range(len(L_trans)): + for j in range(16): + jpk=j^L_trans_ind[k] + if HSL_trans[k][j]!=0 and O_coords[jpk]!=0: + L_lambda_inv_num.append(HS_P[jpk]*O_coords[j]) + L_lambda_inv_denom.append(HSL_trans[k][j]*O_coords[jpk]) + break + L_lambda_inv_denom=batch_inversion(L_lambda_inv_denom) + L_lambda_inv=[L_lambda_inv_num[i]*L_lambda_inv_denom[i] for i in range(len(L_trans))] + + for k in range(len(L_trans)): + for j in range(16): + jpk=j^L_trans_ind[k] + assert HS_P[jpk]*O_coords[j]==L_lambda_inv[k]*HSL_trans[k][j]*O_coords[jpk] + + U_fP=[] + for i in range(16): + if self._precomputation[i]!=None: + U_fP.append(self._precomputation[i]*HS_P[i]) + else: + for k in range(len(L_trans)): + ipk=i^L_trans_ind[k] + if self._precomputation[ipk]!=None: + U_fP.append(self._precomputation[ipk]*HSL_trans[k][ipk]*L_lambda_inv[k]) + break + + fP=hadamard(U_fP) + if self._coerce: + fP=[self._coerce(x) for x in fP] + + return self._codomain(fP) + + def dual(self): + return DualIsogenyDim4(self._codomain,self._domain, hadamard=False) diff --git a/theta_lib/isogenies/isogeny_chain_dim4.py b/theta_lib/isogenies/isogeny_chain_dim4.py new file mode 100644 index 0000000..8c0b78e --- /dev/null +++ b/theta_lib/isogenies/isogeny_chain_dim4.py @@ -0,0 +1,114 @@ +from sage.all import * +from ..utilities.strategy import precompute_strategy_with_first_eval, precompute_strategy_with_first_eval_and_splitting +from .isogeny_dim4 import IsogenyDim4 + + +class IsogenyChainDim4: + def __init__(self, B_K, first_isogenies, e, m, splitting=True, strategy = None): + self.e=e + self.m=m + + if strategy == None: + strategy = self.get_strategy(splitting) + self.strategy = strategy + + self._isogenies=self.isogeny_chain(B_K, first_isogenies) + + + def get_strategy(self,splitting): + if splitting: + return precompute_strategy_with_first_eval_and_splitting(self.e,self.m,M=1,S=0.8,I=100) + else: + return precompute_strategy_with_first_eval(self.e,self.m,M=1,S=0.8,I=100) + + def isogeny_chain(self, B_K, first_isogenies): + """ + Compute the isogeny chain and store intermediate isogenies for evaluation + """ + # Store chain of (2,2)-isogenies + isogeny_chain = [] + + # Bookkeeping for optimal strategy + strat_idx = 0 + level = [0] + ker = B_K + kernel_elements = [ker] + + # Length of the chain + n=self.e-self.m + + for k in range(n): + prev = sum(level) + ker = kernel_elements[-1] + + while prev != (n - 1 - k): + level.append(self.strategy[strat_idx]) + prev += self.strategy[strat_idx] + + # Perform the doublings and update kernel elements + # Prevent the last unnecessary doublings for first isogeny computation + if k>0 or prev!=n-1: + ker = [ker[i].double_iter(self.strategy[strat_idx]) for i in range(4)] + kernel_elements.append(ker) + + # Update bookkeeping variable + strat_idx += 1 + + # Compute the codomain from the 8-torsion + if k==0: + phi = first_isogenies + else: + phi = IsogenyDim4(Th,ker) + + # Update the chain of isogenies + Th = phi._codomain + # print(parent(Th.null_point().coords()[0])) + isogeny_chain.append(phi) + + # Remove elements from list + if k>0: + kernel_elements.pop() + level.pop() + + # Push through points for the next step + kernel_elements = [[phi(T) for T in kernel] for kernel in kernel_elements] + # print([[parent(T.coords()[0]) for T in kernel] for kernel in kernel_elements]) + + return isogeny_chain + + def evaluate_isogeny(self,P): + Q=P + for f in self._isogenies: + Q=f(Q) + return Q + + def __call__(self,P): + return self.evaluate_isogeny(P) + + def dual(self): + n=len(self._isogenies) + isogenies=[] + for i in range(n): + isogenies.append(self._isogenies[n-1-i].dual()) + return DualIsogenyChainDim4(isogenies) + + +class DualIsogenyChainDim4: + def __init__(self,isogenies): + self._isogenies=isogenies + + def evaluate_isogeny(self,P): + n=len(self._isogenies) + Q=P + for j in range(n): + Q=self._isogenies[j](Q) + return Q + + def __call__(self,P): + return self.evaluate_isogeny(P) + + + + + + diff --git a/theta_lib/isogenies/isogeny_dim4.py b/theta_lib/isogenies/isogeny_dim4.py new file mode 100644 index 0000000..2f483bf --- /dev/null +++ b/theta_lib/isogenies/isogeny_dim4.py @@ -0,0 +1,162 @@ +from sage.all import * + +from ..theta_structures.Theta_dim4 import ThetaStructureDim4 +from ..theta_structures.theta_helpers_dim4 import hadamard, squared, batch_inversion +from .tree import Tree + +class IsogenyDim4: + def __init__(self,domain,K_8,codomain=None,precomputation=None): + r""" + Input: + - domain: a ThetaStructureDim4. + - K_8: a list of 4 points of 8-torision (such that 4*K_8 is a kernel basis), used to compute the codomain. + - codomain: a ThetaStructureDim4 (for the codomain, used only when K_8 is None). + - precomputation: list of inverse of dual theta constants of the codomain, used to compute the image. + """ + + if not isinstance(domain, ThetaStructureDim4): + raise ValueError("Argument domain should be a ThetaStructureDim4 object.") + self._domain = domain + self._precomputation=None + if K_8!=None: + self._compute_codomain(K_8) + else: + self._codomain=codomain + self._precomputation=precomputation + + def _compute_codomain(self,K_8): + r""" + Input: + - K_8: a list of 4 points of 8-torision (such that 4*K_8 is a kernel basis). + + Output: + - codomain of the isogeny. + Also initializes self._precomputation, containing the inverse of theta-constants. + """ + HSK_8=[hadamard(squared(P.coords())) for P in K_8] + + # Choice of reference index j_0<->chi_0 corresponding to a non-vanishing theta-constant. + found_tree=False + j_0=0 + while not found_tree: + found_k0=False + for k in range(4): + if j_0>15: + raise NotImplementedError("The codomain of this 2-isogeny could not be computed.\nWe may have encountered a product of abelian varieties\nsomewhere unexpected along the chain.\nThis is exceptionnal and should not happen in larger characteristic.") + if HSK_8[k][j_0]!=0: + k_0=k + found_k0=True + break + if not found_k0: + j_0+=1 + else: + j0pk0=j_0^(2**k_0) + # List of tuples of indices (index chi of the denominator: HS(f(P_k))_chi, + #index chi.chi_k of the numerator: HS(f(P_k))_chi.chi_k, index k). + L_ratios_ind=[(j_0,j0pk0,k_0)] + L_covered_ind=[j_0,j0pk0] + + # Tree containing the the theta-null points indices as nodes and the L_ratios_ind reference indices as edges. + tree_ratios=Tree(j_0) + tree_ratios.add_child(Tree(j0pk0),k_0) + + # Filling in the tree + tree_filled=False + while not tree_filled: + found_j=False + for j in L_covered_ind: + for k in range(4): + jpk=j^(2**k) + if jpk not in L_covered_ind and HSK_8[k][j]!=0: + L_covered_ind.append(jpk) + L_ratios_ind.append((j,jpk,k)) + tree_j=tree_ratios.look_node(j) + tree_j.add_child(Tree(jpk),len(L_ratios_ind)-1) + found_j=True + break + if found_j: + break + if not found_j or len(L_covered_ind)==16: + tree_filled=True + if len(L_covered_ind)!=16: + j_0+=1 + else: + found_tree=True + + L_denom=[HSK_8[t[2]][t[0]] for t in L_ratios_ind] + L_denom_inv=batch_inversion(L_denom) + L_num=[HSK_8[t[2]][t[1]] for t in L_ratios_ind] + L_ratios=[L_num[i]*L_denom_inv[i] for i in range(15)] + + L_coords_ind=tree_ratios.edge_product(L_ratios) + + O_coords=[ZZ(0) for i in range(16)] + for t in L_coords_ind: + O_coords[t[1]]=t[0] + + # Precomputation + # TODO: optimize inversions + L_prec=[] + L_prec_ind=[] + for i in range(16): + if O_coords[i]!=0: + L_prec.append(O_coords[i]) + L_prec_ind.append(i) + L_prec_inv=batch_inversion(L_prec) + precomputation=[None for i in range(16)] + for i in range(len(L_prec)): + precomputation[L_prec_ind[i]]=L_prec_inv[i] + + self._precomputation=precomputation + # Assumes there is no zero theta constant. Otherwise, squared(precomputation) will raise an error (None**2 does not exist) + self._codomain=ThetaStructureDim4(hadamard(O_coords),null_point_dual=O_coords) + + def codomain(self): + return self._codomain + + def domain(self): + return self._domain + + def image(self,P): + HS_P=list(hadamard(squared(P.coords()))) + + for i in range(16): + HS_P[i] *=self._precomputation[i] + + return self._codomain(hadamard(HS_P)) + + def dual(self): + return DualIsogenyDim4(self._codomain,self._domain, hadamard=True) + + def __call__(self,P): + return self.image(P) + + +class DualIsogenyDim4: + def __init__(self,domain,codomain,hadamard=True): + # domain and codomain are respectively the domain and codomain of \tilde{f}: domain-->codomain, + # so respectively the codomain and domain of f: codomain-->domain. + # By convention, domain input is given in usual coordinates (ker(\tilde{f})=K_2). + # codomain is in usual coordinates if hadamard, in dual coordinates otherwise. + self._domain=domain.hadamard() + self._hadamard=hadamard + if hadamard: + self._codomain=codomain.hadamard() + self._precomputation=batch_inversion(codomain.zero().coords()) + else: + self._codomain=codomain + self._precomputation=batch_inversion(codomain.zero().coords()) + + def image(self,P): + # When ker(f)=K_2, ker(\tilde{f})=K_1 so ker(\tilde{f})=K_2 after hadamard transformation of the + # new domain (ex codomain) + HS_P=list(hadamard(squared(P.coords()))) + for i in range(16): + HS_P[i] *=self._precomputation[i] + if self._hadamard: + return self._codomain(hadamard(HS_P)) + else: + return self._codomain(HS_P) + + def __call__(self,P): + return self.image(P) diff --git a/theta_lib/isogenies/tree.py b/theta_lib/isogenies/tree.py new file mode 100644 index 0000000..a6e3da3 --- /dev/null +++ b/theta_lib/isogenies/tree.py @@ -0,0 +1,28 @@ +from sage.all import * + +class Tree: + def __init__(self,node): + self._node=node + self._edges=[] + self._children=[] + + def add_child(self,child,edge): + self._children.append(child) + self._edges.append(edge) + + def look_node(self,node): + if self._node==node: + return self + elif len(self._children)>0: + for child in self._children: + t_node=child.look_node(node) + if t_node!=None: + return t_node + + def edge_product(self,L_factors,factor_node=ZZ(1)): + n=len(self._children) + L_prod=[(factor_node,self._node)] + for i in range(n): + L_prod+=self._children[i].edge_product(L_factors,factor_node*L_factors[self._edges[i]]) + return L_prod + diff --git a/theta_lib/isogenies_dim2/gluing_isogeny_dim2.py b/theta_lib/isogenies_dim2/gluing_isogeny_dim2.py new file mode 100644 index 0000000..2dac387 --- /dev/null +++ b/theta_lib/isogenies_dim2/gluing_isogeny_dim2.py @@ -0,0 +1,292 @@ +""" +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 new file mode 100644 index 0000000..23adb23 --- /dev/null +++ b/theta_lib/isogenies_dim2/isogeny_chain_dim2.py @@ -0,0 +1,178 @@ +""" +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 new file mode 100644 index 0000000..b740ab1 --- /dev/null +++ b/theta_lib/isogenies_dim2/isogeny_dim2.py @@ -0,0 +1,229 @@ +""" +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 new file mode 100644 index 0000000..e8e94dd --- /dev/null +++ b/theta_lib/theta_structures/Theta_dim1.py @@ -0,0 +1,98 @@ +from sage.all import * +from ..utilities.supersingular import compute_linearly_independent_point +from .montgomery_theta import ( + montgomery_point_to_theta_point, + theta_point_to_montgomery_point, + torsion_to_theta_null_point, +) + + +class ThetaStructureDim1: + def __init__(self,E,P=None,Q=None): + self.E=E + + a_inv=E.a_invariants() + + A =a_inv[1] + if a_inv != (0,A,0,1,0): + raise ValueError("The elliptic curve E is not in the Montgomery model.") + + if Q==None: + y2=A-2 + y=y2.sqrt() + Q=E([-1,y,1]) + P=compute_linearly_independent_point(E,Q,4) + else: + if Q[0]!=-1: + raise ValueError("You should enter a canonical 4-torsion basis. Q[0] should be -1.") + + self.P=P + self.Q=Q + + self._base_ring=E.base_ring() + + self._point=ThetaPointDim1 + self._null_point=self._point(self,torsion_to_theta_null_point(P)) + + def null_point(self): + """ + """ + return self._null_point + + def base_ring(self): + """ + """ + return self._base_ring + + def zero(self): + """ + """ + return self.null_point() + + def elliptic_curve(self): + return self.E + + def torsion_basis(self): + return (self.P,self.Q) + + def __call__(self,coords): + r""" + Input: either a tuple or list of 2 coordinates or an elliptic curve point. + + Output: the corresponding theta point for the self theta structure. + """ + if isinstance(coords,tuple): + return self._point(self,coords) + elif isinstance(coords,list): + return self._point(self,coords) + else: + return self._point(self,montgomery_point_to_theta_point(self.null_point().coords(),coords)) + + def __repr__(self): + return f"Theta structure on {self.elliptic_curve()} with null point: {self.null_point()} induced by the 4-torsion basis {self.torsion_basis()}" + + def to_montgomery_point(self,P): + return theta_point_to_montgomery_point(self.E,self.null_point().coords(),P.coords()) + + +class ThetaPointDim1: + def __init__(self, parent, coords): + """ + """ + if not isinstance(parent, ThetaStructureDim1): + raise ValueError("Entry parent should be a ThetaStructureDim1 object.") + + self._parent = parent + self._coords = tuple(coords) + + + def coords(self): + return self._coords + + def __repr__(self): + return f"Theta point with coordinates: {self.coords()}" + + + + + diff --git a/theta_lib/theta_structures/Theta_dim2.py b/theta_lib/theta_structures/Theta_dim2.py new file mode 100644 index 0000000..0d9955a --- /dev/null +++ b/theta_lib/theta_structures/Theta_dim2.py @@ -0,0 +1,440 @@ +# Sage Imports +from sage.all import Integer +from sage.structure.element import get_coercion_model, RingElement + +cm = get_coercion_model() + +from .theta_helpers_dim2 import batch_inversion, product_theta_point +from .Theta_dim1 import ThetaStructureDim1 +from .Tuple_point import TuplePoint +from ..basis_change.base_change_dim2 import apply_base_change_theta_dim2 + +# =========================================== # +# Class for Theta Structure (level-2) # +# =========================================== # + + +class ThetaStructureDim2: + """ + Class for the Theta Structure in dimension 2, defined by its theta null point. This type + represents the generic domain/codomain of the (2,2)-isogeny in the theta model. + """ + + def __init__(self, null_point, null_point_dual=None): + if not len(null_point) == 4: + raise ValueError + + self._base_ring = cm.common_parent(*(c.parent() for c in null_point)) + self._point = ThetaPointDim2 + self._precomputation = None + + self._null_point = self._point(self, null_point) + self._null_point_dual = null_point_dual + + def null_point(self): + """ + Return the null point of the given theta structure + """ + return self._null_point + + def null_point_dual(self): + if self._null_point_dual==None: + self._null_point_dual = self._point.to_hadamard(*self.coords()) + return self._null_point_dual + + def base_ring(self): + """ + Return the base ring of the common parent of the coordinates of the null point + """ + return self._base_ring + + def zero(self): + """ + The additive identity is the theta null point + """ + return self.null_point() + + def zero_dual(self): + return self.null_point_dual() + + def __repr__(self): + return f"Theta structure over {self.base_ring()} with null point: {self.null_point()}" + + def coords(self): + """ + Return the coordinates of the theta null point of the theta structure + """ + return self.null_point().coords() + + def hadamard(self): + """ + Compute the Hadamard transformation of the theta structure + """ + return ThetaStructureDim2(self.null_point_dual(),null_point_dual=self.coords()) + + def squared_theta(self): + """ + Square the coefficients and then compute the Hadamard transformation of + the theta null point of the theta structure + """ + return self.null_point().squared_theta() + + def _arithmetic_precomputation(self): + """ + Precompute 6 field elements used in arithmetic and isogeny computations + """ + if self._precomputation is None: + a, b, c, d = self.null_point().coords() + + # Technically this computes 4A^2, 4B^2, ... + # but as we take quotients this doesnt matter + # Cost: 4S + AA, BB, CC, DD = self.squared_theta() + + # Precomputed constants for addition and doubling + b_inv, c_inv, d_inv, BB_inv, CC_inv, DD_inv = batch_inversion([ + b, c, d, BB, CC, DD] + ) + + y0 = a * b_inv + z0 = a * c_inv + t0 = a * d_inv + + Y0 = AA * BB_inv + Z0 = AA * CC_inv + T0 = AA * DD_inv + + self._precomputation = (y0, z0, t0, Y0, Z0, T0) + return self._precomputation + + def __call__(self, coords): + return self._point(self, coords) + + def base_change_struct(self,N): + null_coords=self.null_point().coords() + new_null_coords=apply_base_change_theta_dim2(N,null_coords) + return ThetaStructure(new_null_coords) + + def base_change_coords(self,N,P): + coords=P.coords() + new_coords=apply_base_change_theta_dim2(N,coords) + return self.__call__(new_coords) + + +# =================================================== # +# Class for Product Theta Structure (level-2) # +# =================================================== # + + +class ProductThetaStructureDim2(ThetaStructureDim2): + def __init__(self,*args): + r"""Defines the product theta structure at level 2 of 2 elliptic curves. + + Input: Either + - 2 theta structures of dimension 1: T0, T1; + - 2 elliptic curves: E0, E1. + - 2 elliptic curves E0, E1 and their respective canonical 4-torsion basis B0, B1. + """ + if len(args)==2: + theta_structures=list(args) + for k in range(2): + if not isinstance(theta_structures[k],ThetaStructureDim1): + theta_structures[k]=ThetaStructureDim1(theta_structures[k]) + elif len(args)==4: + theta_structures=[ThetaStructureDim1(args[k],args[2+k][0],args[2+k][1]) for k in range(2)] + else: + raise ValueError("2 or 4 arguments expected but {} were given.\nYou should enter a list of 2 elliptic curves or ThetaStructureDim1\nor a list of 2 elliptic curves with a 4-torsion basis for each of them.".format(len(args))) + + self._theta_structures=theta_structures + + null_point=product_theta_point(theta_structures[0].zero().coords(),theta_structures[1].zero().coords()) + + ThetaStructureDim2.__init__(self,null_point) + + def product_theta_point(self,theta_points): + t0,t1=theta_points[0].coords() + u0,u1=theta_points[1].coords() + return self._point(self,[t0*u0,t1*u0,t0*u1,t1*u1]) + + def __call__(self,point): + if isinstance(point,TuplePoint): + theta_points=[] + theta_structures=self._theta_structures + for i in range(2): + theta_points.append(theta_structures[i](point[i])) + return self.product_theta_point(theta_points) + else: + return self._point(self,point) + + def to_theta_points(self,P): + coords=P.coords() + theta_coords=[(coords[0],coords[1]),(coords[1],coords[3])] + theta_points=[self._theta_structures[i](theta_coords[i]) for i in range(2)] + return theta_points + + def to_tuple_point(self,P): + theta_points=self.to_theta_points(P) + montgomery_points=[self._theta_structures[i].to_montgomery_point(theta_points[i]) for i in range(2)] + return TuplePoint(montgomery_points) + + +# ======================================= # +# Class for Theta Point (level-2) # +# ======================================= # + + +class ThetaPointDim2: + """ + A Theta Point in the level-2 Theta Structure is defined with four projective + coordinates + + We cannot perform arbitrary arithmetic, but we can compute doubles and + differential addition, which like x-only points on the Kummer line, allows + for scalar multiplication + """ + + def __init__(self, parent, coords): + if not isinstance(parent, ThetaStructureDim2) and not isinstance(parent, ProductThetaStructureDim2): + raise ValueError + + self._parent = parent + self._coords = tuple(coords) + + self._hadamard = None + self._squared_theta = None + + def parent(self): + """ + Return the parent of the element, of type ThetaStructureDim2 + """ + return self._parent + + def theta(self): + """ + Return the parent theta structure of this ThetaPointDim2""" + return self.parent() + + def coords(self): + """ + Return the projective coordinates of the ThetaPointDim2 + """ + return self._coords + + def is_zero(self): + """ + An element is zero if it is equivalent to the null point of the parent + ThetaStrcuture + """ + return self == self.parent().zero() + + @staticmethod + def to_hadamard(x_00, x_10, x_01, x_11): + """ + Compute the Hadamard transformation of four coordinates, using recursive + formula. + """ + x_00, x_10 = (x_00 + x_10, x_00 - x_10) + x_01, x_11 = (x_01 + x_11, x_01 - x_11) + return x_00 + x_01, x_10 + x_11, x_00 - x_01, x_10 - x_11 + + def hadamard(self): + """ + Compute the Hadamard transformation of this element + """ + if self._hadamard is None: + self._hadamard = self.to_hadamard(*self.coords()) + return self._hadamard + + @staticmethod + def to_squared_theta(x, y, z, t): + """ + Square the coordinates and then compute the Hadamard transform of the + input + """ + return ThetaPointDim2.to_hadamard(x * x, y * y, z * z, t * t) + + def squared_theta(self): + """ + Compute the Squared Theta transformation of this element + which is the square operator followed by Hadamard. + """ + if self._squared_theta is None: + self._squared_theta = self.to_squared_theta(*self.coords()) + return self._squared_theta + + def double(self): + """ + Computes [2]*self + + NOTE: Assumes that no coordinate is zero + + Cost: 8S 6M + """ + # If a,b,c,d = 0, then the codomain of A->A/K_2 is a product of + # elliptic curves with a non product theta structure. + # Unless we are very unlucky, A/K_1 will not be in this case, so we + # just need to Hadamard, double, and Hadamard inverse + # If A,B,C,D=0 then the domain itself is a product of elliptic + # curves with a non product theta structure. The Hadamard transform + # will not change this, we need a symplectic change of variable + # that puts us back in a product theta structure + y0, z0, t0, Y0, Z0, T0 = self.parent()._arithmetic_precomputation() + + # Temp coordinates + # Cost 8S 3M + xp, yp, zp, tp = self.squared_theta() + xp = xp**2 + yp = Y0 * yp**2 + zp = Z0 * zp**2 + tp = T0 * tp**2 + + # Final coordinates + # Cost 3M + X, Y, Z, T = self.to_hadamard(xp, yp, zp, tp) + X = X + Y = y0 * Y + Z = z0 * Z + T = t0 * T + + coords = (X, Y, Z, T) + return self._parent(coords) + + def diff_addition(P, Q, PQ): + """ + Given the theta points of P, Q and P-Q computes the theta point of + P + Q. + + NOTE: Assumes that no coordinate is zero + + Cost: 8S 17M + """ + # Extract out the precomputations + Y0, Z0, T0 = P.parent()._arithmetic_precomputation()[-3:] + + # Transform with the Hadamard matrix and multiply + # Cost: 8S 7M + p1, p2, p3, p4 = P.squared_theta() + q1, q2, q3, q4 = Q.squared_theta() + + xp = p1 * q1 + yp = Y0 * p2 * q2 + zp = Z0 * p3 * q3 + tp = T0 * p4 * q4 + + # Final coordinates + PQx, PQy, PQz, PQt = PQ.coords() + + # Note: + # We replace the four divisions by + # PQx, PQy, PQz, PQt by 10 multiplications + # Cost: 10M + PQxy = PQx * PQy + PQzt = PQz * PQt + + X, Y, Z, T = P.to_hadamard(xp, yp, zp, tp) + X = X * PQzt * PQy + Y = Y * PQzt * PQx + Z = Z * PQxy * PQt + T = T * PQxy * PQz + + coords = (X, Y, Z, T) + return P.parent()(coords) + + def scale(self, n): + """ + Scale all coordinates of the ThetaPointDim2 by `n` + """ + x, y, z, t = self.coords() + if not isinstance(n, RingElement): + raise ValueError(f"Cannot scale by element {n} of type {type(n)}") + scaled_coords = (n * x, n * y, n * z, n * t) + return self._parent(scaled_coords) + + def double_iter(self, m): + """ + Compute [2^n] Self + + NOTE: Assumes that no coordinate is zero at any point during the doubling + """ + if not isinstance(m, Integer): + try: + m = Integer(m) + except: + raise TypeError(f"Cannot coerce input scalar {m = } to an integer") + + if m.is_zero(): + return self.parent().zero() + + P1 = self + for _ in range(m): + P1 = P1.double() + return P1 + + def __mul__(self, m): + """ + Uses Montgomery ladder to compute [m] Self + + NOTE: Assumes that no coordinate is zero at any point during the doubling + """ + # Make sure we're multiplying by something value + if not isinstance(m, (int, Integer)): + try: + m = Integer(m) + except: + raise TypeError(f"Cannot coerce input scalar {m = } to an integer") + + # If m is zero, return the null point + if not m: + return self.parent().zero() + + # We are with ±1 identified, so we take the absolute value of m + m = abs(m) + + P0, P1 = self, self + P2 = P1.double() + # If we are multiplying by two, the chain stops here + if m == 2: + return P2 + + # Montgomery double and add. + for bit in bin(m)[3:]: + Q = P2.diff_addition(P1, P0) + if bit == "1": + P2 = P2.double() + P1 = Q + else: + P1 = P1.double() + P2 = Q + + return P1 + + def __rmul__(self, m): + return self * m + + def __imul__(self, m): + self = self * m + return self + + def __eq__(self, other): + """ + Check the quality of two ThetaPoints. Note that as this is a + projective equality, we must be careful for when certain coefficients may + be zero. + """ + if not isinstance(other, ThetaPointDim2): + return False + + a1, b1, c1, d1 = self.coords() + a2, b2, c2, d2 = other.coords() + + if d1 != 0 or d2 != 0: + return all([a1 * d2 == a2 * d1, b1 * d2 == b2 * d1, c1 * d2 == c2 * d1]) + elif c1 != 0 or c2 != 0: + return all([a1 * c2 == a2 * c1, b1 * c2 == b2 * c1]) + elif b1 != 0 or b2 != 0: + return a1 * b2 == a2 * b1 + else: + return True + + def __repr__(self): + return f"Theta point with coordinates: {self.coords()}" diff --git a/theta_lib/theta_structures/Theta_dim4.py b/theta_lib/theta_structures/Theta_dim4.py new file mode 100644 index 0000000..7851c56 --- /dev/null +++ b/theta_lib/theta_structures/Theta_dim4.py @@ -0,0 +1,351 @@ +from sage.all import * +from sage.structure.element import get_coercion_model + +from .theta_helpers_dim4 import ( + hadamard, + batch_inversion, + product_theta_point_dim4, + product_to_theta_points_dim4, + product_theta_point_dim2_dim4, + product_to_theta_points_dim4_dim2, + act_point, + squared, +) +from .Theta_dim1 import ThetaStructureDim1 +from .Tuple_point import TuplePoint +from ..basis_change.base_change_dim4 import ( + apply_base_change_theta_dim4, + random_symplectic_matrix, + base_change_theta_dim4, +) + +cm = get_coercion_model() + + +class ThetaStructureDim4: + def __init__(self,null_point,null_point_dual=None,inv_null_point_dual_sq=None): + r""" + INPUT: + - null_point: theta-constants. + - inv_null_point_dual_sq: inverse of the squares of dual theta-constants, if provided + (meant to prevent duplicate computation, since this data is already computed when the + codomain of an isogeny is computed). + """ + if not len(null_point) == 16: + raise ValueError("Entry null_point should have 16 coordinates.") + + self._base_ring = cm.common_parent(*(c.parent() for c in null_point)) + self._point = ThetaPointDim4 + self._null_point = self._point(self, null_point) + self._null_point_dual=null_point_dual + self._inv_null_point=None + self._inv_null_point_dual_sq=inv_null_point_dual_sq + + def null_point(self): + """ + """ + return self._null_point + + def null_point_dual(self): + if self._null_point_dual==None: + self._null_point_dual=hadamard(self._null_point.coords()) + return self._null_point_dual + + def base_ring(self): + """ + """ + return self._base_ring + + def zero(self): + """ + """ + return self.null_point() + + def zero_dual(self): + return self.null_point_dual() + + def __repr__(self): + return f"Theta structure over {self.base_ring()} with null point: {self.null_point()}" + + def __call__(self,coords): + return self._point(self,coords) + + def act_null(self,I,J): + r""" + Point of 2-torsion. + + INPUT: + - I, J: two 4-tuples of indices in {0,1}. + + OUTPUT: the action of (I,\chi_J) on the theta null point given by: + (I,\chi_J)*P=(\chi_J(I+K)^{-1}P[I+K])_K + """ + return self.null_point().act_point(I,J) + + def is_K2(self,B): + r""" + Given a symplectic decomposition A[2]=K_1\oplus K_2 canonically + induced by the theta-null point, determines if B is the canonical + basis of K_2 given by act_nul(0,\delta_i)_{1\leq i\leq 4}. + + INPUT: + - B: Basis of 4 points of 2-torsion. + + OUTPUT: Boolean True if and only if B is the canonical basis of K_2. + """ + I0=(0,0,0,0) + if B[0]!=self.act_null(I0,(1,0,0,0)): + return False + if B[1]!=self.act_null(I0,(0,1,0,0)): + return False + if B[2]!=self.act_null(I0,(0,0,1,0)): + return False + if B[3]!=self.act_null(I0,(0,0,0,1)): + return False + return True + + def base_change_struct(self,N): + null_coords=self.null_point().coords() + new_null_coords=apply_base_change_theta_dim4(N,null_coords) + return ThetaStructureDim4(new_null_coords) + + def base_change_coords(self,N,P): + coords=P.coords() + new_coords=apply_base_change_theta_dim4(N,coords) + return self.__call__(new_coords) + + #@cached_method + def _arithmetic_precomputation(self): + r""" + Initializes the precomputation containing the inverse of the theta-constants in standard + and dual (Hadamard transformed) theta-coordinates. Assumes no theta-constant is zero. + """ + O=self.null_point() + if all([O[k]!=0 for k in range(16)]): + self._inv_null_point=batch_inversion(O.coords()) + if self._inv_null_point_dual_sq==None: + U_chi_0_sq=hadamard(squared(O.coords())) + if all([U_chi_0_sq[k]!=0 for k in range(16)]): + self._inv_null_point_dual_sq=batch_inversion(U_chi_0_sq) + self._arith_base_change=False + else: + self._arith_base_change=True + self._arithmetic_base_change() + #print("Zero dual theta constants.\nRandom symplectic base change is being used for duplication.\nDoublings are more costly than expected.") + else: + self._arith_base_change=True + self._arithmetic_base_change() + #print("Zero theta constants.\nRandom symplectic base change is being used for duplication.\nDoublings are more costly than expected.") + + def _arithmetic_base_change(self,max_iter=50): + F=self._base_ring + if F.degree() == 2: + i=self._base_ring.gen() + else: + assert(F.degree() == 1) + Fp2 = GF((F.characteristic(), 2), name='i', modulus=var('x')**2 + 1) + i=Fp2.gen() + + count=0 + O=self.null_point() + while count<max_iter: + count+=1 + M=random_symplectic_matrix(4) + N=base_change_theta_dim4(M,i) + + NO=apply_base_change_theta_dim4(N,O) + NU_chi_0_sq=hadamard(squared(NO)) + + if all([NU_chi_0_sq[k]!=0 for k in range(16)]) and all([NO[k]!=0 for k in range(16)]): + self._arith_base_change_matrix=N + self._arith_base_change_matrix_inv=N.inverse() + self._inv_null_point=batch_inversion(NO) + self._inv_null_point_dual_sq=batch_inversion(NU_chi_0_sq) + break + + def has_suitable_doubling(self): + O=self.null_point() + UO=hadamard(O.coords()) + if all([O[k]!=0 for k in range(16)]) and all([UO[k]!=0 for k in range(16)]): + return True + else: + return False + + def hadamard(self): + return ThetaStructureDim4(self.null_point_dual()) + + def hadamard_change_coords(self,P): + new_coords=hadamard(P) + return self.__call__(new_coords) + + +class ProductThetaStructureDim1To4(ThetaStructureDim4): + def __init__(self,*args): + r"""Defines the product theta structure at level 2 of 4 elliptic curves. + + Input: Either + - 4 theta structures of dimension 1: T0, T1, T2, T3; + - 4 elliptic curves: E0, E1, E2, E3. + - 4 elliptic curves E0, E1, E2, E3 and their respective canonical 4-torsion basis B0, B1, B2, B3. + """ + if len(args)==4: + theta_structures=list(args) + for k in range(4): + if not isinstance(theta_structures[k],ThetaStructureDim1): + try: + theta_structures[k]=ThetaStructureDim1(theta_structures[k]) + except: + pass + elif len(args)==8: + theta_structures=[ThetaStructureDim1(args[k],args[4+k][0],args[4+k][1]) for k in range(4)] + else: + raise ValueError("4 or 8 arguments expected but {} were given.\nYou should enter a list of 4 elliptic curves or ThetaStructureDim1\nor a list of 4 elliptic curves with a 4-torsion basis for each of them.".format(len(args))) + + self._theta_structures=theta_structures + + null_point=product_theta_point_dim4(theta_structures[0].zero().coords(),theta_structures[1].zero().coords(), + theta_structures[2].zero().coords(),theta_structures[3].zero().coords()) + + ThetaStructureDim4.__init__(self,null_point) + + def product_theta_point(self,theta_points): + return self._point(self,product_theta_point_dim4(theta_points[0].coords(),theta_points[1].coords(), + theta_points[2].coords(),theta_points[3].coords())) + + def __call__(self,point): + if isinstance(point,TuplePoint): + theta_points=[] + theta_structures=self._theta_structures + for i in range(4): + theta_points.append(theta_structures[i](point[i])) + return self.product_theta_point(theta_points) + else: + return self._point(self,point) + + def to_theta_points(self,P): + theta_coords=product_to_theta_points_dim4(P) + theta_points=[self._theta_structures[i](theta_coords[i]) for i in range(4)] + return theta_points + + def to_tuple_point(self,P): + theta_points=self.to_theta_points(P) + montgomery_points=[self._theta_structures[i].to_montgomery_point(theta_points[i]) for i in range(4)] + return TuplePoint(montgomery_points) + +class ProductThetaStructureDim2To4(ThetaStructureDim4): + def __init__(self,theta1,theta2): + self._theta_structures=(theta1,theta2) + + null_point=product_theta_point_dim2_dim4(theta1.zero().coords(),theta2.zero().coords()) + + ThetaStructureDim4.__init__(self,null_point) + + def product_theta_point(self,P1,P2): + return self._point(self,product_theta_point_dim2_dim4(P1.coords(),P2.coords())) + + def __call__(self,point): + return self._point(self,point) + + def to_theta_points(self,P): + theta_coords=product_to_theta_points_dim4_dim2(P) + theta_points=[self._theta_structures[i](theta_coords[i]) for i in range(2)] + return theta_points + +class ThetaPointDim4: + def __init__(self, parent, coords): + """ + """ + if not isinstance(parent, ThetaStructureDim4): + raise ValueError("Entry parent should be a ThetaStructureDim4 object.") + + self._parent = parent + self._coords = tuple(coords) + + def parent(self): + """ + """ + return self._parent + + def theta(self): + """ + """ + return self.parent() + + def coords(self): + """ + """ + return self._coords + + def is_zero(self): + """ + """ + return self == self.parent().zero() + + def __eq__(self, other): + P=self.coords() + Q=other.coords() + + k0=0 + while k0<15 and P[k0]==0: + k0+=1 + + for l in range(16): + if P[l]*Q[k0]!=Q[l]*P[k0]: + return False + return True + + def __repr__(self): + return f"Theta point with coordinates: {self.coords()}" + + def __getitem__(self,i): + return self._coords[i] + + def scale(self,lamb): + if lamb==0: + raise ValueError("Entry lamb should be non-zero.") + + P=self.coords() + return self._parent([lamb*x for x in P]) + + def act_point(self,I,J): + r""" + Translation by a point of 2-torsion. + + Input: + - I, J: two 4-tuples of indices in {0,1}. + + Output: the action of (I,\chi_J) on P given by: + (I,\chi_J)*P=(\chi_J(I+K)^{-1}P[I+K])_K + """ + return self._parent(act_point(self._coords,I,J)) + + def double(self): + ## This formula is projective. + ## Works only when theta constants are non-zero. + P=self.coords() + if self.parent()._inv_null_point==None or self.parent()._inv_null_point_dual_sq==None: + self.parent()._arithmetic_precomputation() + inv_O,inv_U_chi_0_sq=self.parent()._inv_null_point,self.parent()._inv_null_point_dual_sq + + if self.parent()._arith_base_change: + P=apply_base_change_theta_dim4(self.parent()._arith_base_change_matrix,P) + + U_chi_P=squared(hadamard(squared(P))) + for chi in range(16): + U_chi_P[chi]*=inv_U_chi_0_sq[chi] + + theta_2P = list(hadamard(U_chi_P)) + for i in range(16): + theta_2P[i] *= inv_O[i] + + if self.parent()._arith_base_change: + theta_2P=apply_base_change_theta_dim4(self.parent()._arith_base_change_matrix_inv,theta_2P) + + return self._parent(theta_2P) + + def double_iter(self,n): + ## Computes 2**n*self + Q=self + for i in range(n): + Q=Q.double() + return Q diff --git a/theta_lib/theta_structures/Tuple_point.py b/theta_lib/theta_structures/Tuple_point.py new file mode 100644 index 0000000..aac248b --- /dev/null +++ b/theta_lib/theta_structures/Tuple_point.py @@ -0,0 +1,106 @@ +from sage.all import * +from ..utilities.discrete_log import weil_pairing_pari + +class TuplePoint: + def __init__(self,*args): + if len(args)==1: + self._points=list(args[0]) + else: + self._points=list(args) + + def points(self): + return self._points + + def parent_curves(self): + return [x.curve() for x in self._points] + + def parent_curve(self,i): + return self._points[i].curve() + + def n_points(self): + return len(self._points) + + def is_zero(self): + return all([self._points[i]==0 for i in range(self.n_points())]) + + def __repr__(self): + return str(self._points) + + def __getitem__(self,i): + return self._points[i] + + def __setitem__(self,i,P): + self._points[i]=P + + def __eq__(self,other): + n_self=self.n_points() + n_other=self.n_points() + return n_self==n_other and all([self._points[i]==other._points[i] for i in range(n_self)]) + + def __add__(self,other): + n_self=self.n_points() + n_other=self.n_points() + + if n_self!=n_other: + raise ValueError("Cannot add TuplePoint of distinct lengths {} and {}.".format(n_self,n_other)) + + points=[] + for i in range(n_self): + points.append(self._points[i]+other._points[i]) + return self.__class__(points) + + def __sub__(self,other): + n_self=self.n_points() + n_other=self.n_points() + + if n_self!=n_other: + raise ValueError("Cannot substract TuplePoint of distinct lengths {} and {}.".format(n_self,n_other)) + + points=[] + for i in range(n_self): + points.append(self._points[i]-other._points[i]) + return self.__class__(points) + + def __neg__(self): + n_self=self.n_points() + points=[] + for i in range(n_self): + points.append(-self._points[i]) + return self.__class__(points) + + def __mul__(self,m): + n_self=self.n_points() + points=[] + for i in range(n_self): + points.append(m*self._points[i]) + return self.__class__(points) + + def __rmul__(self,m): + return self*m + + def double_iter(self,n): + result=self + for i in range(n): + result=2*result + return result + + def weil_pairing(self,other,n): + n_self=self.n_points() + n_other=self.n_points() + + if n_self!=n_other: + raise ValueError("Cannot compute the Weil pairing of TuplePoint of distinct lengths {} and {}.".format(n_self,n_other)) + + zeta=1 + for i in range(n_self): + zeta*=weil_pairing_pari(self._points[i],other._points[i],n) + + return zeta + + + + + + + + diff --git a/theta_lib/theta_structures/montgomery_theta.py b/theta_lib/theta_structures/montgomery_theta.py new file mode 100644 index 0000000..6dddb2b --- /dev/null +++ b/theta_lib/theta_structures/montgomery_theta.py @@ -0,0 +1,68 @@ +from sage.all import * + +def torsion_to_theta_null_point(P): + r=P[0] + s=P[2] + return (r+s,r-s) + +def montgomery_point_to_theta_point(O,P): + if P[0]==P[2]==0: + return O + else: + a,b=O + return (a*(P[0]-P[2]),b*(P[0]+P[2])) + +def theta_point_to_montgomery_point(E,O,P,twist=False): + a,b=O + + x=a*P[1]+b*P[0] + z=a*P[1]-b*P[0] + + if twist: + x=-x + + if z==0: + return E(0) + else: + x=x/z + + a_inv=E.a_invariants() + + A =a_inv[1] + if a_inv != (0,A,0,1,0): + raise ValueError("The elliptic curve E is not in the Montgomery model.") + + y2=x*(x**2+A*x+1) + if not is_square(y2): + raise ValueError("The Montgomery point is not on the base field.") + else: + y=y2.sqrt() + return E([x,y,1]) + +def lift_kummer_montgomery_point(E,x,z=1): + if z==0: + return E(0) + elif z!=1: + x=x/z + + a_inv=E.a_invariants() + + A =a_inv[1] + if a_inv != (0,A,0,1,0): + raise ValueError("The elliptic curve E is not in the Montgomery model.") + + y2=x*(x**2+A*x+1) + if not is_square(y2): + raise ValueError("The Montgomery point is not on the base field.") + else: + y=y2.sqrt() + return E([x,y,1]) + +def null_point_to_montgomery_coeff(a,b): + return -2*(a**4+b**4)/(a**4-b**4) + + + + + + diff --git a/theta_lib/theta_structures/theta_helpers_dim2.py b/theta_lib/theta_structures/theta_helpers_dim2.py new file mode 100644 index 0000000..a57789d --- /dev/null +++ b/theta_lib/theta_structures/theta_helpers_dim2.py @@ -0,0 +1,32 @@ +from sage.all import * + +def batch_inversion(L): + r"""Does n inversions in 3(n-1)M+1I. + + Input: + - L: list of elements to invert. + + Output: + - [1/x for x in L] + """ + # Given L=[a0,...,an] + # Computes multiples=[a0, a0.a1, ..., a0...an] + multiples=[L[0]] + for ai in L[1:]: + multiples.append(multiples[-1]*ai) + + # Computes inverses=[1/(a0...an),...,1/a0] + inverses=[1/multiples[-1]] + for i in range(1,len(L)): + inverses.append(inverses[-1]*L[-i]) + + # Finally computes [1/a0,...,1/an] + result=[inverses[-1]] + for i in range(2,len(L)+1): + result.append(inverses[-i]*multiples[i-2]) + return result + +def product_theta_point(theta1,theta2): + a,b=theta1 + c,d=theta2 + return [a*c,b*c,a*d,b*d] diff --git a/theta_lib/theta_structures/theta_helpers_dim4.py b/theta_lib/theta_structures/theta_helpers_dim4.py new file mode 100644 index 0000000..29599b0 --- /dev/null +++ b/theta_lib/theta_structures/theta_helpers_dim4.py @@ -0,0 +1,259 @@ +from sage.all import * +import itertools + +## Index management +@cached_function +def index_to_multindex(k): + r""" + Input: + - k: integer between 0 and 15. + + Output: binary decomposition of k. + """ + L_ind=[] + l=k + for i in range(4): + L_ind.append(l%2) + l=l//2 + return tuple(L_ind) + +@cached_function +def multindex_to_index(*args): + r""" + Input: 4 elements i0,i1,i2,i3 in {0,1}. + + Output: k=i0+2*i1+4*i2+8*i3. + """ + if len(args)==4: + i0,i1,i2,i3=args + else: + i0,i1,i2,i3=args[0] + return i0+2*i1+4*i2+8*i3 + +@cached_function +def scal_prod(i,j): + r""" + Input: Two integers i and j in {0,...,15}. + + Output: Scalar product of the bits of i and j mod 2. + """ + return (int(i)&int(j)).bit_count()%2 + +def act_point(P,I,J): + r""" + Input: + - P: a point with 16 coordinates. + - I, J: two 4-tuples of indices in {0,1}. + + Output: the action of (I,\chi_J) on P given by: + (I,\chi_J)*P=(\chi_J(I+K)^{-1}P[I+K])_K + """ + Q=[] + i=multindex_to_index(I) + j=multindex_to_index(J) + for k in range(16): + ipk=i^k + Q.append((-1)**scal_prod(ipk,j)*P[ipk]) + return Q + +## Product of theta points +def product_theta_point_dim4(P0,P1,P2,P3): + # Computes the product theta coordinates of a product of 4 elliptic curves. + P=[0 for k in range(16)] + for i0, i1, i2, i3 in itertools.product([0,1],repeat=4): + P[multindex_to_index(i0,i1,i2,i3)]=P0[i0]*P1[i1]*P2[i2]*P3[i3] + return P + +def product_theta_point_dim2_dim4(P0,P1): + # Computes the product theta coordinates of a product of 2 abelian surfaces. + P=[0 for k in range(16)] + for i0, i1, i2, i3 in itertools.product([0,1],repeat=4): + P[multindex_to_index(i0,i1,i2,i3)]=P0[i0+2*i1]*P1[i2+2*i3] + return P + +## 4-dimensional product Theta point to 1-dimensional Theta points +def product_to_theta_points_dim4(P): + Fp2=P[0].parent() + d_Pi={0:None,1:None,2:None,3:None} + d_index_ratios={0:None,1:None,2:None,3:None}# Index of numertors/denominators to compute the theta points. + for k in range(4): + is_zero=True# theta_1(Pk)=0 + for I in itertools.product([0,1],repeat=3): + J=list(I) + J.insert(k,1) + j=multindex_to_index(*J) + if P[j]!=0: + is_zero=False + d_index_ratios[k]=[j^(2**k),j] + break + if is_zero: + d_Pi[k]=(Fp2(1),Fp2(0)) + L_num=[] + L_denom=[] + d_index_num_denom={} + for k in range(4): + if d_Pi[k]==None:# Point has non-zero coordinate theta_1 + d_index_num_denom[k]=len(L_num) + L_num.append(P[d_index_ratios[k][0]]) + L_denom.append(P[d_index_ratios[k][1]]) + L_denom=batch_inversion(L_denom) + for k in range(4): + if d_Pi[k]==None: + d_Pi[k]=(L_num[d_index_num_denom[k]]*L_denom[d_index_num_denom[k]],Fp2(1)) + return d_Pi[0],d_Pi[1],d_Pi[2],d_Pi[3] + +## 4-dimensional product Theta point to 2-dimensional Theta points +def product_to_theta_points_dim4_dim2(P): + Fp2=P[0].parent() + k0=0# Index of the denominator k0=multindex_to_index(I0,J0) and + # we divide by theta_{I0,J0}=theta_{I0}*theta_{J0}!=0 + while k0<=15 and P[k0]==0: + k0+=1 + i0, j0 = k0%4, k0//4 + inv_theta_k0=1/P[k0] + theta_P1=[] + theta_P2=[] + for i in range(4): + if i==i0: + theta_P1.append(1) + else: + theta_P1.append(P[i+4*j0]*inv_theta_k0) + for j in range(4): + if j==j0: + theta_P2.append(1) + else: + theta_P2.append(P[i0+4*j]*inv_theta_k0) + return theta_P1, theta_P2 + + + + + +## Usual theta transformations +@cached_function +def inv_16(F): + return 1/F(16) + +def hadamard2(x,y): + return (x+y, x-y) + +def hadamard4(x,y,z,t): + x,y=hadamard2(x,y) + z,t=hadamard2(z,t) + return (x+z, y+t, x-z, y-t) + +def hadamard8(a,b,c,d,e,f,g,h): + a,b,c,d=hadamard4(a,b,c,d) + e,f,g,h=hadamard4(e,f,g,h) + return (a+e, b+f, c+g, d+h, a-e, b-f, c-g, d-h) + +def hadamard16(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p): + a,b,c,d,e,f,g,h=hadamard8(a,b,c,d,e,f,g,h) + i,j,k,l,m,n,o,p=hadamard8(i,j,k,l,m,n,o,p) + return (a+i, b+j, c+k, d+l, e+m, f+n, g+o, h+p, a-i, b-j, c-k, d-l, e-m, f-n, g-o, h-p) + +def hadamard_inline(a,b,c,d,e,f,g,h,i,j,k,l,m,n,o,p): + a,b=a+b,a-b + c,d=c+d,c-d + e,f=e+f,e-f + g,h=g+h,g-h + i,j=i+j,i-j + k,l=k+l,k-l + m,n=m+n,m-n + o,p=o+p,o-p + a,b,c,d=a+c,b+d,a-c,b-d + e,f,g,h=e+g,f+h,e-g,f-h + i,j,k,l=i+k,j+l,i-k,j-l + m,n,o,p=m+o,n+p,m-o,n-p + a,b,c,d,e,f,g,h=a+e, b+f, c+g, d+h, a-e, b-f, c-g, d-h + i,j,k,l,m,n,o,p=i+m,j+n,k+o,l+p,i-m,j-n,k-o,l-p + return (a+i, b+j, c+k, d+l, e+m, f+n, g+o, h+p, a-i, b-j, c-k, d-l, e-m, f-n, g-o, h-p) + +def hadamard(P): + return hadamard16(*P) + #return hadamard_inline(*P) + +def hadamard_inverse(P): + H_inv_P=[] + C=inv_16(P[0].parent()) + for j in range(16): + HP.append(0) + for k in range(16): + if scal_prod(k,j)==0: + H_inv_P[j]+=P[k] + else: + H_inv_P[j]-=P[k] + H_inv_P[j]=H_inv_P[j]*C + return H_inv_P + +def squared(P): + return [x**2 for x in P] + +def batch_inversion(L): + r"""Does n inversions in 3(n-1)M+1I. + + Input: + - L: list of elements to invert. + + Output: + - [1/x for x in L] + """ + # Given L=[a0,...,an] + # Computes multiples=[a0, a0.a1, ..., a0...an] + multiples=[L[0]] + for ai in L[1:]: + multiples.append(multiples[-1]*ai) + + # Computes inverses=[1/(a0...an),...,1/a0] + inverses=[1/multiples[-1]] + for i in range(1,len(L)): + inverses.append(inverses[-1]*L[-i]) + + # Finally computes [1/a0,...,1/an] + result=[inverses[-1]] + for i in range(2,len(L)+1): + result.append(inverses[-i]*multiples[i-2]) + return result + + +## Functions to handle zero theta coordinates +def find_zeros(P): + L_ind_zeros=[] + for i in range(16): + if P[i]==0: + L_ind_zeros.append(i) + return L_ind_zeros + +def find_translates(L_ind_zeros): + L_ind_non_zero=[] + L_ind_origin=L_ind_zeros.copy() + + for i in range(16): + if i not in L_ind_zeros: + L_ind_non_zero.append(i) + + L_trans=[] + while L_ind_origin!=[]: + n_target_max=0 + ind_trans_max=0 + for k in range(16): + trans=[x^k for x in L_ind_origin] + n_target=0 + for y in trans: + if y in L_ind_non_zero: + n_target+=1 + if n_target>n_target_max: + n_target_max=n_target + ind_trans_max=k + L_trans.append(ind_trans_max) + L_ind_remove=[] + for x in L_ind_origin: + if x^ind_trans_max in L_ind_non_zero: + L_ind_remove.append(x) + for x in L_ind_remove: + L_ind_origin.remove(x) + return L_trans + + + + diff --git a/theta_lib/utilities/discrete_log.py b/theta_lib/utilities/discrete_log.py new file mode 100644 index 0000000..dff04dc --- /dev/null +++ b/theta_lib/utilities/discrete_log.py @@ -0,0 +1,273 @@ +# ===================================== # +# 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 new file mode 100644 index 0000000..0609fe5 --- /dev/null +++ b/theta_lib/utilities/fast_sqrt.py @@ -0,0 +1,55 @@ + +# ============================================ # +# 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 new file mode 100644 index 0000000..223f568 --- /dev/null +++ b/theta_lib/utilities/order.py @@ -0,0 +1,117 @@ +# ================================================== # +# 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 new file mode 100644 index 0000000..550ef09 --- /dev/null +++ b/theta_lib/utilities/strategy.py @@ -0,0 +1,229 @@ +# ============================================================================ # +# 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 new file mode 100644 index 0000000..e760c1b --- /dev/null +++ b/theta_lib/utilities/supersingular.py @@ -0,0 +1,413 @@ +""" +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 |