""" NOTE: Follow SIKE's spec in the generation of E[2^a] torsion basis. It will allow you to compute the isogeny faster. """
sqrtof2 = Fp2(2).sqrt() f = x**3 + A0 * x**2 + x
Pa = E0(0) Qa = E0(0) Pa_done = False Qa_done = False d = 0 for c inrange(0, p): Rx = ii + c Ry_square = f(ii + c) ifnot Ry_square.is_square(): continue Ra = E0.lift_x(Rx) Pa = 3**b * Ra
Ta = 2 ** (a - 1) * Pa if Ta.is_zero(): continue Tax_plus_3 = Ta.xy()[0] + 3 if Tax_plus_3 == 2 * sqrtof2 or Tax_plus_3 == -2 * sqrtof2: Pa_done = True elif Tax_plus_3 == 3andnot Qa_done: Qa = Pa Qa_done = True else: raise ValueError('Unexcepted order 2 point.')
if Pa_done and Qa_done: break
assert Pa.order() == 2**a and Qa.order() == 2**a assert Pa.weil_pairing(Qa, 2**a) ** (2 ** (a - 1)) != 1
assert jab == jba h = bytes_to_long(hashlib.sha256(str(jab).encode()).digest()) enc = h ^ bytes_to_long(flag) print(f'enc = {enc}')
""" Pa = (199176096138773310217268*ii + 230014803812894614137371, 21529721453350773259901*ii + 106703903226801547853572) Qa = (8838268627404727894538*ii + 42671830598079803454272, 232086518469911650058383*ii + 166016721414371687782077) Pb = (200990566762780078867585*ii + 156748548599313956052974, 124844788269234253758677*ii + 161705339892396558058330) Qb = (39182754631675399496884*ii + 97444897625640048145787, 80099047967631528928295*ii + 178693902138964187125027) phib_Pa = (149703758091223422379828*ii + 52711226604051274601866, 112079580687990456923625*ii + 147229726400811363889895) phib_Qa = (181275595028116997198711*ii + 186563896197914896999639, 181395845909382894304538*ii + 69293294106635311075792) Ea: Elliptic Curve defined by y^2 = x^3 + (11731710804095179287932*ii+170364860453198752624563)*x^2 + x over Finite Field in ii of size 232900919541184113672191^2 Eb: Elliptic Curve defined by y^2 = x^3 + (191884939246592021710422*ii+96782382528277357218650)*x^2 + x over Finite Field in ii of size 232900919541184113672191^2 enc = 48739425383997297710665612312049549178322149326453305960348697253918290539788 """
The Flag
d3ctf{is0geny_gr4ph_m33t_1n_7he_m1ddl3}
Intended Solution
The main idea
Meet in the middle attack has time and space complexity $O(2^{a/2}) = O(p^{1/4})$. Using the given isogeny computing function, it takes about 1h to solve Alice’s secret isogeny $\phi_a$.
MITM
Let $K_a = 2^{a/2}P_a + s_a \cdot 2^{a/2}Q_a$, where $\left \langle P_a + s_a \cdot Q_a\right\rangle$ is the kernel of $\phi_a$. As a path in the 2-isogeny graph, $\phi_a$ passes through the node $E_{mid1} := E_0/ \left\langle{K_a}\right\rangle$.
Therefore we can compute all possible $2^{a/2}$-isogenies and construct a hash table. The keys are the j-invariant of codomains and the values are ${s_a}\pmod{2^{a/2}}$.
Remark: In my script, I compute $2^{a/2+1}$-isogenies and construct the hash table. Because I planned to set a=43 at first.
Use the j-invariant as the key of hash table
Note that each node of isogeny graph is an isomorphism class (over $\mathbb{F_{p^2}}$) of elliptic curve, so we need to compute the j-invariant (an isomorphic invariant) and use it as the key of our hash table.
The j-invariant of a untwisted Montgomery curve $E_A$ is
To avoid unnecessary computation of inversion, write $A = A_x/A_z$ and we get
See [1] for the j-invariant.
So we can use 1 inversion to compute the $j_A$ from the fractional representation of $A$.
Working with the given isogeny formula
Note that computing $2^a$-isogeny requires us to compute a chain a 2-isogeny and 4-isogeny. As SIKE’s spec pointed out in page 6~7, the input kernel point of isog2 cannot be $(0,0)$, and the input kernel point of isog4 cannot be $(1, \cdot)$ or $(-1, \cdot)$.
So how can we make use of this formula to compute isogenies starting from $E_a$? The generation process of $P_a$ and $Q_a$ in challenge.py suggests that we can choose $Q_a’$ such that $2^{a-1} \cdot Q_a’ = (0,0)$, and choose a $P_a’$ such that $e_{2^{a}}(P_a’, Q_a’)$ is a primitive $2^a$-th root.
Searching from $E_a$
According to the structure of 2-isogeny graph, there are $3 \cdot 2^{a/2-1} $ isogenies to try. But taking $2^{a/2}P_a’ + s_a’ \cdot 2^{a/2}Q_a’$ as the kernel only makes $2^{a/2}$ isogenies.
In fact, the other isogenies correspond to $\left\langle{t_a’P_a’ + Q_a’}\right\rangle$, where $t_a’=0,2, …, 2^{a/2}$.
After MITM
The collision gives us ${s_a}\pmod{2^a}$ and $t_a’, s_a’$ , where $E_{mid2} := E_a/\left\langle K_a’\right\rangle \cong E_{mid1}$ and $K_a’ = 2^{a/2}(t_a’P_a’ + s_a’Q_a’)$
Denote $\phi_1: E_0 \to E_{mid1}$, $\phi_2: E_a \to E_{mid2}$, $\sigma: E_{mid1} \xrightarrow{\sim} E_{mid2}$ . We can reconstruct the isogeny $\phi_a = \hat{\phi_2} \circ \sigma \circ \phi_1$ with the help of sagemath.
Since $\phi_a(P_a+s_aQ_a) = O$, we have $\phi_a(P_a) = -s_a\phi_a(Q_a)$. So we can solve the discrete logarithm to recover $s_a$
Now it is sufficient to follow SIDH protocol to find $j_{AB}$.
Some Further Optimization
In fact, there are some further optimization for the attack [3] and the computation of isogeny chain (optimal strategy) [2, 4].
# Get j-invariant from coeff (coeff in the form of fraction) defget_j_invariant(A: tuple): Ax, Az = A Ax_square = Ax**2 Az_square = Az**2 Jx = 256 * (Ax_square - 3 * Az_square) ** 3# Jx = 256(Ax^2 - 3Az^2)^3 Jz = Ax_square * Az_square**2 - 4 * Az_square**3# Jz = Ax^2*Az^4 - 4Az^6 if Jz == 0: print(f"A = {A}") return Jx * Jz ** (-1)
# phi: E -> E/<P> =: E2 # All points of input and output are projective x-coordinate, # i.e. if P = (XP: YP: ZP) ∈ E, then P is (XP, ZP) in the code # Input: E's A24, kernel point K of order 2, T the point to push # Output: E2's A24, phi(T) defisog2(K, T): if K[0] == 0: print("Meet the case we cannot use isog2.") returnNone
defxmul_2k(A24: tuple, k: int, T: tuple): assert k >= 0
defxdbl(P: tuple, A24: tuple) -> tuple: XP, ZP = P # assert XP != 0 and ZP != 0 if XP == 0or ZP == 0: print(f"k={k}") # print(f'XP={XP}') # print(f'ZP={ZP}') V1 = XP + ZP # line 1 of my pseudo code V1 **= 2 V2 = XP - ZP V2 **= 2 Z2P = A24[1] * V2 X2P = Z2P * V1 # line 6 of my pseudo code V1 -= V2 Z2P += A24[0] * V1 Z2P *= V1 return X2P, Z2P
for _ inrange(k): T = xdbl(T, A24) return T
# A: Estart affine coeff, T: kernel point (sagemath EC point object), k: T has order 2**k (same as the degree of this isogeny) defisog_2k(A, T, k): T = T[0], T[2] A24 = xA24(A) if k % 2 == 1: # T has order 2^a P = T P = xmul_2k(A24, k - 1, P) # now P has order 2 result_isog2 = isog2(P, T) if result_isog2 isNone: returnNone A24, T = result_isog2 k -= 1 assert k % 2 == 0 for i inrange(2, k + 2, 2): # T has order 2^(a-i+1) P = T P = xmul_2k(A24, k - i, P) # now P has order 4 result_isog4 = isog4(P, T) if result_isog4 isNone: returnNone A24, T = result_isog4
Pa_new = 2 ** (a - steps) * Pa Qa_new = 2 ** (a - steps) * Qa sa = sa_start
Ra = Pa_new + sa * Qa_new while sa < sa_stop: # Using the basis in SIKE spec, our formula never fail. Jnew, Anew = isog_2k(A0, Ra, steps) # table[Jnew] = (Anew, sa) table[Jnew] = sa # table[str(Jnew)] Ra += Qa_new sa += 1 return table
time1 = time.time()
steps = a // 2 + 1
print("Start computing table.")
# NOTE: The following code for parallel will raise segmentation fault on my PC... # Maybe the problem has something to do with thread lock of finite field element?
## multiproc # num_proc = 12 # # assert num_proc % 2 == 0 # size_per_proc = 2**steps // num_proc # parameters = [[i*size_per_proc, (i+1)*size_per_proc, steps] for i in range(0, num_proc)] # parameters[-1][1] = 2**steps # print(f'parameters = {parameters}') # with mp.Pool(num_proc) as p: # tables = p.starmap(compute_hash_table, parameters)
# NOTE: You need to convert finite field type to str, Integer type to int before json.dump # with open('hash_table.json', 'w') as f: # json.dump(table2)
### from another side (start from Ea)
Pa_E0 = deepcopy(Pa) Qa_E0 = deepcopy(Qa)
# First choose a nice basis for Ea[2^b] Pa = Ea(0) Qa = Ea(0) whileTrue: Qa = Ea.random_point() Qa = 3**b * Qa Ra = 2 ** (a - 1) * Qa if Ra.is_zero(): continue assert Ra.order() == 2 if Ra.xy()[0] == 0: break else: print(f"Ra = {Ra}") print("Qa found.") whileTrue: Pa = Ea.random_point() Pa = 3**b * Pa if Pa.weil_pairing(Qa, 2**a) ** (a - 1) != 1: break print("Pa found.")
assert Pa.order() == 2**a and Qa.order() == 2**a assert Pa.weil_pairing(Qa, 2**a) ** (2 ** (a - 1)) != 1
steps = a - steps
defsearch_from_Ea(steps=steps): # Save the collision with hashtable where the kernel is taPa + saQa collision_isogenies = [] # Save failed isogeny(ta, sa) where the kernel is taPa + saQa Pa_new = 2 ** (a - steps) * Pa Qa_new = 2 ** (a - steps) * Qa failed_isogenies = []
defsearch_one_isogeny(ta, sa, Ra): result_isog2k = isog_2k(ka, Ra, steps) if result_isog2k isNone: print(f"Isogeny failed for ta={ta}, sa={sa}") failed_isogenies.append((ta, sa)) return Jnew, Anew = result_isog2k if Jnew in table2: collision_isogenies.append((ta, sa, table2[Jnew])) print("Found one collision.") print(f"Start from Ea, sa={sa}, ta={ta}") print(f"Start from E0, sa={table2[Jnew]}") print(f"Meet at A={Anew}") return
# First try these 2**steps SIDH-like isogeny sa = 0 Ra = Pa_new while sa < 2**steps: search_one_isogeny(ta=1, sa=sa, Ra=Ra) Ra += Qa_new sa += 1 # Note that we also need to try these isogenies correspond to <(0, 1)>, <(2, 1)>, ..., <(2**steps, 1)> ta = 0 Ra = Qa_new two_Pa_new = 2 * Pa_new iflen(collision_isogenies) == 0: print("Collision not found in SIDH like isogenies") print("Start searching the second part!") while ta < 2**steps: search_one_isogeny(ta=ta, sa=1, Ra=Ra) Ra += two_Pa_new ta += 2 return collision_isogenies, failed_isogenies
## decryption h = bytes_to_long(hashlib.sha256(str(jab).encode()).digest()) flag = long_to_bytes(h ^ enc) print(f'flag = {flag}')
Unintended Solution
Actually, there are some alternative ways to find a collision.
In [5] there’s an attack leverages modular polynomial, and it compute the j-invariant only when computing 2-isogeny, instead of computing coefficient of the codomain. That’s really cool!
Another way is to get order-two points with sagemath’s division_points when going one step in DFS. (I thank hash_hash for telling me that) It seems that division_points make use of division polynomial to get these order-two points with low cost. Usually, when the degree of isogeny $l$ is a larger prime, we need to sample a random point $R$ and compute $(p+1)/l * R$. This is one of the main reasons why isogeny is slow! Another way is to push some points through the isogeny, but I think this is complicated in the scenario of MITM attack, perhaps not doable.
Here’s a MITM snippet using division_points and DFS from hash_hash.
from Crypto.Util.number import * defdp(way, E): iflen(way) > deepth: return0 ker = E(0).division_points(2)[1:] j_next = [E.isogeny_codomain(k).j_invariant() for k in ker] for j in j_next: if j notin forward.keys(): forward[j] = way+[j] dp(way+[j], EllipticCurve(j=Fp2(j)))
defdp_find(way, E): iflen(way) > deepth: return0 ker = E(0).division_points(2)[1:] j_next = [E.isogeny_codomain(k).j_invariant() for k in ker] for j in j_next: if j in forward.keys(): col = forward[j]+way[::-1] print("FIND!", len(col)) elif j notin back.keys(): back[j] = way+[j] dp_find(way+[j], EllipticCurve(j=Fp2(j)))
a = 38 b = 25 p = 2**a * 3**b - 1
assert is_prime(p) Fp = GF(p)
Fpx = PolynomialRing(Fp, "x") x = Fpx.gen() Fp2 = Fp.extension(x**2 + 1, "ii") ii = Fp2.gen()
By the way, a similar challenge had already appeared in the race SEETF 2023. Only after the race did I recognize that… See maple’s writeup for IsogenyMaze [6].
