Schnorr 署名に使われる数学 - 初等整数論

前置き

本記事は Nostr Advent Calendar 第一会場 12月3日 分の記事として書きました。

Nostr Advent Calener 第一会場

序論

Schnorr 署名にかかわる演算のほとんどの部分は整数座標上で定義された楕円曲線上の演算で実行されます。

その楕円曲線上の演算を更に要素に分けていくと、結局のところ、整数の演算の組み合わせとなります。

幸いなことに、Schnorr 署名で使われる範囲の演算ならば、そこさえ見通せば署名のアルゴリズムを概観するのには十分…というのが私の所感です。

具体的には以下を把握できれば、楕円曲線上の計算を具体的なプログラムとして実装できると考えます。

• 整数の表現を把握する
• 四則演算を定義
• 拡張ユークリッド互除法による割り算（乗法逆元）の計算
• 巨大な指数による累乗計算（バイナリ法）
• 平方剰余の計算

以下、一つ一つ見ていきます。

整数表現

素数を法とする合同式による表現

楕円曲線に関わる計算は、ほぼ全てが、特定の素数を法として合同な整数上で行われます。

上記の定義は、ちょっとわかりにくいかもしれないので、例として、素数 $p=13$ を法として合同な整数を例として挙げます。

$13$ を法として合同な整数は、13で割った数の余りが等しいものを同じ数としてみなします。つまり、

• $0 = 13 \pmod {13} $
• $1 = 14 \pmod {13} $
• $2 = 15 \pmod {13} $
• $3 = 16 \pmod {13} $
• …
• $11 = 24 \pmod {13} $
• $12 = 25 \pmod {13} $
• $13 = 26 = 0 \pmod {13}$

となります。

ここで、$\pmod {13}$ は 左式が13を法として合同であることを示す表記方法です。

ところで、この記事を見るのはプログラマーの方やソフトウェア開発に興味のある方が多いと思われます。

そのようなバックグラウンドのある方にとって、$p$ を法として合同の数とは「足し算や掛け算でp以上になるとオーバーフローしてそれ以上の数が無視される」といった法がわかりやすいかもしれません。

つまり、

$7 + 8 = 15 = 13 + 2 = 2 \pmod {13}$

と、計算結果が13以上にあふれた場合は、ゼロに戻るということです。

引き算に関しても同じで、0を下回るとアンダーフローし最大値の12に飛んでいくイメージです。

$7 - 8 = -1 = 13 -1 = 12 \pmod p$

必然的に、$p$ を法として合同とする数は、その範囲が $0 .. p-1$ となります。

整数のサイズ

Schnorr 署名では、ほとんどの計算にて 64バイトの整数が用いられます。

使われている法としての素数のサイズも64バイトです。

例えば、ビットコインやNostrで署名として使われる scpe256k1 という体系では、

$p$ = 0xFFFF FFFF FFFF FFFF FFFF FFFF FFFF FFFF FFFF FFFF FFFF FFFF FFFF FFFE FFFF FFC2F

という巨大な素数が使われています。

この巨大な素数を法として用いることで、公開鍵から秘密鍵を探索することを難しくしています。

足し算/引き算/掛け算

先ほど述べた通り、剰余が同じ整数は、同じものとみなします。

よって、整数による足し算、引き算、掛け算を行った後、法とする整数で割った余りをもちいることで、合同式上の加減算・乗算を実現できます。

$z = x + y \pmod p$

は、計算機上では、

$z = r$, ($r$ は $x + y$ を$p$で割った余り)

とすればOKです。

この時注意すべきは、$x + y$, $x - y$, $x \times y$ という演算はどのパターンも 64バイトの整数が扱える範囲をオーバーする可能性があることです。

それを考慮に入れた計算を実装するか、多倍長整数の計算ライブラリを使う必用があります。

割り算と拡張ユークリッド互除法による乗法逆元の計算

楕円曲線上で用いる合同式の四則演算の中で、割り算だけが特殊です。

通常、プログラミング言語で標準的に実装されいている整数の割り算は、計算結果の小数点以下を切り捨てることで実現します。

ですが、ここでの割り算 $z = x \div y \pmod p$ は、$ y \times w = 1 \pmod p$ となる $w$ を求めて、$z = x \times w \pmod p$ を求める計算を行う必要があります。

つまり、$y$の 逆数に相当する整数をめます。その数は乗法逆元と呼ばれています。

乗法逆元の求め方 - 拡張ユークリッド互除法による解放

ここでは、拡張ユークリッド互除法による乗法逆元の求め方を紹介します。

拡張ユークリッド互除法は、$a$,$b$を整数、$g$を$a$,$b$の最大公約数として、以下の方程式、

$ax + by = g$

を満たす$x$ と $y$ を求める解法です。

ここで $p$ を素数として、$a=p$ とします。$b$ を $b<p$ なる数とし、、$b$の乗法逆元を求めたいとします。

この場合、$p$ と $b$の最大公約数は明らかに $1$ なので、

$px + by = 1$

です。

ここで、両辺の素数pを法とした合同式を考えると、$py \pmod p = 0$ であるため、

$by = 1 \pmod p$

が得られます。よって、$y$ が拡張ユークリッド互除法で求められていたとしたら、定義より、$y$ は $b$ の乗法逆元ということになります。

拡張ユークリッド互除法のアルゴリズム

拡張ユークリッド互除法は以下の手順によって、方程式をより小さい数によるものに帰着していくものです。

巨大な数を使ったとしても、場合でもそれほど多くない繰り返しで解にたどり着くことが知られています。

また冠している名前「ユークリッド」から推察できるように、古代から知られ、隅から隅まで研究しつくされている解法でもあります。 具体的な手順以下の通りです。

手順を見ることで「方程式をより小さい数を使ったものに帰着している」の感じがつかめるのではないかと思います。

• $m_1 = p$, $n_1 = b$ とする
• $m_1x_1 + n_1y_1 = 1$
• $m_1$ を $n_1$ で割った商を $q_1$, 余りを$r_1$とする。
• 書き換えて、$(n_1q_1+r_1)x_1 + n_1y_1 = n_1(q_1x_1+y_1) + r_1 x_1=1$
• 新たな変数 $x_2$, $y_2$を導入し、逐次的により小さい問題に変換することを考える。
• $x_2 = q_1 x_1 + y_1$, $y_2=x_1$, $m_2=n_1$, $n_2=r_1$, とし、新たな方程式、$m_2x_2 + n_2y_2=1$ を得る。
• このより小さい係数からなる新しい方程式が解けた場合、与式より、$x_1 = y_2$, $y_1=x_2-q_1 y_2$ として求められる。
• これを、最終的に $r_n = 1$ となるまで繰り返せばよい。このとき $m_n x_n + 0 y_n = 1$ なので、$m_n=1$, $x_n = 1$ が最終的に到達する点であり、$x_{n+1}=1$,$y_{n+1}=0$ と定めれば、一連の計算を逆にたどることで解にたどり着けることが分かる。

一連のステップは行列として書くことができて、繰り返し計算に還元するにはこちらを使うのが簡単です。

つまり、

$x_1 = y_2$, $y_1=x_2-q_1 y_2$ より、

$$ \begin{pmatrix} x_1 \\ y_1 \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ 1 & -q_1 \end{pmatrix} \begin{pmatrix} x_2 \\ y_2 \end{pmatrix} $$

と書けることから、$x_n$,$y_n$ まで来た段階で計算をストップできると過程すると、

$$\begin{pmatrix} x_1 \\ y_1 \end{pmatrix} = \begin{pmatrix} 0 & 1 \\ 1 & -q_1 \end{pmatrix} \begin{pmatrix} 0 & 1 \\ 1 & -q_2 \end{pmatrix} \begin{pmatrix} 0 & 1 \\ 1 & -q_3 \end{pmatrix} \cdots \begin{pmatrix} 0 & 1 \\ 1 & -q_{n-1} \end{pmatrix} \begin{pmatrix} 0 & 1 \\ 1 & -1\end{pmatrix} \begin{pmatrix} 1 \\ 0 \end{pmatrix} $$

です。行列の積の性質 $(AB)C = A(BC)$ から反復的に解いていくには行列の先頭から順に計算していけばよいことが分かります。

このようにして乗法逆元を求め、それをもとに割り算を実行しています。

巨大な数を法とする合同式上の割り算は、拡張ユークリッド互除法等で高速に計算できるのが実用に供することのできる一つのポイントになっています。

巨大な指数による累乗計算（バイナリ法）

累乗は要素としては掛け算の範疇ではありますが、大事なポイントのひとつとして、巨大な指数による整数の累乗計算の高速計算方法を紹介します。

巨大な指数による累乗計算は、合同式上の平方剰余の計算を高速に行う際に必要になります。 この方法のエッセンスは、途中の計算過程をうまく利用することにより、本来なら$O(n)$の計算を$O(\log(n))$で実行することです。

この方法は、高速化の良く見るパターンのため、もしこれに馴染みがあまり無い場合には、学習しておく価値があると信じます。

以下概要です。

今ある整数 $x$ の正の整数 $m$ による累乗を計算するとします。

$S = x^{m}$

ナイーブな計算方法は、$x$ を $m$ 回掛けることです。$m$ が例えば、前述の scep256k1 の $p$ のような場合、計算に途方もない時間が掛かることがわかると思います。

高速化のためまず、指数の表現を書き換えます。

$m$をN桁の二進数表現、つまり、

$m = \sum_{i=0}^{N} d_i 2^i$ $(d_i = 0 , or , 1)$

と書くこととします。こちらを使ってあらためて累乗の式を書き直すと、

$S = \prod_{i=0}^N\left( x^{2^i d_i} \right)$

となり、N回の剰積の中に $x^{2^i d_i}$ が含まれている形式となりました。$N$ は $m$ を二進数表現した時の桁数になります。

例えば、64バイトの整数なら $N=8 \times 64=512$ です。

$d_i$ は 二進数表現に応じて$0$か$1$です。書き換えると、

$x^{2^i d_i} = x^{2i}$ if $d_i=1$ else $d_i = 1$

です。

ここで、工夫するの $x^{2i}$ の計算です。

i から N までの反復中、すでに、 $x^{2^i}$が計算されていれば、 $x^{2(i+1)}$ は、

$x^{2(i+1)} = x^{2}x^{2i}$

として計算できます。つまり、

 // 疑似コード

 x // 累乗したい数
 d(0 to N) // 指数の二進数表現

 S = 1 // 計算結果
 X = 1 // xの累乗途中結果保存

 for int i = 0 to Length(d)
  if d(i) = 1 then
    S = S*X
  end if
  X = X * x^2
 Next i

と、途中の乗積結果を再利用する形で計算すれば、たかだか $x$ の桁数のオーダー、つまり $O(log(x))$ で計算を実行することができます。

平方剰余の計算

平方剰余とはいわゆる合同式上の平方根のことです。 実数の場合とはことなり、どんな場合でも平方根に相当するものが存在するわけではありません。

ある数$p$を法として合同のな整数の平方剰余は以下の定義となります。

$y^2 = x \pmod p$

このとき、y を x の平方剰余と呼びます。例として、$p=7$の場合の各元の自乗から、平方剰余をリストアップしてみます。全ての数に平方剰余が存在するわけではないことが概観できるはずです。

このように、1,2,4 にはそれぞれ 1 or 6, 3 or 4, 2 or 5 が平方剰余として定義できる半面、3,5, には平方剰余、つまり二乗して 3, 5 になる数 は存在しないことになります。

実数の計算と異なり、平方剰余では、実数でのニュートン法のような収束による高速計算が存在しないことが知られています。

今回は法とする整数$p$を4で割った際に余りが3となる場合に適用できる、フェルマーの小定理に基づいた平方剰余の計算方法を紹介します。 これは、平方剰余の計算を累乗計算に帰着させる方法で、前述のバイナリ法によって高速計算が可能です。

また、secp256k1 に用いられている素数はまさに、$4$ で割ったとき $3$ あまるケースに該当しています（と言うよりそうなるように選んでいる）。

フェルマーの小定理

フェルマーの小定理は、素数 $p$ を法として合同な整数 $x$ が以下の性質を持つという定理です。

$ x^{p-1} = 1 \pmod p$

平方剰余の計算

$p$ は素数であり奇数と考えられる $(p-1)/2$ は整数です。よって $x^{(p-1)/2}$ も合同式上の整数ですが、$p$ を法とする場合、二乗して $1$ になる整数は $p-1$ もしくは $1$ です。

$x^{(p-1)/2} = 1 , or , p-1 \pmod p$

両辺に $x$ を掛けます。

$x^{(p+1)/2} = x , or , xp-x = x , or , p - x \pmod p$

ここで $p$は $4$ で割って $3$ 余ることから、$(p+1)/4$ も整数なので、

$(x^{(p+1)/4})^2 = x , or , p-x \pmod p$

となり、$x$ および、$p-x$ の平方剰余が $x^{(p+1)/4}$ であることが示されました。

累乗計算に、先述のバイナリ法にを用いることで、平方剰余が比較的低い計算コストで求めることができます。

結論

以上です。

ここまで、Schnorr 署名に必用な楕円曲線上の計算に必要な要素に関して、初等整数論に属している（と私が考える部分）に関して述べました。

計算の枠組みとしての合同式は、日頃から整数のデータ形式に馴染んでいるプログラマーの方々には、特に馴染みやすい概念なのではないかと思います。

拡張ユークリッド互除法、バイナリ法は、それぞれ逐次計算、計算過程で得られる途中結果を上手く利用することで、ナイーブなアルゴリズムから劇的に計算方法を抑えることができる例でした。

フーリエ変換や畳み込み計算の高速計算方法へとつながる知識でもあります。このようなシンプルなアルゴリズムに立ち戻ることによって、これらのもう少し複雑な例の理解の役に立つと考えます。

上記例は、歴史も長く、標準の整数計算ライブラリには組み込まれていることもあると思います。

ですが、もし、多倍長計算で速度に問題を感じた場合、上記のような高速計算が適用されているか、確認してみるとよいのではないかと思います。

次回で、署名に必要な楕円曲線上の計算のに進みたいと思います。

Reference

• Nostr Implementation Posibilities #01 - NIP-01
• Bitcoin Improvement Proposals #0340 - BIP0340