Broadcasting NTRU writeup ~NTRU暗号の理論を添えて~

:::note info
こちらは[脆弱エンジニアの Advent Calendar 2024](https://qiita.com/advent-calendar/2024/full-weak-engineer)の14日目の記事です。
:::

こんにちは。
ぽちゃま[@PotyaExe](https://x.com/PotyaExe)と申します。
現在、修士2年で博士課程への進学を予定しています。ただし、受験予定の大学の募集要項の関係で前期日程での受験ができないため、まだ確定していません。

さて、[脆弱さん](https://x.com/asusn_online)にCODE BLUEのアフターパーティーでお会いできなかったため、こちらに参加させていただきました。
今後お会いする機会がありましたら、どうぞよろしくお願いいたします。

現在、CTFで主にCrypto分野の問題を解いており、絶賛Cryptoについて勉強中です。
そんな中、AlpacaHackでNTRU暗号に関する問題が出題されました。

高専在学時代に量子暗号をゼミで学んでいたこともあり、ぜひ再勉強を兼ねて解いてみたいと思います。
なお、今回の問題の作者である[チョコラスクさん](https://x.com/nuo_chocorusk)が作問者Writeupも公開されていますので、ぜひご覧ください。

https://chocorusk.hatenablog.com/entry/2024/10/12/180717
# 今回解く問題

https://alpacahack.com/challenges/broadcasting-ntru

# 1. NTRU暗号
まず、NTRU暗号が何かをきちんと理解する必要があります。
本セクションでは、「NTRU: A Ring-Based Public Key Cryptosystem」を参考に執筆しています。
そのため、最新のアルゴリズムとは違いますが、許してください。
## 1.1 表記の定義
NTRUは3つの整数パラメータ$(N, p, q)$と整数係数を持つ次数が$N-1$の4つの多項式集合$\mathcal{L}_f, \mathcal{L}_g, \mathcal{L}_Φ, \mathcal{L}_m$に依存しています。$\mathcal{L}_f, \mathcal{L}_g$は後述するRの係数が小さい部分集合とする。

ここで、$p,q$が素数である必要はないですが、$gcd(p, q) = 1$、つまり互いに素である必要があります。また、常に$p \ll q$である必要があります。

また、以降の説明では、環$R = \mathbb{Z}[X]/(X^N -1)$上であることにします。環の要素$F \in R$は多項式もしくはベクトルで表されます。
```math
F = \sum_ {i=0}^{N-1}F_ix^i = [F_0, F_1, \ldots, F_{N-1}]
```

$R$における乗算を$\circledast$で表すことにします。この乗算は循環畳み込み積として定義されています。
```math
F \circledast G = H\quad with\quad H_k = \sum_ {i=0}^{k}F_iG_{k-i} + \sum_ {i=k+1}^{k}F_i G_{N+k-i}=\sum_ {i+j\equiv k\mod N}F_iG_j
```
基本的に$F \circledast G$は$N^2$回の乗算が必要になりますが、NTRUで使用される典型的な積では、$F$もしくは$G$の一方が小さな係数を持つため、積の計算は非常に高速です。もし、$N$が大きい値の場合は高速フーリエ変換を使うと$O(N\log N)$となり、高速化できます。
## 1.2 鍵生成
Bobがランダムに多項式$f, g$を$\mathcal{L}_g$から選びます。
ここで、$f$は以下を満たす必要があります。
```math
F_q \circledast f \equiv1\mod q \quad and \quad F_p \circledast f \equiv1\mod p
```
準備ができたら、Bobは次の値を計算します。
```math
h \equiv F_q \circledast g \mod q
```
この多項式$h$が公開鍵となり、多項式$f$が秘密鍵となります。また、$F_p$もどこかに保存して必要があることに注意してください。

## 1.3 暗号化
Aliceは平文$m$を$\mathcal{L}_m$から選びます。
次にランダムな多項式$Φ$を$\mathcal{L}_Φ$から選びます。
最後にBobの公開鍵$h$を使って暗号文$e$を計算し、Bobに送信します。
```math
e \equiv pΦ \circledast h + m \mod q
```

## 1.4 復号
1.2で用意した秘密鍵$f$,多項式$F_p$を用いて復号します
まず、秘密鍵$f$を用いて暗号文$e$を掛け、多項式$a$を計算します。
```math
a = f \circledast e \mod q
```
ここで、多項式$a$の係数の範囲を$-\frac{q}{2}$から$\frac{q}{2}$の範囲となるように調整する。
次に、aを整数係数の多項式として扱い、平文$m$を復号する。
```math
m = F_p \circledast a \mod q
```
なお、パラメータ$N, p, q$が適切に設定されている場合は非常に高い確率で元の平文に復元されますが、一部パラメータ設定では計算の誤差などによって復号に失敗する可能性があります。そのため、各メッセージブロックにいくつかのチェックビットを含めるのが良いとされています。また、失敗の原因として、多項式$a$の係数が$[-\frac{q}{2}, \frac{q}{2}]$から外れていることが考えられます。そのため、範囲を少し広げ、$[-\frac{q}{2} +x, \frac{q}{2} +x]$とすると成功するかもしれません。

# 2. Broadcasting Attack
ここからは、Broadcasting Attackについてです。
以下は、CTF問題の説明文に記載されているものとなります。
そのため、問題を解くためには、以下の論文を理解する必要があります。

https://eprint.iacr.org/2011/590

## 2.1 多項式の行列表現
環の要素$F \in R$は以下のように表すことができました。
```math
F = \sum_ {i=0}^{N-1}F_ix^i = [F_0, F_1, \ldots, F_{N-1}]
```
また、$R$における乗算を$\circledast$で表し、この乗算は循環畳み込み積として定義されています。
そのため、$F,G\in R$のとき、$F \circledast G$は以下のように表すことができます。
```math
F * G = \begin{bmatrix}
F_0 & F_{N-1} & \cdots & F_1 \\
F_1 & F_0 & \cdots & F_2 \\
\vdots & \vdots & \ddots & \vdots \\
F_{N-1} & F_{N-2} & \cdots & F_0
\end{bmatrix}
\begin{bmatrix}
G_0 \\
G_1 \\
\vdots \\
G_{N-1}
\end{bmatrix}
```
この形式を使って、NTRUの暗号文$e$は
```math
e \equiv pΦ \circledast h + m \mod q
```
は次のような行列方程式に変換されます。ただし、論文の表記合わせのため、$r=pΦ$としています。
```math
e = \textbf{Hr} + \textbf{m} \mod q
```

ここで、行列 $\textbf{H}$ は公開鍵 $h$に基づく循環行列となります。

```math
\textbf{H} = \begin{bmatrix}
h_0 & h_{N-1} & \cdots & h_1 \\
h_1 & h_0 & \cdots & h_2 \\
\vdots & \vdots & \ddots & \vdots \\
h_{N-1} & h_{N-2} & \cdots & h_0
\end{bmatrix}
```
今回の(論文内で説明している)Broadcasting Attackでは$\mathbb{Z}_q^{N\times N}$内の循環行列のみを考えます。
## 2.2 Hrからrを取り出す
今回の攻撃はNTRU-1998、NTRU-2001、NTRU-2005の3つに対して有効だと述べられており、特にNTRU-1998は任意の平文$m$に対して攻撃が成立します。
以降はNTRU-1998と仮定して説明を行います(1で説明した理論が1998年なので、多分関係あると思います)。
NTRU-2001とNTRU-2005の攻撃はいつか実装してみたいですね。
さて、暗号文は行列を用いて以下で表すことができました。
```math
e = \textbf{Hr} + \textbf{m} \mod q
```
を用いて、$\textbf{H}$の可逆行列$\hat{\textbf{H}}$を定義します。
```math
b = \hat{\textbf{H}}e \mod q
```
このとき、以下が成立します。
```math
b = \hat{\textbf{H}}\textbf{m} + \textbf{r} \mod q
```
よって、以下の式を得られます。
```math
\textbf{r}  = b - \hat{\textbf{H}}\textbf{m} \mod q
```
ここで内積を考えます。
```math
  (b - \hat{\textbf{H}}\textbf{m})^T(b - \hat{\textbf{H}}\textbf{m}) = \textbf{r}^T \textbf{r} \mod q
```
$\textbf{r}^T \textbf{r} = d$とすると、
```math
  \textbf{m}^T\hat{\textbf{H}^T}\textbf{H}^T - 2\textbf{b}^T\hat{\textbf{H}}\textbf{m}=d- \textbf{b}^T \textbf{b} \mod q
```
ここで、$\hat{\textbf{H}^T}\textbf{H}^T$は対称行列となります(補題は省略。論文内のLemma(2.2)参照)。
## 2.3 線形化
簡略化のために、以下の置換を行います。
```math
d- \textbf{b}^T \textbf{b} = s, \quad \textbf{b}^T \textbf{H}^T  = (w_0, w_1, \cdots, w_{N-1})
```
また、行列積は次のようになります。

```math
\hat{\textbf{H}^T}\textbf{H}^T = \begin{bmatrix}
a_0 & a_{N-1} & \ldots & a_1 \\
a_1 & a_0 & \ldots & a_2 \\
\vdots & \vdots & \ddots & \vdots \\
a_{N-1} & a_{N-2} & \ldots & a_0
\end{bmatrix}
```

さて、2.2の最後の式を次のように書き換えます。

```math
a_0 \sum_{i=0}^{N-1} m_i^2 + a_1(m_1m_0 + m_2m_1 +\cdots +m_0m_{N-1}) +  \cdots+ a_{N-1}(m{N-1}m_0 + m_0m_1 + \cdots + m{N-2}m{N-1}) - 2 \sum_{i=0}^{N-1} w_i m_i = s \mod q
```

ここで、以下の置換を行います。
```math
x_i = \sum m_i m_{i+j}
```
これにより、次の線形方程式が得られます。

```math
a_0 x_0 + 2a_1 x_1 + \ldots + 2a_{\lfloor N/2 \rfloor} x_{\lfloor N/2 \rfloor} - 2 \sum_{i=0}^{N-1} w_i m_i = s \mod q
```

さらに、NTRUは特定の条件下で意味的安全性を失ってしまいます。
暗号文多項式を評価すると次のようになります:

```math
h(1)r(1) + m(1) = c(1) \mod q
``` 
NTRU-1998のとき$r(1) = 0$なので、暗号文から平文が漏れてしまいます。

よって、NTRU-1998のとき、以下の方程式を得られます。
```math
m_0 + m_1 + \cdots + m_{N-1} = c(1) \mod q  
``` 
この方程式の平方を考えましょう。すると、以下のようにまります。
```math
(m_0 + m_1 + \cdots + m_{N-1})^2 = (c(1))^2 \mod q  
``` 
ここで、以下の変数で置換することにしましょう。
```math
x_i = \sum m_j m_{j+i} \quad \text{for} \quad i = 0, 1, \ldots, \lfloor N/2 \rfloor
```
すると、以下のように表せます。
```math
x_0 + 2x_1 + \ldots + 2x_{\lfloor N/2 \rfloor} = c(1)^2 \mod q \quad
```
さて、2.3で得られた線型方程式と組み合わせることで以下のように書くことができます。
```math
2(a_1 - a_0)x_1 + \ldots + 2(a_{\lfloor N/2 \rfloor} - a_0)x_{\lfloor N/2 \rfloor} - 2 \sum w_i m_i = s - a_0 c(1)^2 \mod 2^k \quad 
```
2で割ると、
```math
(a_1 - a_0)x_1 + \ldots + (a_{\lfloor N/2 \rfloor} - a_0)x_{\lfloor N/2 \rfloor} - \sum w_i m_i = \frac{s - a_0 c(1)^2}{2} \mod 2^{k-1}
```

以上を行うことで、次の線形合同方程式が得られます。

```math
L \times Y = S \mod 2^{k-1}
```

ここで
- $ L $: $(N + \lfloor N/2 \rfloor) \times (N + \lfloor N/2 \rfloor)$ 行列
- $Y$ : $(x_1, x_2, \ldots, x_{\lfloor N/2 \rfloor}, m_0, m_1, \ldots, m_{N-1})^T$
- $ S $: 定数ベクトル

論文では、($L$|$S$)としてガウスの消去法で解くと書かれていました。

# 3 Broadcasting NTRU
## 3.1 問題概要
chall.sageとoutput.txtが渡されます。
```sagemath:chall.sage
import os
from Crypto.Cipher import AES
from hashlib import sha256

FLAG = os.environ.get("FLAG", "fakeflag")

N = 509
q = 2048
p = 3
d = 253

Zx.<x> = ZZ[]

def invertmodprime(f, p):
    T = Zx.change_ring(Integers(p)).quotient(x ^ N - 1)
    return Zx(lift(1 / T(f)))

def invertmodpowerof2(f, q):
    assert q.is_power_of(2)
    h = invertmodprime(f, 2)
    while True:
        r = balancedmod(convolution(h, f), q)
        if r == 1:
            return h
        h = balancedmod(convolution(h, 2 - r), q)

def balancedmod(f, q):
    g = list(((f[i] + q // 2) % q) - q // 2 for i in range(N))
    return Zx(g)

def convolution(f, g):
    return (f * g) % (x ^ N - 1)

def generate_polynomial(d):
    coeffs = [1] * d + [0] * (N - d)
    shuffle(coeffs)
    return Zx(coeffs)

def generate_keys():
    while True:
        try:
            f = generate_polynomial(d)
            g = generate_polynomial(d)
            f_p = invertmodprime(f, p)
            f_q = invertmodpowerof2(f, q)
            break
        except:
            pass
    public_key = balancedmod(p * convolution(f_q, g), q)
    secret_key = f, f_p
    return public_key, secret_key

def generate_message():
    result = list(randrange(2) for j in range(N))
    return Zx(result)

def encrypt(message, public_key):
    r = Zx(list(randrange(2) for j in range(N)))
    return balancedmod(convolution(public_key, r) + message, q)


msg = generate_message()

public_keys = []
ciphertexts = []
for _ in range(777):
    public_key, secret_key = generate_keys()
    ct = encrypt(msg, public_key)
    public_keys.append(public_key)
    ciphertexts.append(ct)

print("public keys:", public_keys)
print("ciphertexts:", ciphertexts)

key = sha256(str(msg).encode()).digest()[:16]
cipher = AES.new(key=key, mode=AES.MODE_CTR)
enc_flag = cipher.encrypt(FLAG.encode())
print("encrypted flag:", (cipher.nonce + enc_flag).hex())
```
重要なところだけをピックアップして説明します。
```sage
for _ in range(777):
    public_key, secret_key = generate_keys()
    ct = encrypt(msg, public_key)
    public_keys.append(public_key)
    ciphertexts.append(ct)
```
ここで、暗号文が777個生成されています。
さて、行列$Y$に存在する未知変数の数は$x$で約$N/2$個あり、$m$は$N-1$個あることを思い出しましょう。
- $Y$ : $(x_1, x_2, \ldots, x_{\lfloor N/2 \rfloor}, m_0, m_1, \ldots, m_{N-1})^T$

約$3N/2$個の未知変数があることになります。
つまり、今回約3N/2個の暗号文が存在すれば$m$が復元できそうです。パラメータを見ると、N = 509なので、$3\times 509/2 = 763.5 < 777$となり、無事足りそうです。

そして、以下の処理についてです。
```sage
def encrypt(message, public_key):
    r = Zx(list(randrange(2) for j in range(N)))
    return balancedmod(convolution(public_key, r) + message, q)
```
2.2 にて、$\textbf{r}^T \textbf{r} = d$としています。しかし、上記の処理を見ると、ランダムなバイナリ列となっています。
つまり、係数1が$d_r$個、係数-1が$d_r$個となり、$\textbf{r}^T \textbf{r} = 1^2 d_r + (-1)^2 d_r = 2d_r$となるように修正すべきです。
これは、2倍して1引くことで解決できます。1のときは、1のままになり、0の時は-1に変更することができます。
以上の情報を使って実装を行っていきましょう。
## 3.2 攻撃の実装
基本的には、上のコードとチョコラスクさんのコードを参考にして作成します。
```sage
from Crypto.Cipher import AES
from hashlib import sha256

N = 509
q = 2048
p = 3
d = 253

Zx.<x> = ZZ[]


def invertmodprime(f, p):
    T = Zx.change_ring(Integers(p)).quotient(x ^ N - 1)
    return Zx(lift(1 / T(f)))

def invertmodpowerof2(f, q):
    assert q.is_power_of(2)
    h = invertmodprime(f, 2)
    while True:
        r = balancedmod(convolution(h, f), q)
        if r == 1:
            return h
        h = balancedmod(convolution(h, 2 - r), q)

def balancedmod(f, q):
    g = list(((f[i] + q // 2) % q) - q // 2 for i in range(N))
    return Zx(g)

def convolution(f, g):
    return (f * g) % (x ^ N - 1)

def calc(pubkeys, cipher):
    pol_one = sum(x^i for i in range(N))
    mat = []
    vec = []

    for h, c in zip(pubkeys, cipher):
        try:
            h1 = invertmodpowerof2(h, q)
        except:
            continue

        h1 *= 2
        b = balancedmod(convolution(h1,c)-pol_one, q)
        s = (N-sum(v*v for v in list(b))) % q
        vec.append(s)

        ht = list(h1)
        ht += [0]*(N-len(ht))
        ht = ht[::-1]
        ht = Zx([ht[-1]]+ht[:-1])
        a = balancedmod(convolution(ht,h1), q)
        w = balancedmod(convolution(ht,b), q)

        a = list(a)
        w = list(w)

        a += [0]*(N-len(a))
        w += [0]*(N-len(w))

        mat_row = [a[0]]+[2*a[i] for i in range(1,N//2+1)]+[-2*w[i] for i in range(N)]
        mat.append(mat_row)

    return mat, vec

with open('output.txt') as f:
    pubkeys = sage_eval(f.readline()[len("public keys: "):].strip(), locals={'x':x})
    cipher = sage_eval(f.readline()[len("ciphertexts: "):].strip(), locals={'x':x})
    enc_flag = bytes.fromhex(f.readline()[len("encrypted flag: "):].strip())

mat, vec = calc(pubkeys, cipher)
mat = matrix(Zmod(q), mat)
vec = vector(Zmod(q), vec)
res = mat.solve_right(vec)
res = [int(v) for v in list(res)[-N:]]

msg = balancedmod(res, 4)
key = sha256(str(msg).encode()).digest()[:16]
cipher = AES.new(key=key, mode=AES.MODE_CTR, nonce=enc_flag[:8])
flag = cipher.decrypt(enc_flag[8:])
print(flag)
```
説明していきます。
```
def calc(pubkeys, cipher):
    pol_one = sum(x^i for i in range(N))
    mat = []
    vec = []
```
pol_oneは係数が1の多項式を用意しています。matは線形合同方程式の行列を、vecはベクトルを格納します。
```
for h, c in zip(pubkeys, cipher):
    try:
        h1 = invertmodpowerof2(h, q)
    except:
        continue
```
公開鍵$h$の逆元を計算します。
逆元が存在しない場合は次に進みます。
```
h1 *= 2
b = balancedmod(convolution(h1, c) - pol_one, q)
s = (N - sum(v*v for v in list(b))) % q
vec.append(s)
```
まず、以下の式
```math
b = \hat{\textbf{H}}e \mod q
```
を計算しています。その後、内積を計算し、vec配列に格納しています。ここで、3.1で説明した2倍することを忘れないようにしてください。

```
        ht = list(h1)
        ht += [0]*(N-len(ht))
        ht = ht[::-1]
        ht = Zx([ht[-1]]+ht[:-1])
        a = balancedmod(convolution(ht,h1), q)
        w = balancedmod(convolution(ht,b), q)

        a = list(a)
        w = list(w)

        a += [0]*(N-len(a))
        w += [0]*(N-len(w))

        mat_row = [a[0]] + [2 * a[i] for i in range(1, N//2+1)] + [-2 * w[i] for i in range(N)]
        mat.append(mat_row)
```
ここらの計算うまく動かなかったので、少しチョコラスクさんのwriteupからお借りしました。
a,wの行列要素を計算後、行を構築、mat配列に格納をしています。

後は、ガウスの消去法で解けば結構早い段階でフラグを獲得することができます。

# 終わりに
NTRU暗号僕が学んだ時より理解できたのでかなり楽しかったです。学びの歩みを止めなくてよかったと思います。
後、新しいことを学ぶたびに頭が締め付けられるような痛みに襲われたのがとてもよかったです。学んでるって感じがするので。
もし、今回の説明で間違いがありましたら**優しく**教えていただけると助かります。
最後まで読んでいただきありがとうございました。
良いお年をお迎えください。僕は修論完成させてきます。
# 補足
##  A. 環とは？
環は集合$R$とその上の2つの二項演算$(+, \cdot)$の組$(R, +, \cdot)$であり、以下が成立することをいいます。ただし、$a,b \in R, a+b, ab\in R$です。

1. $(R, +)$が可換群である。
    a. **結合法則**: $(a + b) + c = a + (b + c)$
    b. **可換性**: $a + b = b + a$
    c. **単位元の存在**: $ a + 0 = 0 + a = a$
    d. **逆元の存在**: $a + (-a) = (-a) + a = 0$

2. $(R, \cdot)$がモノイドである。
    a. **結合法則**: $(ab)c = a + (bc)$
    b. **単位元の存在**: $ a1 = 1a = a$

3. $+, \cdot$間の分配法則が成立する。
    a. $a(b + c) = ab + ac$
    b. $(a + b)c = ac + ab$