Supersingular isogeny key exchange for beginners

Supersingular isogeny key exchange for beginners

Companion notebook to Craig Costello’s article
by Wojciech Nawrocki

This SageMath notebook is about RSA. Sike! It’s actually about Diffie-Hellman key exchange using isogenies of supersingular elliptic curves. Being a companion to Craig Costello’s article, the notebook will not make any sense without it. A reasonable way to read it would be to, for each chapter, first read Costello’s text and then the corresponding code here. As with any bit of math, the only way to learn this is to play around with it. So download the notebook, run it, mess around with the values, try out other functionality, etc.

For API reference, see Sage docs: Elliptic curves and the industrial-strength spec.

2 The set of supersingular \(j\)-invariants

Sage docs: Elliptic curves over finite fields

# Choose a prime. In real crypto, this prime would be 256-bit or so.
p = 431
# It will come up later why p is of this form.
assert p == 2^4*3^3 - 1 and is_prime(p)

# We work in a field F - the quadratic extension of F_p.
# Elements of F are given by `u + v*i` with u,v ∈ F_p.
# (Using `i` overrides Sage's notation for the imaginary unit,
# but we don't need it anyway.)
F.<i> = GF(p^2, modulus=x^2+1)
Finite Field in i of size 431^2

\(J\)-invariants of supersingular elliptic curves are in one-to-one correspondence with isomorphism classes of such curves. \(E_{a_1}\) and \(E_{a_2}\) have the same \(j\)-invariant, so they are isomorphic.

def E_a(a):
    """Returns a curve given by Montgomery form y² = x³ + a*x² + x."""
    return EllipticCurve(F, [0, a, 0, 1, 0])

E_a1 = E_a(208*i + 161)
E_a2 = E_a(172*i + 162)
assert E_a1.is_supersingular() and E_a2.is_supersingular()

assert E_a1.j_invariant() == E_a2.j_invariant()
assert E_a1.is_isomorphic(E_a2)
364*i + 304

3 Isogenies

Sage docs: Isogenies

The doubling and tripling maps are separable isogenies from \(E_{a_1}\) to itself. The former has \(4\) elements in the kernel, while the latter has \(9\).

doublingMap = E_a1.multiplication_by_m_isogeny(2)
assert == 4

triplingMap = E_a1.multiplication_by_m_isogeny(3)
assert == 9

# Their kernels are the 2- and 3-torsions, respectively.
for P in E_a1(0).division_points(2):
    assert doublingMap(P) == E_a1(0)

for P in E_a1(0).division_points(3):
    assert triplingMap(P) == E_a1(0)

every subgroup \(G\) of points on an elliptic curve \(E\) gives rise to a unique isogeny \(\phi: E \to E'\) whose kernel is \(G\), and vice versa

But only up to isomorphism. Computing an isogeny from the doubling map’s kernel, we actually get a different doubling map whose codomain is isomorphic to \(E_{a_1}\) but is not \(E_{a_1}\).

G = E_a1(0).division_points(2)
ϕ = E_a1.isogeny(G)
assert ϕ.codomain().is_isomorphic(E_a1)
assert ϕ.codomain() != E_a1

# We can transport points along the isomorphism
# to recover the doubling map. Let's test this on
# some points of order 3.
iso = ϕ.codomain().isomorphism_to(E_a1)
for P in E_a1(0).division_points(3):
    assert iso(ϕ(P)) == doublingMap(P)

Now make an isogeny \(\phi: E \to E/G\) with the cyclic subgroup \(G = \{\mathcal{O}, (\alpha,0)\}\) of order \(2\) as its kernel.

α = E_a1(350*i + 68, 0)
assert 2*α == E_a1(0)

# It suffices to specify α here since 𝒪 must be in the kernel anyway.
ϕ = E_a1.isogeny(α)
# Its codomain is isomorphic to the one in the text.
assert ϕ.codomain().is_isomorphic(E_a(102*i + 423))
# Curves connected by an isogeny are isogenous.
assert E_a1.is_isogenous(ϕ.codomain())

Every isogeny \(\phi: E \to E'\) of degree \(d\) has a dual isogeny \(\phi': E' \to E\) s.t. \(\phi' \circ \phi = [d]\), the multiply-by-\(d\) map on \(E\), and \(\phi \circ \phi' = [d]\) on \(E'\). Composing the isogeny with kernel \(G\) with its dual gives the doubling map.

ϕDual = ϕ.dual()
for P in E_a1(0).division_points(2):
    assert (ϕDual * ϕ)(P) == E_a1(0)
for P in ϕ.codomain()(0).division_points(2):
    assert (ϕ * ϕDual)(P) == ϕ.codomain()(0)

The point \(P\) below has order \(8\), but is s.t. \([4]P = (\alpha,0)\). So \(\phi([4]P) = \mathcal{O}\). Since isogenies are group homomorphisms, \(\phi([4]P) = [4]\phi(P) = \mathcal{O}\). Thus, the order of \(P\) drops to \(4\) on passing through \(\phi\). \(Q\) is just an ordinary point without this property.

P = E_a1(390*i + 23, 104*i + 7)
Q = E_a1(151*i + 140, 110*i + 136)
assert P.order() == Q.order() == 8
assert ϕ(P).order() == 4
assert ϕ(Q).order() == 8

Two curves over a finite field \(\mathbb{F}_q\) are isogenous over it iff they have the same number of points over \(\mathbb{F}_q\). Additionally, in the case here the group \(E(\mathbb{F}_{p^2})\) is isomorphic to \(\mathbb{Z}_{p+1} \times \mathbb{Z}_{p+1}\).

assert E_a1.cardinality() == ϕ.codomain().cardinality()
assert E_a1.cardinality() == (p+1)^2

4 Isogeny graphs

First set up some graph-drawing code with cytoscape.js since matplotlib looks awful.

import ipycytoscape
import networkx as nx

def drawGraph(edges):
    """Draws a graph using cytoscape.js and networkx."""
    edges = { str(j) : { str(e) for e in es } for (j,es) in edges.items() }
    Gnx = nx.Graph(edges)
    gWidget = ipycytoscape.CytoscapeWidget()
        'selector': 'node',
        'css': {
        'content': 'data(label)',
        'text-valign': 'center',
        'color': 'white',
        'text-outline-width': 2,
        'text-outline-color': 'blue',
        'background-color': 'blue' }}])
    #Uncomment for bouncy graph
    #gWidget.set_layout(animate=True, infinite=True, fit=False)
    for v, dat in Gnx.nodes(data=True):
        dat['label'] = v
    return gWidget

Now let’s reproduce the \(2\)- and \(3\)-isogeny graphs!

def mkLIsogenyGraph(l):
    # Start at E_a1
    curvesToAdd = [E_a1]
    j = E_a1.j_invariant()
    jInvariants = { j: set() }

    while len(curvesToAdd) != 0:
        E = curvesToAdd.pop()
        for P in E(0).division_points(l):
            if P != E(0):
                ϕ = E.isogeny(P)
                j = ϕ.codomain().j_invariant()
                if not j in jInvariants:
                    jInvariants[j] = set()

    return jInvariants


5 Walking through the protocol

Recall that \(p=2^{e_A}3^{e_B}-1\) with \(e_A=4\), \(e_B=3\). Let’s confirm that the \(2^{4}\)-torsion \(E_{a_1}[2^4] \cong \mathbb{Z}_{2^4} \times \mathbb{Z}_{2^4}\) has \((2^4)^2\) points and that \(E_{a_1}[3^3]\) has \((3^3)^2\).

assert len(E_a1(0).division_points(2^4)) == (2^4)^2
assert len(E_a1(0).division_points(3^3)) == (3^3)^2

Now let’s do a key exchange. Fix the public protocol parameters - the starting curve \(E\) as well as Alice’s basis points \(P_A,Q_A\) and Bob’s basis \(P_B,Q_B\).

E = E_a(329*i + 423)
P_A = E(100*i + 248, 304*i + 199)
Q_A = E(426*i + 394, 51*i + 79)
P_B = E(358*i + 275, 410*i + 104)
Q_B = E(20*i + 185, 281*i + 239)

Confirm that \(E\) is supersingular and that the basis points are of the right order.

assert E.is_supersingular()
assert P_A.order() == Q_A.order() == 2^4
assert P_B.order() == Q_B.order() == 3^3

Alice generates her secret value \(k_A\) and subgroup generator \(S_A=P_A+[k_A]Q_A\).

k_A = 11 # Chosen by fair dice roll
S_A = P_A + k_A*Q_A
# S_A remains in E[2^4]
assert S_A.order() == 2^4

She then computes her secret \(2^4\)-isogeny \(\phi_A: E \to E/\langle S_A \rangle\). In Sage, we can actually just do this directly using Vélu’s formulas.

ϕ_A_direct = E.isogeny(S_A)
assert ϕ_A_direct.codomain().j_invariant() == 222*i + 118
assert ϕ == 2^4

[But since] algorithms for general isogeny computation are linear in the degree of the isogeny, this would be out of the question if it was not for our being able to instead compute it as the composition of \(e\) individual \(2\)-isogenies.

So instead we compute the secret isogenies in \(e_A\) \(2\)-isogeny steps for Alice and \(e_B\) \(3\)-isogeny steps for Bob. Let’s not bother with unrolling Alice and Bob’s key generation loops like the text does and instead define a generic function to compute an \(l^e\)-isogeny \(E \to E/\langle S \rangle\) as a composition of \(e\) \(l\)-isogenies.

def computeIsogenies(E, S, l, e):
    """Returns a list of e l-isogenies making up the target isogeny."""
    S_tmp = S
    E_tmp = E
    ϕs = []
    for k in range(e):
        # Generate a point in the l-torsion via repeated
        # doubling (l=2) or tripling (l=3) starting from S_tmp.
        R_tmp = S_tmp
        for _ in range(e-k-1):
            R_tmp = l*R_tmp
        assert R_tmp.order() == l
        # Generate an l-isogeny to take a step in the graph.
        ϕ_k = E_tmp.isogeny(R_tmp)
        assert ϕ == l
        S_tmp = ϕ_k(S_tmp)
        # Each S_tmp lies above a point in the l-isogeny's kernel,
        # so its order decreases by a factor of l per step.
        assert S_tmp.order() == l^(e-k-1)
        E_tmp = ϕ_k.codomain()
    return ϕs

Why does this algorithm work in the first place? Consider Alice’s example and the fact that, for \(\phi: E \to E/\langle P \rangle\), \(E/\langle P,Q \rangle \cong (E/\langle P \rangle)/\langle \phi(Q) \rangle \hspace{4em} (1)\)

We want to get to \(E/\langle S_A\rangle\). Since \(S_A\) has order \(2^4\), \(\langle S_A \rangle = \langle 8S_A, 4S_A, 2S_A, S_A \rangle\). Using a subset of this and fact \((1)\), our target is \(E/\langle S_A \rangle = E/\langle 2S_A, S_A \rangle \cong (E/\langle 2S_A \rangle)/\langle \phi'(S_A) \rangle\) for \(\phi': E \to E/\langle 2S_A \rangle\).

But then we can iterate the expansion because \(\langle 2S_A \rangle = \langle 4S_A, 2S_A \rangle\) and so on up to \(8S_A\). This gives \(E/\langle S_A \rangle \cong (((E/\langle 8S_A \rangle)/\langle \phi'''(4S_A) \rangle)/\langle \phi''(2S_A) \rangle)/\langle \phi'(S_A) \rangle\) for \(\phi''': E \to E/\langle 8S_A \rangle\) and \(\phi'': E \to E/\langle 4S_A \rangle\).

The above sequence of isogenies almost looks like the loop we have. Looking at the algorithm we can identify \(\phi'''=\phi_0\) which is already a \(2\)-isogeny and hence fast to compute, but we need to break the others into sub-steps. In the end we get \(E/\langle S_A \rangle \cong (((E/\langle 8S_A \rangle)/\langle \phi_0(4S_A) \rangle)/\langle \phi_1(\phi_0((2S_A)) \rangle)/\langle \phi_2(\phi_1(\phi_0(S_A))) \rangle\) where:

  • \(\phi_1: E/\langle 8S_A \rangle \to E/\langle 4S_A \rangle\)
  • \(\phi_2: E/\langle 4S_A \rangle \to E/\langle 2S_A \rangle\)
  • \(\phi_3: E/\langle 2S_A \rangle \to E/\langle S_A \rangle\)

and the secret isogeny is \(\phi_A = \phi_3 \circ \phi_2 \circ \phi_1 \circ \phi_0\).

Now Alice can generate her information to be exchanged publicly.

ϕs_A = computeIsogenies(E, S_A, 2, 4)
# Check that we took the expected path.
assert [ ϕ.codomain().j_invariant() for ϕ in ϕs_A ] == \
       [ 107, 344*i + 190, 350*i + 65, 222*i + 118 ]
# As of 9.1, Sage doesn't implement isogeny composition so the best
# we can get is a "composite map".
ϕ_A = ϕs_A[3] * ϕs_A[2] * ϕs_A[1] * ϕs_A[0]

# Alice's public key.
# We can't just store a₄ here like Costello does
# because our curves aren't in Montgomery form.
# Note also that we didn't carry the basis of Bob's torsion subgroup
# through the isogeny in the generic function, so we do that here.
Alice = (ϕ_A.codomain(), ϕ_A(P_B), ϕ_A(Q_B))

Bob does the same.

k_B = 2 # Bob's secret dice roll
S_B = P_B + k_B*Q_B
assert S_B.order() == 3^3

ϕs_B = computeIsogenies(E, S_B, 3, 3)
assert [ ϕ.codomain().j_invariant() for ϕ in ϕs_B ] == \
       [ 106*i + 379, 325*i + 379, 344*i + 190 ]
ϕ_B = ϕs_B[2] * ϕs_B[1] * ϕs_B[0]
Bob = (ϕ_B.codomain(), ϕ_B(P_A), ϕ_B(Q_A))

Finally they both compute the shared secret \(j\)-invariant by starting from the other party’s curve, shifting their own secret generator to lie on said curve, and computing the same sequence of small isogenies as before. This time on Bob’s side for variety, we start from the curve \(E/\langle S_A \rangle\) and use the secret generator \(_AS_B = \phi_A(S_B) = \phi_A(P_B) + [k_B]\phi_A(Q_B)\) to get to \(E/\langle S_A, S_B \rangle\). The validity of this is justifiable using fact \((1)\) from before.

# Alice
(E_B, B_P_A, B_Q_A) = Bob
B_S_A = B_P_A + k_A*B_Q_A
ϕs_A = computeIsogenies(E_B, B_S_A, 2, 4)
assert [ ϕ.codomain().j_invariant() for ϕ in ϕs_A ] == \
       [ 364*i + 304, 67, 242, 234 ]
jAlice = ϕs_A[-1].codomain().j_invariant()

# Bob
(E_A, A_P_B, A_Q_B) = Alice
A_S_B = A_P_B + k_B*A_Q_B
ϕs_B = computeIsogenies(E_A, A_S_B, 3, 3)
assert [ ϕ.codomain().j_invariant() for ϕ in ϕs_B ] == \
       [ 299*i + 315, 61, 234 ]
jBob = ϕs_B[-1].codomain().j_invariant()

# Great success!
assert jAlice == jBob == 234

6 SIDH in practice

Since Sage does most of the heavy lifting, we can carry out a full key exchange with cryptographically-sized parameters in a rather short snippet. Let’s instantiate public parameters from the SIKEp434 specification.

p = 2^216*3^137 - 1
F.<i> = GF(p^2, modulus=x^2+1)
E = EllipticCurve(F, [0, 6, 0, 1, 0])


Q_A = E(xQ20 + xQ21*i, yQ20 + yQ21*i)
P_A = E(xP20 + xP21*i, yP20 + yP21*i)
Q_B = E(xQ30 + xQ31*i, yQ30 + yQ31*i)
P_B = E(xP30 + xP31*i, yP30 + yP31*i)

And then carry out the computations! Unlike an optimised algorithm which would work almost instantaneously, this takes a few seconds, but that’s still pretty good for what’s mostly Python.

def computeIsogeny(E, S, l, e):
    S_tmp = S
    E_tmp = E
    ϕ = None
    for k in range(e):
        R_tmp = S_tmp
        for _ in range(e-k-1):
            R_tmp = l*R_tmp
        ϕ_k = E_tmp.isogeny(R_tmp)
        S_tmp = ϕ_k(S_tmp)
        E_tmp = ϕ_k.codomain()
        if ϕ is None:
            ϕ = ϕ_k
            ϕ = ϕ_k * ϕ
    return ϕ

# Alice generates public information
k_A = randint(0, 2^216-1)
S_A = P_A + k_A*Q_A
ϕ_A = computeIsogeny(E, S_A, 2, 216)
Alice = (ϕ_A.codomain(), ϕ_A(P_B), ϕ_A(Q_B))

# Bob
k_B = randint(0, 3^137-1) 
S_B = P_B + k_B*Q_B
ϕ_B = computeIsogeny(E, S_B, 3, 137)
Bob = (ϕ_B.codomain(), ϕ_B(P_A), ϕ_B(Q_A))

# Alice computes shared secret
(E_B, B_P_A, B_Q_A) = Bob
B_S_A = B_P_A + k_A*B_Q_A
jAlice = computeIsogeny(E_B, B_S_A, 2, 216).codomain().j_invariant()

# Bob
(E_A, A_P_B, A_Q_B) = Alice
A_S_B = A_P_B + k_B*A_Q_B
jBob = computeIsogeny(E_A, A_S_B, 3, 137).codomain().j_invariant()

assert jAlice == jBob

7 Security and cryptanalysis

Sage provides an algorithm for computing the isogeny from a domain and codomain. The hardness assumption at the heart of SIDH is that any such algorithm takes infeasibly long to finish on parameters of cryptographic size.

def thisWontWork():
    ϕ_A_direct = E.isogeny(kernel=None, codomain=ϕ_A.codomain(), degree=2^216)

7.1 Attack on static keys

Unlike SIKE, the original SIDH cryptoscheme is vulnerable to leakage of static secret keys. If Alice reuses her secret key \(k_A\) for multiple key exchanges, a malicious Bob can leak it bit-by-bit by sending \(log_2(k_A)\) specially crafted requests.

Reference: On the Security of Supersingular Isogeny Cryptosystems, ch. 3

The idea is that we can use the success or failure of key exchange as an oracle which determines whether Alice recovered the correct curve \(E_{AB} = E/\langle S_A, S_B \rangle\) or not. From the isogeny-subgroup bijection we know she will do so if and only if she computes an isogeny from \(E_B=E/\langle S_B \rangle\) with the correct subgroup \(\langle _BS_A \rangle = \langle \phi_B(P_A) + [k_A]\phi_B(Q_A) \rangle\) as kernel. But the curve points \(\phi_B(P_A)\) and \(\phi_B(Q_A)\) are sent by Bob, and as it turns out there is no way for Alice to reliably verify whether they are legit. Consequently, Bob can manipulate them for nefarious purposes.

For brevity, let’s use notation from On the Security... Let \(R \triangleq \phi_B(P_A)\) and \(S \triangleq \phi_B(Q_A)\). Let Alice’s private key be \(\alpha \triangleq k_A\) and her subgroup order exponent be \(n \triangleq e_A\). We can expand the key in binary as \(\alpha = \alpha_0 + 2^1\alpha_1 + 2^2\alpha_2 + \ldots + 2^{n-1}\alpha_{n-1}\).

Now suppose we have already leaked \(i\) bits and denote that part of the key by \(K_i\). Initially we know nothing, so \(K_0 = 0\). When we are interested in the \(i\)-th bit, we write the key as \(\alpha = K_i + 2^i\alpha_i + 2^{i+1}\alpha'\) with \(\alpha'\) denoting all higher bits.

Now consider what happens when Bob sends the following points:

\[(R - [2^{n-i-1}K_i]S, [1 + 2^{n-i-1}]S)\]

Thinking that these are her basis on \(E_B\), Alice will construct the subgroup

\[\begin{aligned} &\langle R - [2^{n-i-1}K_i]S + [\alpha][1 + 2^{n-i-1}]S \rangle \\ = &\langle R + [\alpha]S + [(\alpha-K_i)2^{n-i-1}]S \rangle \\ = &\langle R + [\alpha]S + [(2^i\alpha_i + 2^{i+1}\alpha')2^{n-i-1}]S \rangle \\ = &\langle R + [\alpha]S + [2^{n-1}\alpha_i]S + [2^n\alpha']S \rangle \\ = &\langle R + [\alpha]S + [2^{n-1}\alpha_i]S \rangle \\ \stackrel{?}{=} &\langle R + [\alpha]S \rangle \end{aligned}\]

Where \([2^n\alpha']S\) vanishes because \(S\) is in the \(2^n\)-torsion. Now some lemmas from the paper tell us that the last question-equality \(\stackrel{?}{=}\) holds if and only if \(\alpha_i = 0\).

So, when that bit of Alice’s secret key is \(0\), key exchange succeeds and we get the same \(j\)-invariant, while when it’s \(1\) it fails! By sending this construction for all \(i \in [0, n)\) we can recover the key. Read the paper for additional details about how to avoid certain forms of detection by bruteforcing the last two or so bits.

We’re approaching the limits of Sage’s speed, so I only run the attack below for the first 11 bits.

def simulateAlice(E_B, B_P_A, B_Q_A):
    """Pretends to be Alice computing an isogeny from a subgroup
    on E_B and returning the resulting j-invariant."""
    B_S_A = B_P_A + k_A*B_Q_A
    return computeIsogeny(E_B, B_S_A, 2, 216).codomain().j_invariant()

(E_B, R, S) = Bob
Ki = 0
n = 216
bitsToLeak = 11
for i in range(bitsToLeak):
    b = -Ki*2^(n-i-1)
    d = 1 + 2^(n-i-1)
    EvilBob = (E_B, R + b*S, d*S)
    if simulateAlice(*EvilBob) != jBob:
        Ki += 2^i
    print(f"k_A = {bin(k_A)[-bitsToLeak:]}")
    print(f"Ki  = {bin(Ki)[2:].rjust(bitsToLeak,'0')}")

assert Ki & (2^bitsToLeak - 1) == k_A & (2^bitsToLeak -1)
k_A = 11101000001
Ki  = 00000000001
k_A = 11101000001
Ki  = 00000000001
k_A = 11101000001
Ki  = 00000000001
k_A = 11101000001
Ki  = 00000000001
k_A = 11101000001
Ki  = 00000000001
k_A = 11101000001
Ki  = 00000000001
k_A = 11101000001
Ki  = 00001000001
k_A = 11101000001
Ki  = 00001000001
k_A = 11101000001
Ki  = 00101000001
k_A = 11101000001
Ki  = 01101000001
k_A = 11101000001
Ki  = 11101000001

7.2 Supersingular Isogeny Key Encapsulation

SIKE doesn’t really solve the leakage problem. Instead, it requires that Alice always learns Bob’s secret key \(k_B\). So while Bob cannot carry out the above attack anymore, the protocol becomes asymmetric in that Alice’s secret key can safely be static, but Bob’s must be ephemeral.

# TODO encapsulate key