From d40de259097c5e8d8fd35539560ca7c3d47523e7 Mon Sep 17 00:00:00 2001 From: Ryan Rueger Date: Sat, 1 Mar 2025 20:25:41 +0100 Subject: Initial Commit MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-Authored-By: Damien Robert Co-Authored-By: Frederik Vercauteren Co-Authored-By: Jonathan Komada Eriksen Co-Authored-By: Pierrick Dartois Co-Authored-By: Riccardo Invernizzi Co-Authored-By: Ryan Rueger [0.01s] Co-Authored-By: Benjamin Wesolowski Co-Authored-By: Arthur Herlédan Le Merdy Co-Authored-By: Boris Fouotsa --- elkies.py | 228 ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 228 insertions(+) create mode 100644 elkies.py (limited to 'elkies.py') diff --git a/elkies.py b/elkies.py new file mode 100644 index 0000000..aa6e8d3 --- /dev/null +++ b/elkies.py @@ -0,0 +1,228 @@ +from sage.all import * +import random +from xonly import xPoint + +def Elkies(E, ell, Fp, lam = None, prev_j = None): + """ + Elkies algorithm for computing isogenies as detailed in [BSS - Elliptic + Curves in Cryptography, Chapter VII]. + Input: + - E: starting curve + - ell: a split prime + - Fp: base field + - lam: if looking for a specific Elkies isogeny, optionally provide the + corresponding eigenvalue + - prev_j: if the curve is the codomain of another ell-isogeny, + optionally provide the j invariant of the previous curve to avoid + backtracking + Output: + - the kernel polynomial hc + """ + #TODO maybe we want montomery curves later + p = Fp.characteristic() + assert ell > 2 + d = (ell-1)/2 + j = Fp(E.j_invariant()) + _, _, _, a, b = [Fp(c) for c in E.a_invariants()] + + R = Fp['X,Y'] + X, Y = R._first_ngens(2) + + E4 = -48*a + E6 = 864*b + + jp = -E6*j/E4 + + phi_ell = R(classical_modular_polynomial(ell)) + X, Y = phi_ell.variables() + phi_ellx = derivative(phi_ell, X) + phi_ellxx = derivative(phi_ellx, X) + phi_ellxy = derivative(phi_ellx, Y) + + phi_elly = derivative(phi_ell, Y) + phi_ellyy = derivative(phi_elly, Y) + + phi_ell_eval = phi_ell(X, j).univariate_polynomial() + if prev_j: + x = phi_ell_eval.variables()[0] + phi_ell_eval = phi_ell_eval // (x - prev_j) + # print(f"Finding j: {time()-tstart}") + if prev_j: + j_ells = [phi_ell_eval.any_root()] #There is only one root in this case! + #assert len(j_ells) == 1, f"j_ells was {j_ells}, and prev_j was {prev_j}" + else: + j_ells = [r for r,e in phi_ell_eval.roots()] + assert len(j_ells) == 2, f"j_ells was {j_ells}" + + #print(j_ells) + jt = j_ells[0] + + assert jt != 0 and jt != 1728, "encountered E0" + + hc = _derive_hc(Fp, a, b, ell, d, j, jt, jp, E4, E6, phi_ellx, phi_ellxx, + phi_ellxy, phi_elly, phi_ellyy) + + if lam and not prev_j: + if not CheckElkies(E, ell, hc, lam): + #Do the other one + assert j_ells[1] != 0 and j_ells[1] != 1728, "encountered E0" + hc = _derive_hc(Fp, a, b, ell, d, j, j_ells[1], jp, E4, E6, phi_ellx, phi_ellxx, phi_ellxy, phi_elly, phi_ellyy) + return hc + +def WeierstrassP(a,b,d): + # Weierstrass p up to degree d in w=z^2 + coefs = [-a/5, -b/7] + for k in range(2, d): + ck = 0 + for j in range (k-1): + ck = ck + coefs[j]*coefs[k-2-j] + ck = ck*3/((k-1)*(2*(k+1)+3)) + coefs.append(ck) + return coefs + + +def _derive_hc(Fp, a, b, ell, d, j, jt, jp, E4, E6, phi_ellx, phi_ellxx, phi_ellxy, phi_elly, phi_ellyy): + jtp = -jp*phi_ellx(j, jt)/(ell*phi_elly(j, jt)) + at = -jtp**2/(48*jt*(jt-1728)) + bt = -jtp**3/(864*jt**2*(jt-1728)) + + E4t = -48*at + E6t = 864*bt + + fracjs = -(jp**2*phi_ellxx(j, jt) + 2*ell*jp*jtp*phi_ellxy(j, jt) + (ell*jtp)**2*phi_ellyy(j, jt))/ (jp*phi_ellx(j, jt)) + p1 = ell*fracjs/2 + (Fp(ell)/4)*(E4**2/E6 - ell*E4t**2/E6t) + (Fp(ell)/3)*(E6/E4 - ell*E6t/E4t) + + c = WeierstrassP(a,b,d) + ckst = WeierstrassP(ell**4*at, ell**6*bt, d) + ct = [ckst[i] - ell*c[i] for i in range(len(c))] # difference in formula VII.23 + + # Computing the coefficients of VII.23 and store as A[i] + Fpw, w = PowerSeriesRing(Fp, 'w', default_prec=d+1).objgen() + evp = -(p1/2)*w - sum(w**(i+1)*ct[i-1] / ((2*i+1)*(2*i+2)) for i in range(1,d)) + exp_evp = evp.exp(d+1) + A = exp_evp.coefficients() + C = sum(c[i-1]*w**i for i in range(1, d+1)) + + # Computing all powers of C starting with zeroth power + Cpow = [Fpw(1), C]; + for i in range(2, d+1): + Cpow.append(C*Cpow[-1]) + + # Now doing recurrence relation VII.24 + hc = [1, -p1/2] + for i in range(2, d+1): + newcoeff = A[i] + for k in range(1, i+1): + insum = 0 + for j in range(k+1): + insum += Fp(binomial(d-i+k, k-j))*Cpow[k-j][j] + newcoeff -= insum*hc[i-k] + hc.append(newcoeff) + + Rx, x = PolynomialRing(Fp, 'x').objgen() + hc = Rx(hc[::-1]) + return hc + +def CheckElkies(E, ell, h, lam): + p = E.base_field().characteristic() + if kronecker_symbol(-p, ell) != 1: + assert False, "not Elkies" + + _, _, _, a, b = E.a_invariants() + f = h.parent()([b, a, 0, 1]) + + B = pow(f, (p-1)/2, h) + check_wrong = False + if lam > ell//2: + lam = ell-lam + check_wrong = True + + if lam == 1: + if check_wrong: + return B != 1 + else: + return B == 1 + + # Stupid way for now, no sage function for directly computing mod h + RR = h.parent().quotient_ring(h) + y_coord = E.multiplication_by_m(lam)[1] + y_coord_num = y_coord.numerator() + x, y = y_coord_num.variables() + y_coord_den = y_coord.denominator() + + if check_wrong: + return RR(B)*RR(y_coord_den) != RR(y_coord_num(x, 1).univariate_polynomial()) + else: + return RR(B)*RR(y_coord_den) == RR(y_coord_num(x, 1).univariate_polynomial()) + + +def NonElkiesIsogeny(E, ell, Fp2): + # For u, v, when not elkies prime. + # Make that either ell divides p+1, the other way works... + + E = E.change_ring(Fp2) + #print(factor(E.division_polynomial(ell))) + # Make sure these are the only cases... + h_ell = factor(E.division_polynomial(ell))[0][0] + + m = (ell-1)//(2*h_ell.degree()) + if m == 1: + return h_ell + + #Deuring for the people <3 + from sage.schemes.elliptic_curves.isogeny_small_degree import _least_semi_primitive + a = _least_semi_primitive(ell) + fs = [h_ell] + Fbig = Fp2.extension(h_ell) + x = Fbig.gens()[0] + xi = xPoint(x, E.change_ring(Fbig)) + for _ in range(1, m): + xi = xi.mul(a) + fs.append(xi.X.minpoly()) + + h_ell = prod(fs) + assert h_ell.degree() == (ell-1)/2, f"Degree : {h_ell.degree()}, shouldve been {(ell-1)/2}" + return h_ell + + + +### This is way slower for now :(( +def ElkiesDirect(E, ell, lam = None): + psi_x = E.division_polynomial(ell) + Fp = E.base_field() + p = Fp.characteristic() + if not lam: + F_ell = GF(ell) + lam = F_ell(-p).sqrt() + + R = psi_x.parent() + x = psi_x.variables()[0] + RR = R.quotient_ring(psi_x) + + X = RR._first_ngens(1)[0] + + check_wrong = False + if lam > ell//2: + lam = ell-lam + check_wrong = True + + if lam == 1: + h_sq_bar = X**p - X + else: + mult_by_lam = E.multiplication_by_m(lam)[0] + mult_by_lam_num = mult_by_lam.numerator().univariate_polynomial() + mult_by_lam_denom = mult_by_lam.denominator().univariate_polynomial() + h_sq_bar = X**p*mult_by_lam_denom(X) - mult_by_lam_num(X) + + h_sq = 0 + x_pow = x**0 + for ai in list(h_sq_bar): + h_sq += ai*x_pow + x_pow *= x + h_sq = gcd(psi_x, h_sq) + #print(h_sq) + #print(factor(h_sq)) + + #print(Elkies(E, ell, Fp)) + + -- cgit v1.2.3-70-g09d2