Solution de n3ige86 pour Unknown Public Key

hardware attaque side channel

14 avril 2024

Origine

Une recherche perdue sur le net de nuit, et on tombe sur un article de 2004 Key Recovery Method for RSA CRT Implementation qui permet de retrouver une clef privée sans connaître la clef publique. Comme on a déjà une épreuve sans connaître la clef publique où il faut s’occuper de clef privée, autant en faire une où on demande la clef publique. De plus en lisant l’article, on constate une petite erreur dans l’algorithme. On tombe sur un dépôt Git implémentant l’article et il y a toujours l’erreur. Ce serait vicieux 😈 d’exploiter cette erreur dans une épreuve pour le FCSC…

Énoncé de l’épreuve

Une histoire à dormir debout sert de prétexte pour justifier qu’on recherche une clef publique qui n’est pas publique. A-t-on idée? On nous parle également d’une exécution en temps non constant.

Le matériel fourni pour cette épreuve est minimaliste. Un seul couple (message,signature) et une seule trace d’exécution 😨. Ça va se jouer entre cette trace et nous !

(Il s’agit de la seule épreuve du FCSC2024 n’ayant pas été résolue.)

Analyse

On commence par afficher cette fameuse trace à l’aide par exemple de matplotlib.

import numpy as np
import matplotlib.pyplot as plt
trace = np.load("trace_file.npy")
plt.plot(trace)
plt.show()

trace

On aperçoit distinctement 2 gros pâtés, ce qui fait penser à un code exploitant le théorème des restes chinois pour calculer la signature RSA: un calcul de $S_q=m^{d_q}\mod q$ et un calcul de $S_p=m^{d_p}\mod p$.

Si on zoom sur le début de l’exponentiation, on a l’image suivante:

debut_exponentiation

On peut voir qu’on monte en puissance, avant d’arriver sur une série stable à partir des environs du point (85000,50). Cette (longue) série stable représente des multiplications modulaires, dont on n’arrive pas vraiment à distinguer s’il s’agit de carré modulaire ou de simple multiplication modulaire. Pour expliquer ce qui précède, il faut savoir qu’avant d’effectuer l’exponentiation, il faut préparer le calcul modulaire, et également faire la réduction modulaire $m\mod q$. En effet, le message $m$ est de la taille du module public $n$ alors qu’on veut calculer modulo $q$ ou $p$ pour des raisons de performance. On en profite d’ailleurs pour deviner la taille de $n$ à partir du fichier output.txt. Ce fichier contient deux grands nombres message et signature de tailles respectivement 1020 et 1023 bits. Selon toute vraisemblance, on a donc $2^{1023} < n <2^{1024}$.

Il existe plusieurs algorithmes pour effectuer une exponentiation. L’énoncé indique toutefois que l’algorithme ne s’effectue pas en temps constant, il s’agit potentiellement d’un algorithme simple. Commençons par cela donc:

#exponentiation from most significant bit
def exp_msb(m,e,q):
    r = 1
    exp = bin(e)[2:]
    for i in exp:
        r = gmpy2.mod(r*r,q)
        if i=='1'
            r = gmpy2.mod(r*m,q)
    return r
#exponentiation from least significant bit
def exp_lsb(m,e,q):
    r = 1
    a = m
    exp = bin(e)[2:]
    for i in exp[::-1]:
        if i=='1'
            r = gmpy2.mod(r*a,q)
        a = gmpy2.mod(a*a,q)
    return r

Sur de tels algorithmes, il faut un certain temps pour évaluer le if, ainsi que pour boucler, ce qui donne une légère différence de temps d’exécution entre les multiplications. Par exemple pour l’algorithme exp_msb, sur un tour de boucle, on a les séquences suivantes:

  • si exp[i]=='0': carré modulaire, test if, boucle suivante, carré modulaire
  • si exp[i]=='1': carré modulaire, test if, multiplication modulaire, boucle suivante

Pour l’algorithme exp_lsb, sur un enchaînement de boucle, on a les séquences suivantes:

  • si exp[i]=='0': carré modulaire, boucle suivante, test if, carré modulaire
  • si exp[i]=='1': carré modulaire, boucle suivante, test if, multiplication modulaire, carré modulaire

On peut noter une différence de séquence entre les deux parcours d’exposant. Dans un cas il y a toujours du temps entre deux multiplications (test if ou test de boucle), dans l’autre cas il arrive que deux multiplications s’enchaînent très rapidement. Puisqu’on observe une nette spération entre les multiplications, on peut supposer que nous sommes sur un parcours partant des poids fort en premier.

Il s’agit donc de déterminer la durée entre 2 multiplications, afin de déterminer si on est en présence d’un carré (cas test if+boucle) ou d’une multiplication (cas test if ou boucle, suivant si on regarde la durée avant la multiplication ou après). Si l’on zoom encore plus sur les multiplications, on s’intéresse donc à la durée indiquée par le cercle rouge de l’image suivante: timing_entre_multiplication

Pour cela, on peut récupérer les abscisses de tous les points au-dessus de l’ordonnée 50 (j’ai utilisé 51.5 pour ma part) grâce à la fonction find_peaks:

from scipy.signal import find_peaks
trace = np.load("../public/trace.npy")
trace_q = trace[ :5633000].astype(np.int16)
trace_p = trace[ 5633000:].astype(np.int16)

plt.plot(trace)
plt.show()

l = find_peaks(trace_q, height=51.5, distance = 340)[0] #340 indique une plage pour laquelle on ignore les valeurs
d = l[1:]-l[:-1]           #d[i] = l[i+1] - l[i]
plt.scatter(range(len(d)),d)
plt.show()

l = find_peaks(trace_p, height=51.5, distance = 340)[0]
d = l[1:]-l[:-1]            #d[i] = l[i+1] - l[i]
plt.scatter(range(len(d)),d)
plt.show()

On obtient alors la figure suivante pour le premier exposant: lecture_exposant

On voit clairement 3 bandes. La bande du dessous correspond aux durées entre les pics au sein d’une multiplication, la bande du dessus correspond à la durée entre 2 multiplications. Il est plus difficile d’identifier la bande du milieu, mais elle ressemble à la bande du dessus puisqu’elle contient également 3 valeurs possibles.

On a maintenant le matériel pour retrouver des bits des exposants privés, encore faut-il deviner l’algorithme utilisé.

Mais même en admettant qu’on retrouve les exposants privés, comment obtenir le reste de la clef? Sans la connaissance de l’article à l’origine de l’épreuve, ce n’est pas si facile.

Intuitions nécessaires

Il faut deviner l’algorithme utilisé pour calculer les exponentiations. Pour cela, il faut une bonne connaissance des implémentations bas niveau, et tester les cas possibles jusqu’à tomber sur celui qui permet de finir l’attaque. Suivant l’algorithme considéré, on pourra associer les mesures de temps observés entre les multiplications à des bits d’exposant.

Il faut avoir l’intuition de l’article, ou tomber sur l’article via une recherche internet, à savoir partir sur une recherche exhaustive des valeurs possibles de l’exposant public $e$, en supposant que cette valeur soit suffisamment petite pour être trouvée en temps raisonnable. À l’aide de la connaissance de $e$, on peut en effet faire une recherche exhaustive pour retrouver $p$ et $q$ à l’aide des exposants privés $d_p$ et $d_q$. Cette banque étrangère et bien étrange avait évité de prendre un $e$ classique.

Pour finir il faut trouver une explication aux tailles des exposants privés obtenus. En effet, les exposants $d_p$ et $d_q$ lus à partir de la trace font chacun 528 bits. Comme le module public $n$ a une taille de 1024 bits, les premiers $p$ et $q$ qui le composent ont une taille de 512 bits chacun. Les exposants privés sont donc étendus de 16 bits chacun, sans que cela ne produise une signature erronnée (on peut le supposer, même si on n’a pas la clef publique). C’est classiquement ce qui est fait par le masquage, où un multiple aléatoire de l’ordre est ajouté à chaque calcul de signature:

  • $d_p + m_p\cdot (p-1)$
  • $d_q + m_q\cdot (q-1)$

Si c’est le cas, il nous faudra parcourir exhaustivement ces valeurs de masque.

Solution

L’algorithme utilisé en pratique fait l’hypothèse d’un exposant non nul et initialise le résultat à $m$.

#exponentiation from most significant bit
def exp_msb(m,e,q):
    r = m
    exp = bin(e)[3:]
    for i in exp:
        r = gmpy2.mod(r*r,q)
        if i=='1'
            r = gmpy2.mod(r*m,q)
    return r

Il faut donc prendre cela en compte dans la reconstruction des exposants privés:

def build_exponent(d):
    res = 1
    for i in d:
        if i == 935 or i == 542: #select a timing, between square and multiplication or between multiplication and square
            res <<= 1
            res += 1
        if i > 935 or i == 550: #highest timings correspond to squares
            res <<= 1
    return res

Par définition, on a $e\cdot d_p=1+k_p\cdot (p-1)$, donc $e\cdot (d_p+m_p\cdot(p-1))=1+(k_p+e\cdot m_p)\cdot(p-1)$ et alors:

$p=1+\frac{e\cdot (d_p+m_p\cdot(p-1))-1}{k_p+e\cdot m_p}$

avec $1 < k_p < e$ et $2^{15}-1 < m_p < 2^{16}$. On a les mêmes calculs pour retrouver $q$ avec $k_q$, $m_q$ et $d_q$.

On en arrive enfin à l’erreur de l’article, à savoir oublier le cas $k_p=2$ ou $k_q=2$. Modulo cette petite erreur, et la recherche des valeurs de masques, la résolution est identique à l’article. Le message et la signature donnés permettent de distinguer la solution attendue parmi les solutions possibles.

def RSATestKey(p,q,e,m,s):
    return m == gmpy2.powmod(s,e,p*q)
#dp and dq acutally denote masked dp and dq
def RSAKeyRecoveryMaskedExponent(dp, dq, msg, sig, maxe):
    for e in range(3,maxe+1,2):
        kp = e*dp-1         #kp = k.(p-1) for unknown k if the guess on e is correct
        for k in range(2,e):
            for m in range(2**16):
                quotient, remainder = gmpy2.t_divmod(kp, k+m*e)
                p = quotient+1
                if 0 == remainder and gmpy2.is_prime(p):
                    lq = e*dq-1
                    for l in range(2,e):
                        for r in range(2**16):
                            quotient, remainder = gmpy2.t_divmod(lq, l+r*e)
                            q = quotient+1
                            if 0 == remainder and gmpy2.is_prime(q):
                                if RSATestKey(p,q,e,msg,sig):
                                    return p,q,e
    return 1,1,1

Les valeurs sont:

$$ d_q = \begin{array}{c} \mathtt{0xa37ee92c376626d6e970f150214ab70bd051fa5fa2bd034795e33582d610feb75b} \\ \mathtt{c16bc6d66c7111a67af9876d4c4c676fa9b275e305d22b75ccedc699c5ece53f49} \end{array} $$

$$ d_p = \begin{array}{c} \mathtt{0x98aa6b5995e6a3590064bbe414f533fbbe4d882d2fbfb6bceacc3f2f8f00462ec4} \\ \mathtt{dd248c176b1642eacaba112b9b4b663fb1f5d3ef7e9172ab42c58bead025253ced} \end{array} $$

$$ n = \begin{array}{c} \mathtt{0xed124c18e32d106c7fdb6a12aff7d1ab6a9ca7a649efb256cfe1eb85f5258a93} \\ \mathtt{5b3e067ff1fd8f4b9632f8cf76fe6486114d0ebc13b0296b683095aa91cb3f5a} \\ \mathtt{c542e1e23eae7e2d29915b977683a5bc967fc7dd8d3fae1b2a7983115974b453} \\ \mathtt{efb169985318d8776440b3d88d308b98b7d65d7fa63b34042415aa499328c8d1} \end{array} $$

$$ e = 255 $$

Anecdotes de conception

Il s’agit d’une trace simulée avec QEMU d’un code personnel de librairie des grands nombres. Le code original est en C, mais la non-gestion des retenues au niveau du C impose des tests conditionnels pour chaque multiplication ce qui empêchait de retrouver les exposants via une analyse temporelle sur une seule trace. Les calculs ont donc été faits en assembleur sur cortex-M3 afin de corriger ce problème de retenue.

Pour ceux que ça intéresse, au début de la trace, on voit les calculs suivants:

  • Calcul de $v = -q^{-1} \mod 2^{32}$
  • Calcul de $R = 2^{512}\mod q$
  • Calcul de $R^2 \mod q$ avec 5 multiplications de Montgomery
  • Réduction modulaire de $m$: $m\star 1$ où $\star$ désigne la multiplication de Montgomery
  • Calcul de $(m\star 1)\star R^2 \star R^2 = m\cdot R\mod q$

init_expo