🧮

暗号で使われる標数2の体の実装

2023/07/06に公開

aes

拡大体

tech

初めに

今回は標数2の体の基本的な実装方法をいくつか紹介します。標数2の体自体の話は暗号で使われる標数2の体を参照ください。
サンプルコードはgf256にあります。

多項式のビット表現

AES本体で使われる8次拡大体 $K=𝔽_{2^8}$ の実装方法を考えます。
$𝔽_2$ の要素は1ビットで表せるので8次拡大体 $K$ の要素（ $𝔽_2$ 係数7次多項式）は8ビットで表現できます。8ビット整数 $a$ の $i$ ビット目 $a_i$ が多項式の $i$ 次の係数としましょう（ $i=0, \dots, 7$ ）。
たとえば $a$ = 0b01101101なら対応する多項式は $X^6 + X^5 + X^3 + X^2 + 1$ です。

$a$ = 0b01101101に対応する多項式

ビット位置	7	6	5	4	3	2	1	0
$a_i$	0	1	1	0	1	1	0	1
多項式	0	$X^6$	$X^5$	0	$X^3$	$X^2$	0	1

多項式の加算

2個の多項式 $a$ と $b$ の加算を考えます。多項式の加算は同じ次数の係数同士を足したものでした。

a+b=(a_7 X^7 + \cdots + a_0)+(b_7 X^7 + \cdots + b_0) = (a_7 + b_7)X^7 + \cdots + (a_0+b_0).

各係数は

𝔽_2

の要素なので加算は排他的論理和です。それをビットごとにすればよいので、結局整数として排他的論理和を取ればよいです。引き算もaddと同じですね。

typedef uint8_t K;
K add(K a, K b)
{
  return a ^ b;
}
K sub(K a, K b)
{
  return add(a, b);
}

多項式の乗算

まずはナイーブな実装から考えましょう。
多項式 $a$ と $b$ を掛けた結果 $c=ab$ を $f=X^8+X^4+X^3+X+1$ で割ります。 $f$ が $n=8$ 次なので $c$ は（最大） $2n-2=14$ 次多項式です。

c=c_{14} X^{14} + c_{13} X^{13} + \cdots + c_1 X + c_0=ab=(a_7 X^7 + \cdots + a_0)(b_7 X^7 + \cdots + b_0).

\begin{cases} c_0 &=& a_0 b_0,\\ c_1 &=& a_0 b_1 + a_1 b_0,\\ c_2 &=& a_0 b_2 + a_1 b_1 + a_2 b_0,\\ ... &=& ...,\\ c_{13} &=& a_6 b_7 + a_7 b_6,\\ c_{14} &=& a_7 b_7. \end{cases}

$c$ の次数 $i$ について前半 $i=0, \dots, n-1=7$ と後半 $i=n=8, \dots, 2n-2=14$ に分けて $\sum$ 記号を使うと

c_i = \begin{cases} \sum_{j=0}^i a_j b_{i-j} & (i=0, \dots, n-1),\\ \sum_{j=i-n+1}^{n-1} a_j b_{i-j} & (i=n, \dots, 2n-2). \end{cases}

ここの足し算や掛け算も $𝔽_2$ の要素として計算します。記号で書くと複雑に見えますが、単なる積和のループです。C++の例 mulPoly を示します。
まず体 $K$ とその2倍のサイズの16ビット整数型 $K2$ を定義し、 $i$ ビット目の値を取得するgetと設定するsetも準備します。

typdef bool F2;
typedef uint8_t K;
typedef uint16_t K2;
template<class T>
F2 get(T a, int i) { return (a >> i) & 1; }
template<class T>
void set(T& a, int i, F2 b)
{
    a &= ~(1ull<<i);
    if (b) a |= 1ull<<i;
}

これらを用いてmulPolyを実装します。

K2 mulPoly(K a, K b)
{
    const int n = 8;
    K2 ret = 0;
    for (int i = 0; i < n; i++) {
        F2 r = 0;
        for (int j = 0; j <= i; j++) r ^= get(a, j) & get(b, i-j);
        set(ret, i, r);
    }
    for (int i = n; i < 2*n-1; i++) {
        F2 r = 0;
        for (int j = i-n+1; j < n; j++) r ^= get(a, j) & get(b, i-j);
        set(ret, i, r);
    }
    return ret;
}

前半のループが $i=0,\dots, n-1$ , 後半のループが $i=n, \dots, 2n-2$ の計算に相当します。
modPolyにより $2n-2=14$ 次の多項式ができました。これを $f=X^8+X^4+X^3+X+1$ で割る関数 modPoly を作ります。
割り算の筆算と同じく、多項式の割り算は上の桁（上の次数）から計算します。
たとえば $X^{13}$ があれば

X^{13}=X^5 X^8 = X^5(f + (X^4+X^3+X+1)) = X^5 f + X^5(X^4 + X^3 + X + 1) \equiv X^5(X^4+X^3+X+1).

つまり

X^8

を見つけたら

X^4+X^3+X+1

に置き換えます。ロジックとしては

X^{13}

を消して

X^5(X^4+X^3+X+1)

を加算します。

次数	13	12	11	10	9	8	7	6	5	4	3	2	1	0
元の多項式	1	.	.	.	.	.	.	.	.	.	.	.	.	.
mod後の多項式	0	.	.	.	+1	+1	.	+1	+1	.	.	.	.	.

この操作を次数 $2n-2=13, 12, \dots, 8$ までやれば $f$ で割った余りになります。 $X^4+X^3+X+1$ のビット表現は0b11011=0x1bなので次のコードになります。

K modPoly(K2 c)
{
    const int n = 8;
    for (int i = 2*n-2; i >= n; i--) {
        if (get(c, i)) {
            c ^= 0x1b << (i-n);
        }
    }
    return c;
}

速度改善

前節は多項式を掛けてからmodしました。そうではなく、掛け算しながらmodすることを考えます。
そのためにまず多項式 $a$ に $X$ を掛けてみましょう。 $X$ を掛けると次数が1増えます。掛ける前が7次未満なら $f$ で割る必要はありません。
掛ける前が7次なら $X$ を掛けて8次になるので先程と同じく0x1bを足します。

K mulX(K a)
{
    if (a & 0x80) {
        // 7次なら0x1bを足す
        a = (a << 1) ^ 0x1b;
    }else{
        // 7次未満なら単にXを掛けるだけ
        a <<= 1;
    }
    return a;
}

多項式 $a$ と $b$ の積は $b$ だけ展開すると

ab = a(b_0 + b_1 X + b_2 X^2 + \cdots) = ab_0 + (aX)b_1 + (aX^2)b_2 + \cdots.

a, aX, aX^2, \dots

は次々に

X

を掛ければよく、

b_0

b_1

は0か1なので乗算は次のように出来ます。

K mul(K a, K b)
{
    K ret = 0;
    while (b) {
        if (b & 1) ret ^= a; // 最下位ビットがあればret += a;
        a = mulX(a); // a ← aX
        b >>= 1; // 次のビット
    }
    return ret;
}

最初の多項式を掛けてからmodする場合の演算量はざっくり $O(n^2)$ ですが、後者の場合は $O(n)$ です。実際簡単なベンチマークを取ると、370clk→25clkと10倍以上速くなりました。

テーブル引き（その1）

今回入力が8ビット整数なので $a$ , $b$ をインデックスとするテーブル引きを使ってみましょう。
256×256×1バイト=65KiBと少々大きいですが、コードは簡単になります。

static uint8_t g_mulTbl[256 * 256];
// 一度だけテーブル生成
for (uint32_t x = 0; x < 256; x++) {
    for (uint32_t y = 0; y < 256; y++) {
        g_mulTbl[x + y * 256] = gf256_mul(x, y);
    }
}

uint8_t mulTbl(uint8_t x, uint8_t y) // テーブル引きによる乗算
{
    return g_mulTbl[x + y * 256];
}

ベンチマークをとるとコードの単純さの割りに10.3clkと前節の2倍少々しか速くなりませんでした。やはりメモリアクセスは負荷が多いようです。
64KiBはテーブルとして大きすぎる場面もあるでしょう。そこでテーブルサイズを512バイトと128分の1にする方法を紹介します。

テーブル引き（その2）

天下りですが、多項式 $e=1+X$ 、つまり $e=3$ のべき乗を考えます。

e^0=1, e^1=3, e^2=5, \dots

べき乗を次々と計算すると

e^{254}

=0xf6,

e^{255}

=1 と255乗したところで初めて1に戻ります。
0は何乗しても0で変わらないので除外すると、

K=𝔽_{2^8}

の0でない要素

a

はある

A

を使って

a=e^A

と表せることになります（このような

e

を

K^{*}

の生成元といいます）。
これを

A=\log(a)

と書くことにしましょう。

a=e^{\log(a)}

です。
すると

a

と

b

の乗算は

ab = e^{\log(a)} e^{\log(b)} = e^{\log(a)+\log(b)}

と変形できます。

0 \le i < 256

に対する

\{e^i\}

と

\log

のテーブルを用意すると

a, b

に対してlogのテーブル引きをし、加算してから指数のテーブル引きをすれば乗算できます。
実際には次のように書けます。0を除外するので指数の範囲は0以上255未満であることに注意してください。

uint8_t gf256_mul(uint8_t x, uint8_t y)
{
    if (x == 0 || y == 0) return 0;
    // x y = exp(logX) exp(logY) = exp(logX + logY)
    uint8_t logX = logTbl[x - 1];
    uint8_t logY = logTbl[y - 1];
    uint32_t z = logX + logY;
    // mod 255
    if (z >= 255) z -= 255;
    assert(z < 255);
    return expTbl[z];
}

指数と対数のテーブルはそれぞれ256個（厳密には255個）なので、合計512バイトです。ベンチマークを取ると全てテーブル引きした場合の1.6~1.8倍の遅さできた。
テーブルサイズは100分の1になりましたが3回テーブル引きしているので遅くなりました。

除算は、 $a$ に対する $1/a$ のテーブルを用意してテーブルルックアップして計算します。こちらも逆元テーブル256バイト必要です。
しかし、指数と対数のテーブルを用意しておけば、 $1/a=\exp(-\log(a))$ なので逆元テーブルは不要になります（ただしテーブル引きは2回）。このあたりはトレードオフですね。

まとめ

$𝔽_2$ の拡大体、特に $𝔽_{2^8}$ における四則演算の実装方法を紹介しました。ビット演算で処理する方法と、2種類のテーブル引きによる方法です。
アーキテクチャによってどれがよいかは変わるでしょう。次回はx64における専用命令を使った実装を紹介します。

GitHubで編集を提案