楕円曲線のスカラー倍算アルゴリズム Part 1

昨年9月にパテントロックが外れたことで、Bitcoin Coreでもスカラー倍算アルゴリズムの高速化オプションが有効化された↓

もともと、パテントの存在を知らずHal Finneyによって2013年にlibsecp256k1には実装されてたものの無効化されていたものが有効になったらしい↓

楕円曲線暗号においては、楕円曲線のスカラー倍算の計算をすることが多い。

秘密鍵から公開鍵を算出するのも、楕円曲線のベースポイントGに秘密鍵をスカラー値としたスカラー倍算であるし、ECDSAの署名検証アルゴリズムや、ECDHで共有鍵を導出するといったケースでも。

このスカラー倍算は、素直に実装すると乗算するスカラー値分だけポイントを加算することで計算できる。例えば、xPを計算する場合、Pをx回加算することでxPを求めることができる。ただ、素直にこれを計算すると遅く効率が悪いので、スカラー倍算を高速にするためのアルゴリズムがいろいろある。ということで、実際にどんなアルゴリズムがあって、現在Bitcoin Coreではどんなアルゴリズムを実装しているのか調べてみよう（ボリューム的にBitcoin Coreのアルゴリズムは後日別の記事で）。

バイナリ法

楕円曲線のスカラー倍算を多項式時間で実行するための有名なアルゴリズムがバイナリ法。乗算するスカラー値xをバイナリ形式にし、

ビットが1の場合は、加算と2倍算を行う。
ビットが0の場合は、2倍算のみを行う。

ことで、計算をより高速に行う。

例えばx = 5とした場合、それをバイナリ展開すると101（ ${1 * 2^{2} + 0 * 2^{1} + 1 * 2^{0} }$ ）になる。5Pをバイナリ法で計算すると以下のようになる。

計算結果をQとし初期化(無限遠点Oと)する。計算中に加算する点の一時変数v = Pと初期化する。
最下位ビットは1なので、Qにvを加算し、vを2倍算しv = 2v（この時点でv = 2P）とし、ビットシフトする。
続く最下位ビットは0なので、vを2倍算しv = 2v（つまりこの時点でv = 4P）とし、ビットシフトする。
続く最下位ビットは1なので、Qにvを加算し、vを2倍算しv = 2v（この時点でv = 8P）し、ビットシフトする。
処理するビットが無いので、Qが計算結果。

手順2でQ = Pがセットされ、手順4でさらに4Pが加算されQ = 5Pになる。これがバイナリ法のアルゴリズムで、計算処理で実際に実行されるのは2倍算と加算。

Rubyのecdsa gemでも楕円曲線のスカラー倍算ではこのアルゴリズムを使用している。実装もシンプル↓

# Returns the point multiplied by a non-negative integer.
#
# @param i (Integer)
# @return (Point)
def multiply_by_scalar(i)
  raise ArgumentError, 'Scalar is not an integer.' if !i.is_a?(Integer)
  raise ArgumentError, 'Scalar is negative.' if i < 0
  result = group.infinity
  v = self
  while i > 0
    result = result.add_to_point(v) if i.odd?
    v = v.double
    i >>= 1
  end
  result
end

↑のコードのv.doubleが2倍算、result.add_to_point(v)が加算処理だが、これらは有限体上で計算され、具体的には以下の計算が行われている。

有限体上の加算

アフィン座標の楕円曲線の方程式を ${y^{2} ≡ x^{3} + ax + b \mod p}$ とした場合、２つの座標点 ${P = (x_1, y_1), Q = (x_2, y_2) }$ の加算 ${P + Q = (x_1, y_1) + (x_2, y_2) = (x_3, y_3)}$ （P != Q）の計算は、以下の式で計算される。

$x_{3} = \left( \frac{y_{2} - y_{1}}{x_{2} - x_{1}} \right)^{2} - x_{1} - x_{2} \mod p$
$y_{3} = \left( \frac{y_{2} - y_{1}}{x_{2} - x_{1}} \right) (x_{1} - x_{3}) -y_{1} \mod p$

ecdsa gemのadd_to_pointのコードもこの計算をしてる。なお、 $\frac{y_{2} - y_{1}}{x_{2} - x_{1}}$ は点Pと点Qを結ぶ直線の傾き。

有限体上の2倍算

P = Q、つまり2Pの計算は、以下の式になる。

$x_{3} = \left( \frac{3x_{1}^{2} + a}{2y_{1}} \right)^{2} - 2x_{1} \mod p$
$y_{3} = \left( \frac{3x_{1}^{2} + a}{2y_{1}} \right) (x_{1} - x_{3}) -y_{1} \mod p$

ecdsa gemのdoubleのコードもこの計算をしてる。

ウィンドウ法

バイナリ法を少し変更したものがウィンドウ法で、ウィンドウ法ではまずウィンドウサイズwを決め、スカラー値xをm進展開する。ここでmはm = 2^wで、ウィンドウサイズw = 1の場合↑のバイナリ展開になる。

例えばw = 4とした場合のm進展開(16進展開)であれば、スカラー値xは、

${x = x_0 * 2^{4 * 0} + x_1 * 2^{4 * 1} + x_2 * 2^{4* 2} + ... + x_j * 2^{4 * j}}$

と展開され、各 ${x_i}$ は ${x_i \in｛0, 1, ..., m-1(15)｝}$ の値を取る。

例えばx = 515であれば、 ${x = 3 * 2^{4 * 0} + 0 * 2^{4 * 1} + 2 * 2^{4* 2} = 3 + 0 + 512 = 515}$ と展開される。

事前準備

スカラー倍算を行う点をPとした場合、ウィンドウ法では事前に、i = 1, ..., m-1について ${P_i = i * P}$ を計算する。m = 2⁴であれば、i = {1,...,15}について、 ${P_i = P_{i-1} + P}$ を計算する。つまりこの事前準備で2P, 3P, ..., 15Pが計算される。

計算

計算結果をQとし初期化(無限遠点Oと)する。
i = j〜0について以下の計算をループする。
1. Qの2倍算をw回実行し、結果をQにセットする。
2. ${x_i ＞ 0}$ の場合、事前計算したポイント ${x_iP}$ をQに加算し、結果をQにセットする。

w = 4, x = 515の場合、実際以下の計算が行われる。

最初のループ（j = 2）
1. Qの2倍算を4回実行し、結果をQにセット。この段階でQは無限遠点であるため、Qは無限遠点のまま。
2. Qに2Pを加算しQにセット。この段階でQ = 2P。
2回目のループ（j = 1）
1. Qの2倍算を4回実行し、結果をQにセット。この段階でQ = 32P。
3回目のループ（j = 0）
1. Qの2倍算を4回実行し、結果をQにセット。この段階でQ = 512P。
2. Qに3Pを加算しQにセット。この段階でQ = 515P。

ウィンドウ法もバイナリ法と同様、計算は2倍算と加算処理だが、点の加算回数はバイナリ法に比べて少なくなるという特徴があり、点の加算は2倍算に比べて遅いことからバイナリ法より計算コストが小さくなる。

射影座標を使った高速化

バイナリ法やウィンドウ法で計算される加算や2倍算の計算式を↑に記載したけど、どちらの計算でも逆元の計算（ $\left( \frac{y_{2} - y_{1}}{x_{2} - x_{1}} \right)$ および $\left( \frac{3x_{1}^{2} + a}{2y_{1}} \right)$ の部分）が登場する。この逆元の計算は拡張ユークリッド互助法などを使って行われるが、そのコストは他の加算や乗算と比べて非常に大きいので、逆元の計算を減らすことができれば計算速度も速くなる。

この逆元の計算を回避するために使われるのが射影座標になる。楕円曲線上の点は通常（x, y）座標で表すことが多く、これをアフィン座標と呼ぶ。この平面上の点の座標を、３つの要素（X, Y ,Z）で表現する射影座標と呼ばれる座標系に変換することで、計算コストの高い楕円曲線の除算の計算を回避でき、計算の高速化が可能になるという。

アフィン座標（x, y）は、射影座標（プロジェクティブ座標系）を使うと（x = X/Z mod p, Y/Z mod p）と表現でき、この変換を利用すると、アフィン座標系の楕円曲線の方程式

${y^{2} ≡ x^{3} + ax + b \mod p}$

を、射影座標系における楕円曲線の方程式

${Y^{2}Z ≡ X^{3} + aXZ^{2} + bZ^{3} \mod p}$

に変換できる。また射影座標においては以下が成立する。

２つの点(X, Y, Z), (X', Y', Z')について、X' = cX, Y' = cY, Z' = cZが成立する場合、2点は等しいとみなされる。
射影座標における無限遠点はX = 0かつZ = 0の点である（アフィン座標と違って射影座標では無限遠点も座標表記ができる）。

プロジェクティブ座標系における加算

２つのプロジェクティブ座標系の点 ${P = (X_1, Y_1, Z_1)}$ および ${Q = (X_2, Y_2, Z_2)}$ の加算 ${P + Q = (X_3, Y_3, Z_3)}$ （P != Q）の計算は、以下の式で計算される。

${u = Y_{2}Z_{1} - Y_{1}Z_{2}}$
${v = X_{2}Z_{1} - X_{1}Z_{2}}$
${A = u^{2}Z_{1}Z_{2} - v^{3} -2v^{2}X_{1}Z_{2}}$

として、

${X_{3} = vA}$
${Y_{3} = u(v^{2}X_{1}Z_{2} - A) - v^{3}Y_{1}Z_{2}}$
${Z_{3} = v^{3} Z_{1}Z_{2}}$

プロジェクティブ座標系における2倍算

P = Q、つまり2Pの計算は、以下の式になる。

${w = aZ_{1}^{2} + 3X^{2}_{1}}$
${s = Y_{1}Z_{1}}$
${B = sX_{1}Y_{1}}$
${h = w^{2} - 8B}$

として、

${X_{3} = 2hs}$
${Y_{3} = w(4B - h) - 8Y^{2}_{1}s^{2}}$
${Z_{3} = 8s^{3}}$

上記の計算式から逆元の計算が無くなっていることが分かる。そして算出した座標 ${(X_3, Y_3, Z_3)}$ をアフィン座標に変換 ${(x_3, y_3) = (X_3/Z_3, Y_3/Z_3)}$ する際に逆元の計算が発生するだけ。

ちなみに、 ${(x, y) → (X/Z^{2} \mod p, Y/Z^{3} \mod p)}$ の変換を行うのがヤコビ座標系。座標系としては、加算はプロジェクティブ座標系が速く、2倍算はヤコビ座標系が速いという特性があり、楕円曲線のスカラー倍算という意味ではヤコビ座標系を使った方が高速みたい。

参考

https://written.4403.biz/source/ecc_rev30r.pdf

と、長くなってきたので今回はここまで。

追記：

試しに射影座標を使ってスカラー倍算するコードをRubyで書いてみた↓

↑の同じバイナリ法でも、アフィン座標で計算するのに比べてかなり高速。逆元の計算が変換時以外発生いないだけでこうも違うのね。

Develop with pleasure!

福岡でCloudとかBlockchainとか。