nptclのブログ

Common Lisp処理系nptの開発メモです。https://github.com/nptcl/npt

AES-GCMの乗算を実装する

1. はじめに

共通鍵暗号であるAESの、 なんかよさそうな方式にGCMというものがあります。
こいつの乗算だけを実装します。

最初に言っておきますが理論は説明しません、というかできません。 知らないし。
ここで説明するのは、どうやって実装するかです。
特にNISTが配布しているPDFの乗算の最適化を詳しく見ていきます。

あとで詳しく言いますが、先に言っておかないとがっかりする人がいるかもしれないので、 これだけは先行して書いておきます。
作成する乗算の方法は下記の通り。

  • No tables
  • Simple, 8-bit tables
  • Simple, 4-bit tables
  • Simple, 1-bit tables　（PDFにないもの）

次のやつは作りません。

  • Shoup’s, 8-bit tables
  • Shoup’s, 4-bit tables

参考としてCommon Lispの完成品を下記にリンクしておきます。
AES, CCM, GCMが使えます。

https://github.com/nptcl/fixed/blob/main/aes/aes.lisp

1.1. 仕様書

まずは仕様書から。
たぶんこれが元だと思います。

もう一つあり、こいつは何なんだ？
検索で出てくるけどNIST内でリンク張ってるところが見つからなかった。

この下の方のPDFが最高なんです。

実装するという視点で言わせてもらいますが、 上の方は参考程度で、 下の方を見ながら作った方がいいと思います。
下は、encrypt, ghash, decryptをそれぞれ 数式で表したものがあるのですが、 それがすごく良かったです。
乗算以外はその式だけで作れます。
むしろそこだけ見たほうがわかりやすいです。

乗算はどちらのPDFを見ても作れると思います。
おすすめは下の方のPDFです。
さらに下のPDFは、乗算の最適化まで言及されています。
そいつをちゃんとやろう、というのがここでの目的です。

2. ガロア体の演算

ガロア体とは有限体の事だそうです。
GCMでは 2^{128}の有限体を扱うようで、 GF(2^{128})と書くそうです。
GCMではない、素のAESでは GF(2^{8})が使われていたのを覚えています。
楕円曲線ではさんざん素数 pの有限体をやってきました。
しかし今回は 2^{128}であり、 素数ではないので楕円曲線の時とはだいぶ事情が違います。
単に 2^{128}でひたすら割ってやるだけではダメみたいです。
これに気が付くのに時間がかかってしまった。

まず普通の値と、ガロア体の値の対応を考えます。
どちらも128bitの整数を扱いますが、 ビットのオーダーが逆転します。
普通の整数の最下位1bitは、ガロア体では最上位1bitになります。
なんでこんなことになるのかは知りませんが、 実装するならちゃんと覚えておかなければいけません。

覚えておいてほしいので詳しく書きます。
普通の値の、9を二進数で表してみます。

00000000 00000000 00000000 00000000
00000000 00000000 00000000 00000000
00000000 00000000 00000000 00000000
00000000 00000000 00000000 00001001

ちょっと長いですね。
扱う整数が128bitなので、こんなに長くなるわけです。
書きたくないので次のように省略します。

00000000 ... 00001001

左が上位ビットで、右が下位ビットです。
しかし、ガロア体の演算を行う際には逆転するので、 左が下位ビットになります。
これを覚えておいてください。

↓　ガロア体での最下位ビット
00000000 ... 00001001
                    ↑　普通の値の最下位ビット

変数 Yに値が整数が入っているとして、 ガロア体での i番目のビットを Y_{i}と記載します。
このとき、 Y_{0}は整数の最上位ビットであることに注意してください。
つまり、値が9の場合は Y_{0}=0です。
これも気が付くまで全然意味が分からなかった。
整数は128bitなので、 Y_{127}が最下位ビットになります。
つまり、値が9の場合は Y_{127}=1です。

ではガロア体に話を進めます。
ガロア体で使う演算は二種類。
加算と乗算です。

加算は単純なXORです。
普通の加算と違い、繰上りを考慮する必要がありません。
これは通常の加算よりだいぶ楽ですね。
XORは⊕と表します。
つまり、 Xと Yの加算は、 X \oplus Yです。

乗算はかなり面倒です。
乗算は X \cdot Yと表すことにします。
乗算の実装は、2進数のひっ算で計算すると思ってください。
ひっ算に必要な処理は、加算とシフトです。
シフトとはひっ算を進めていくうえで、 足す数の位を上げていく処理です。
このシフトを行った際に128bitを超えた場合は、 ガロア体特有の方法で値を修正する必要があります。
普通に考えれば128bitと128bitの乗算結果は256bit以下ですが、 今はガロア体ですので何とかして128bitに納めなければならないわけです。
素数のガロア体だったら乗算を求めた後に余りを求めればいいだけですが、 今回は素数ではないので簡単にはいきません。
単純に128bit切り捨てではだめです。
これが面倒、というか難しいわけです。

ここまでわかれば乗算のコードを理解できると思います。

3. 乗算コードの説明

乗算 Z = X \cdot Yの実装を示します。

Z ← 0, V ← X
for i = 0 to 127 do
  if Yのi番目のビット = 1 then
    Z ← ZとVの加算
  end if
  V ← Vの2倍
end for
return Z

このコードは、乗算を行うための汎用的なコードです。
ただの整数だろうがガロア体だろうが、これを行うと乗算できます。
やっていることはただの2進数のひっ算です。
Yのビットを一つずつ見ていき、 ビットが立っていたら結果に Xをシフトした値を加算します。

順番に見ていきます。
forの中の加算の部分を取り上げます。

if Yのi番目のビット = 1 then
  Z ← ZとVの加算
end if

先に Yの i番目のビットを Y_{i}と表すと書きました。
ということで、次のように書き直します。

if Yi = 1 then
  Z ← ZとVの加算
end if

for文の最初はi=0ですが、 Y_{0}は Yに格納されている整数の 最上位ビットであることに注意してください。
何度も書いていますが、ガロア体の演算ではビットの順番が逆転します。
Y0は、実際の値Yの最上位ビットです。
Y127は、実際の値Yの最下位ビットです。

コードの内容は、ビットが立っていたら Zに Vを加算するという内容です。
Vはビットの位置に応じて、 Xをシフトした内容です。
Zと Vの加算は、単純なXORですので、 とても簡単に書き直すことができます。

if Yi = 1 then
  Z ← Z ⊕ V
end if

この加算(XOR)は通常の加算とは異なり繰り上げがありませんので、 128bitを超えることがありません。
これで加算の部分は完了です。

問題は次の文です。

V ← Vの2倍

ここはちゃんと理解してください。
PDFには次のように記載されています。

if V127 = 0 then
  V ← rightshift(V)
else
  V ← rightshift(V) ⊕ R
end if

まず最初のV127というのは、 V_{127}のことで、 ガロア体の演算における Vの最上位ビットを表します。
つまり、格納されている整数の最下位ビットです。
最上位だの最下位だの、紛らわしいです。
そもそも何をしようとしているのかを説明します。

ここでは、 Vを2倍するためにシフトしています。
普通のプログラミングであれば左に1つシフトすることで2倍にすることができますが、 今はガロア体での演算のために右にシフトしています。
自分はこれに気が付くのに時間がかかってしまいました。
なんで右にシフトしてるの？ってずっと思ってた。

もしシフトして128bitを超えるような場合には、 ガロア体では単純に切り捨てることができないんです。
そこで、シフトで切り捨てられるビットをあらかじめ判定しておいて、 もしそれが1であったときは、 いきなり出てきた Rという定数のXORを取ります。
Rが何者かというと次の通り。

\displaystyle R = 11100001 \| 0^{120}

0xE1を120bit分だけ左にシフトしてください。
これはガロア体での話ではなく、 0xE1に対して普通に 2^{120}を乗算してくださいという意味です。
そういう定数のXORを取ることで、 はみ出した V_{128}=1のビットを何とかしてくれるんだそうです。
どうしてそうなるかはPDFに説明がありますのでご確認ください。
私にはわかりません。

これらを全部まとめると、 PDFに記載されているコードになります。
次のコードは引用です。

Z ← 0, V ← X
for i = 0 to 127 do
  if Yi = 1 then
    Z ← Z ⊕ V
  end if
  if V127 = 0 then
    V ← rightshift(V)
  else
    V ← rightshift(V) ⊕ R
  end if
end for
return Z

これが完成形です。
上記のコードをそのまま実装すれば正しく動きます。

Vを2倍することろに関してはもう少し説明が必要です。
2倍するという演算を、次のように表すことができるとのこと。

\displaystyle V \leftarrow V \cdot P

Pとの乗算の内容は、先ほど提示したコードにある通り、 右にシフトしてあふれた場合は RでXORを取るというものです。
この Pは、あとでいっぱい使うので覚えておいてください。

コードも次のように書き直せます。

Z ← 0, V ← X
for i = 0 to 127 do
  if Yi = 1 then
    Z ← Z ⊕ V
  end if
  V ← V ⋅ P
end for
return Z

乗算を計算するコードに Pの乗算が現れて混乱しますが、 別に再起呼び出しをするのではなく、 Pが単純に2倍を表すと考えていいと思います。

Pは複数回乗算できることを覚えておきましょう。
例えば 2^{8}倍するときは次のようになります。

\displaystyle V \leftarrow V \cdot P^{8}

4. 配列を用いた実装

乗算は、工夫をすることで早く実行できます。
まずは乗算がどこで使われるか見ていきましょう。
AES-GCMでは、GHASHを算出するときに使用されます。
PDFのGHASHの算出式を見てもらえばわかるのですが、 乗算はすべて次のような形になっています。

\displaystyle X_{i} = (...) \cdot H

カッコの中はいろいろだったので省略しました。
ここでのポイントは、必ず右が「 \cdot H」になっていることです。
乗算の片方は Hで固定なのです。
HはAESの鍵を設定さえすれば求めることができるので、 encrypt/decryptを行う前に Hを用いて最適化の準備をすることができます。

入力を Xとし、固定値 Hとの乗算 X \cdot Hを求めることを考えます。
その時に用いる最適化用の配列を Mとします。
配列にアクセスするための引数を「添字」と呼びます。
例えば、 M \lbrack 123 \rbrackと表したとき、 123が添字です。

4.1. 添字8bitの配列

添字8bitの配列とは、つまり要素数が256個の配列のことです。
こちらを用いて最適化することを考えていきましょう。

まずは Hと乗算したい入力 Xを1byteごとに分割します。
128bit=16byteなので、16個に分割することができます。
例えば次のような値であったとします。

X = 0xFFEEDDCCBBAA99887766554433221100

これを16個に分けます。

FF EE DD CC BB AA 99 88 77 66 55 44 33 22 11 00

これらは、あらかじめ用意しておいた16個の 配列 M_{0}, ..., M_{15}の要素を参照するのに使います。
Xと Hの乗算の結果を、次のように表します。

\displaystyle X \cdot H = M_{0} \lbrack FF \rbrack \oplus M_{1} \lbrack EE \rbrack \oplus ... \oplus M_{15} \lbrack 00 \rbrack

いま、 Xの各要素を FFだの EEだの定数を直打ちしていますが、 そうではなく変数 Xを1byteずつ分割した内容を x_{0} \ x_{1} \ ... \ x_{15}とします。
このとき乗算は次のように表せます。

\displaystyle X \cdot H = \bigoplus_{i=0}^{15} M_{i} \lbrack x_{i} \rbrack

配列 M_{0}, ..., M_{15}さえあれば、 乗算が求まるということです。
こんなことをして何がうれしいのかというと、 乗算をおよそ16回のXORで求めることができるので低コストなのです。

代償もあります。
データの容量、そして初期値の計算です。 いまは一つの配列 M_{0}の要素数は256個であり、 その中に128bitのデータが格納され、 その配列が16個あるという状況です。
事前に256 [個] × 16 [byte] × 16 [個] = 65536 [byte]の容量が必要になります。
さらに Hが決まった後で、 16個の配列を最初に計算しなければなりません。
そのコストは、256 [個] × 16 [個] = 4096 [個]分です。

4.2. 添字4bitの配列

4bitの配列を作ることを考えます。
8bitに比べると次の特徴があります。

  • 配列のサイズが小さい（ひとつあたり256個ではなく16個）
  • 初期値の設定が早い
  • 乗算の算出速度が遅い

いいところも悪いところあります。

まずは Xを4bitごとに分割します。
128bitなので、32個に分割することができます。
例えば次のように分割したとします。

F F E E ... 1 1 0 0

32個は多いので省略しました。
Xを4bitずつ分割した内容を x_{0} \ x_{1} \ ... \ x_{31}とします。
また、配列を M_{0}, ..., M_{31}とします。
このとき乗算は次のように表せます。

\displaystyle X \cdot H = \bigoplus_{i=0}^{31} M_{i} \lbrack x_{i} \rbrack

8bitのときは15回のXORで計算できましたが、 4bitのときは31回になっており、 およそ2倍のコストとなっています。
しかしその分スペースは小さくなります。

配列 M_{0}の要素数は16個であり、 その中に128bitのデータが格納され、 その配列が32個あります。
32 [個] × 16 [byte] × 16 [個] = 8192 [byte]の容量が必要になります。
さらに Hが決まった後で、 16個の配列を事前に計算しなければなりません。
そのコストは、32 [個] × 16 [個] = 512 [個]分です。

4.3. 添字1bitの配列

こちらはPDFにはないものですが、 とても簡単に実装できますし、 配列を用いない場合よりも効率が良くなります。
スペースは2048 [byte]必要になります。
初期値の計算は128個分。
乗算の演算はおよそ128回です。
乗算を求めるコードは配列を用いないときと似ていますが、 さらに単純になっています。

Z ← 0
for i = 0 to 127 do
  if Xi = 1 then
    Z ← Z ⊕ M[i]
  end if
end for
return Z

個人的には手軽でいいのではないかと思います。
合わせて見ていきます。

5. 配列の作り方

配列 M_{0}とそれに続く M_{n}なのですが、 作成方法があんまり詳しく乗ってないので苦労しました。
今回の投稿の主な目的は、この配列を作成する方法を示すことです。

5.1. 添字8bit配列の作り方

添字8bitの配列を作成しましょう。
値 Hをもとに、 M_{0}, ..., M_{15}を作成します。
各配列の要素数は、添字8bitということなので256個です。

まずは添字のbitがたった1つの場合の値を算出します。
値は128bitなので128通りありますが、 最初の配列 M_{0}分の8個を見ていきましょう。
具体的に書くと次の8通り。

10000000  添字128のとき
01000000  添字64のとき
00100000  添字32のとき
00010000  添字16のとき
00001000  添字8のとき
00000100  添字4のとき
00000010  添字2のとき
00000001  添字1のとき

これらの配列 M_{0}の値は、 Hの値をシフトしていくだけです。

  • 添字10000000は、 H
  • 添字01000000は、 H \cdot P
  • 添字00100000は、 H \cdot P^{2}
  • 添字00010000は、 H \cdot P^{3}
  • ...
  • 添字00000001は、 H \cdot P^{7}

H \cdot P^{2} = (H \cdot P) \cdot Pなので、 直前に算出した値を使うことができます。

では、次にbitがたくさんある場合を考えます。
例として24を考えましょう。 24の2進数表記は次のようになります。

00011000  24

つまりは、8と16の加算です。

00010000  16のとき
00001000  8のとき
----------------
00011000  24

値の求め方は簡単であり、8と16のそれぞれの配列の値を足すだけです。
つまり、

\displaystyle H \cdot P^{24} = (H \cdot P^{16}) \oplus (H \cdot P^{8})

右辺の2つの値はすでに算出していますので、 配列に格納されています。
すでに算出した値を用いて、 必要な値を全部求めることができます。

では、他の配列も合わせて見ていきます。
128bitを8bitに分割するので、配列は全部で16個。
配列を M_0, M_1, ..., M_{15}とします。
まずは、すべての配列において、bitが1つの場合を先行して求めてしまいます。

M_{0}の内容は下記の通り。

  • M_{0}
    • 10000000は、 H
    • 01000000は、 H \cdot P = M_0 \lbrack 128 \rbrack \cdot P
    • 00100000は、 H \cdot P^{2} = M_0 \lbrack 64 \rbrack \cdot P
    • 00010000は、 H \cdot P^{3} = M_0 \lbrack 32 \rbrack \cdot P
    • ...
    • 00000001は、 H \cdot P^{7} = M_0 \lbrack 2 \rbrack \cdot P

続いて M_{1}の配列を作成します。
内容は M_{0}の続きです。

  • M_{1}
    • 10000000は、 H \cdot P^{8} = M_{0} \lbrack 1 \rbrack \cdot P
    • 01000000は、 H \cdot P^{9} = M_{1} \lbrack 128 \rbrack \cdot P
    • 00100000は、 H \cdot P^{10} = M_{1} \lbrack 64 \rbrack \cdot P
    • 00010000は、 H \cdot P^{11} = M_{1} \lbrack 32 \rbrack \cdot P
    • ...
    • 00000001は、 H \cdot P^{15} = M_{1} \lbrack 2 \rbrack \cdot P

最初に M_{0}の値を参照している点に注意してください。
実装では直前に算出した値を Pでシフトして求めてください。
最後はこんな感じ。

  • M_{15}
    • 10000000は、 H \cdot P^{120} = M_{14} \lbrack 1 \rbrack \cdot P
    • 01000000は、 H \cdot P^{121} = M_{15} \lbrack 128 \rbrack \cdot P
    • 00100000は、 H \cdot P^{122} = M_{15} \lbrack 64 \rbrack \cdot P
    • 00010000は、 H \cdot P^{123} = M_{15} \lbrack 32 \rbrack \cdot P
    • ...
    • 00000001は、 H \cdot P^{127} = M_{15} \lbrack 2 \rbrack \cdot P

これで、bitが1つだけ立っている128個の要素は算出できました。

こちらを算出するコードですが、 PDFにAlgorithm 3として記載があります。
しかしこれがよくわからなかった。
PDFを参考に似たようなコードを作ってみます。

M0[128] = H
for i = 1 to 128 - 1 do
  j ← trunc(i / 8)
  k ← i % 8
  H ← H ⋅ P
  Mj[1 << (8 - k - 1)] ← H
end for

いくつかの命令を勝手に自作してしまいました。
truncは、小数点を切り捨てる命令です。
i % 8は、iを8で割ったあまりです。
1 << nという部分は1をnbit分左にシフトするという意味です。

他の要素はbitの組み合わせで全部出します。
まずは最初の M_{0}を完成させてみましょう。

i ← 2
while i < 256 do
  for j = 1 to i − 1 do
    M0[i + j] = M0[i] ⊕ M0[j]
  end for
  i ← 2 * i
end while

同じように、ほかの配列も作成できます。
しかし上記のコードを配列ごとに分ける必要はないので、 まとめてやってしまいましょう。
次のように書き換えます。

i ← 2
while i < 256 do
  for j = 1 to i − 1 do
    for k = 0 to 16 - 1 do
      Mk[i + j] ← Mk[i] ⊕ Mk[j]
    end for
  end for
  i ← 2 * i
end while

以上により、添字0以外はすべて完了です。
最後の値は下記の通り。

  • 00000000は、 0

コードは次の通り。

for k = 0 to 16 - 1 do
  Mk[0] ← 0
end for

完成です！

5.2. 添字4bit配列の作り方

添字4bitの配列を作成しましょう。
128bitを4bitに分割するので、配列は全部で32個。
配列を M_{0}, M_{1}, ..., M_{31}とします。
まずは、すべての配列において、bitが1つの場合を先行して求めてしまいます。

  • M_{0}

    • 1000は、 H
    • 0100は、 H \cdot P = M_0 \lbrack 8 \rbrack \cdot P
    • 0010は、 H \cdot P^{2} = M_0 \lbrack 4 \rbrack \cdot P
    • 0001は、 H \cdot P^{3} = M_0 \lbrack 2 \rbrack \cdot P

  • M_{1}

    • 1000は、 H \cdot P^{4} = M_{0} \lbrack 1 \rbrack \cdot P
    • 0100は、 H \cdot P^{5} = M_{1} \lbrack 8 \rbrack \cdot P
    • 0010は、 H \cdot P^{6} = M_{1} \lbrack 4 \rbrack \cdot P
    • 0001は、 H \cdot P^{7} = M_{1} \lbrack 2 \rbrack \cdot P

  • 以降省略

求め方は8bitの時とほぼ同じです。
実装を示します。

M0[8] = H
for i = 1 to 128 - 1 do
  j ← trunc(i / 4)
  k ← i % 4
  H ← H ⋅ P
  Mj[1 << (4 - k - 1)] ← H
end for

他の部分を求めます。

i ← 2
while i < 16 do
  for j = 1 to i − 1 do
    for k = 0 to 32 - 1 do
      Mk[i + j] ← Mk[i] ⊕ Mk[j]
    end for
  end for
  i ← 2 * i
end while

最後は0です。

for k = 0 to 32 - 1 do
  Mk[0] ← 0
end for

5.3. 添字1bit配列の作り方

こちらは8bit, 4bitと同じ内容ではあるものの、 素直に1bitの配列を作るわけではありません。
そんなことしても、どう考えても扱いづらいですから。
代わりに128bit分の、128個の配列を作ります。
各要素には、 Hをシフトした値を入れます。
つまりは次のようになります。

  • M
    • 0は、 H
    • 1は、 H \cdot P = M \lbrack 0 \rbrack \cdot P
    • 2は、 H \cdot P^{2} = M \lbrack 1 \rbrack \cdot P
    • 3は、 H \cdot P^{3} = M \lbrack 2 \rbrack \cdot P
    • ...
    • 127は、 H \cdot P^{127} = M \lbrack 126 \rbrack \cdot P

いくつか注意点があります。
まず、他とは違い、配列は Mひとつだけです。
添字0の値は、0ではなく Hです。
コードは簡単です。

for i = 0 to 128 - 1 do
  M[i] ← H
  H ← H ⋅ P
end for

6. 乗算の算出

では実際に算出していきましょう。

6.1. 添字8bit配列の乗算の算出

算出していきたいわけですが、 まず入力 Xを8bitに分割する必要があります。
これはプログラミング言語によって 実装方法が異なるのではないでしょうか。
疑似コードで無理やり作ってみます。

for i = 0 to 16 - 1 do
  r[i] ← (X >> (i * 8)) % 0x0100
end for

X >> nという部分はXをnbit分右にシフトするという意味です。
入力 Xの代わりに、求めた配列rを使い、乗算を計算します。

Z ← 0
for i = 0 to 16 - 1 do
  k ← r[16 - i - 1]
  Z ← Z ⊕ Mi[k]
end for
return Z

以上で、 X \cdot Hが求まりました。

6.2. 添字4bit配列の乗算の算出

さっそくコードを示します。

for i = 0 to 32 - 1 do
  r[i] ← (X >> (i * 4)) % 0x10
end for

Z ← 0
for i = 0 to 32 - 1 do
  k ← r[32 - i - 1]
  Z ← Z ⊕ Mi[k]
end for
return Z

6.3. 添字1bit配列の乗算の算出

こちらはすでに作成済みです。
同じものを下記に示します。

Z ← 0
for i = 0 to 127 do
  if Xi = 1 then
    Z ← Z ⊕ M[i]
  end if
end for
return Z

7. 乗算の速度

次の表は、AES-GCMのencrypt/decryptの実行結果です。
Common Lispでテストを行いました。
GHASHだけの計測ではないので注意してください。

Type Size [byte] Time [sec] process [cycle] cons [byte]
table0 0 33.7 111 G 14.7 M
table1 2048 27.2 89.7 G 6.68 M
table4 8192 25.0 82.6 G 5.38 M
table8 65536 24.3 80.5 G 4.67 M

Typeに対して、下記のコストを出しました。

  • Sizeは配列のサイズをbyteで表したものです。
  • Timeは、あるテストプログラムの実行時間（秒）です。
  • processはCPUの命令回数です。
  • consはLisp特有の値であり、実行時に生成されたconsの総数です。

上記の表でわかることは、配列 Mを使ったコードは ちゃんと速度において改善されているということです。
これはAES-GCMのencrypt/decryptすべての実行結果ですが、 GHASHだけのものがgcm-spec.pdfにありますので見てみましょう。

Method Storage requirement Throughput (cycles per byte)
No tables 16 bytes/key 119
Simple, 4-bit tables 8,192 bytes/key 17.3
Simple, 8-bit tables 65,535 bytes/key 13.1

この表を見ますと、いわゆるtable0とtable4のコストが119と17.3になっており、 差がとても大きいことがわかります。
配列を使った方が大幅に早いということになります。
この手の最適化はメモリと速度のトレードオフになると思いますが、 それにしても何らかの表を使った方がいいみたいです。

テストに使ったコードは下記の通り。

(defpackage work
  (:use common-lisp aes))
(in-package work)

(defconstant +size+ (* 10 1024 1024))

(defun make-vector8 (size)
  (make-array size :element-type '(unsigned-byte 8)))

(defun integer-big-vector (v size)
  (let ((a (make-array size :element-type '(unsigned-byte 8))))
    (dotimes (i size)
      (let ((k (- size i 1)))
        (setf (aref a i) (ldb (byte 8 (* k 8)) v))))
    a))

(defun main ()
  (let* ((x (make-vector8 +size+))
         (g (make-aes-gcm-128))
         (key (aes-gcm-key g))
         (nonce (aes-gcm-nonce g))
         a b)
    (dotimes (i 32)
      (setf (aref key i) (random #x0100)))
    (dotimes (i 12)
      (setf (aref nonce i) (random #x0100)))
    (dotimes (i +size+)
      (setf (aref x i) (random #x0100)))
    (multiple-value-bind (y a) (aes-gcm-encrypt g x)
      (multiple-value-bind (z b) (aes-gcm-decrypt g y)
        (unless (equalp a b)
          (format t "tab error, ~S /= ~S.~%" a b))
        (unless (equalp x z)
          (format t "decode error.~%"))))))

(let ((aes::*aes-gcm-mode* 'aes::table8))
  (time (main)))

C言語で楕円曲線DSAを実装する

C言語 Common Lisp

タイトルの通りC言語で楕円曲線DSAを作りました。

https://github.com/nptcl/fixed

作成した曲線は次の通り。

  • secp256k1
  • secp256r1
  • ed25519
  • ed448

Cのサンプルコードみたいな扱いだと思ってください。

以前、Common Lispで作成したもをそのままC言語にしています。
シリーズものですので、よろしかったらどうぞ。

Common Lispで楕円曲線DSAを実装する1 - nptclのブログ
Common Lispで楕円曲線DSAを実装する2（加算） - nptclのブログ
Common Lispで楕円曲線DSAを実装する3（乗算など） - nptclのブログ
Common Lispで楕円曲線DSAを実装する4（確認） - nptclのブログ
Common Lispで楕円曲線DSAを実装する5（鍵生成） - nptclのブログ
Common Lispで楕円曲線DSAを実装する6（署名） - nptclのブログ

コンパイルはたぶん簡単にできます。

$ make

あるいは、

$ cc *.c

で行けると思います。

乱数の初期化に/dev/randomを用いるので、 安全に実行したいならFreeBSDかLinux上で実行してください。
一応コンパイルオプションもあります。

$ cc *.c -DFIXED_FREEBSD
$ cc *.c -DFIXED_LINUX

実行すると次のような出力が出ます。

$ ./a.out
*** RSA
....
make_prime: 3
..
make_prime: 1
e: 10001
d: 9FFD3DDCA87DB013AD4163DBD671EF0FC12C46DD31D498315352511390CF8141
n: EEA9EA825D6BDE300EC2A44B8E2C4863AFF28B4ED97BEB34CBAEA1282F56ACE5
p: FA6A10ED0E527D2D5FEA69D823E22055
q: F3FCC03A690A91D03ADCB175C80AAA51
x1: 10
x2: 20
x3: 30
x4: 40
x5: 50
x1.encode: 7C09DE52BE6CDB831CE4718F524E46D338BDBD43D975794AA45860A8D0033436
x2.encode: 2A70611168D18974D12BDB03170BFFD0C9BB42D4016691C4EBAC8CCC122B2ECE
x3.encode: B76B913220960CD7F46940EFB5EFC354B43931CD102C8E11AEDB6E520D4A4275
x4.encode: C0BD68E1C79F60451B5B7A40245B8E793861B682372B365AA776AD6A12E835F
x5.encode: 5F78B535002AE5741BC650F52E22753B8A804F0DE1C590B1A4DECDCDF0A0FC18
x1.decode: 10
x2.decode: 20
x3.decode: 30
x4.decode: 40
x5.decode: 50


*** ECDSA, EdDSA
[private]
secp256k1: "AE988C663941C1BF51CF6FB9389B4ED787D65F33C5B2113F9B841C7F476FE501"
secp256r1: "C7A31EFEC292540504DD7BB2EC17AD5EA00D922AC628B9FD49F5C14C705F84B5"
ed25519  : "8ED9CDAB174D36EF6ABCC187F3EA4A92F390FB3543B180F2880E2659027BE2FF"
ed448    : "0C69F9CFFD0DED14D73A27AE7B31EC440F0E19FE2AD2D2E1B1D959493BB92ED79C08514EEA26752E182012538C21FDA59C45F10A093F222055"

[public]
secp256k1: "0395BD4A04A9FF7AE239FA42B7BA66B2E38E169DCA699AA1E25F6BAD2E5CDBBD33"
secp256r1: "0380FF1091B9436ED1D170B6E0B5F55C52A0525C61C2A487237FEB4AA7149C92A8"
ed25519  : "709F64FD1AC8DD64E5ADDEED6B6172B31D8BF0B589F303ED831FDF853D2C49CB"
ed448    : "27A40DA2BCAE9DCB16151162A8A652137C5C60B56A7BDFD63568AA68831335509D11EAAA131A9C324E95868F870CA0FCE685896BC49E97F000"

[sign] "Hello"
secp256k1.r: "57861B1A29BEB58F52596ABB825B8D3A8BAA638D70A8A5FAC0AA7ABEC16E443E"
secp256k1.s: "61711E5C9354B0113CB1290FB70D24112133038C2BDCE3DC911BAC03C1FB7186"
secp256r1.r: "D647CC95D5658F2AF8408CCE2A1869B4F7F0BD37FE7018DE67C8A9D6BCB68059"
secp256r1.s: "896B2F320EEB0ED2A0E079C215929E5B4E22D07FBD769338046F48C7CCFB50A2"
ed25519.r  : "77AEAFEB4D363B39D71B19A1C9585E78B6DCE1AC73F2EF01BA77E8AF67878B67"
ed25519.s  : "8DF76ADE07066909794E4F603CF8B3A1A12FC9437C868BD970713C7EC2E1B20D"
ed448.r    : "B694F15A8E7F5BFFE134876BEB8E4DE9F47CEB1F299174C1912FE940C6AF8A3908401AB65248E89F0A4C722468C4E508B6141E83C3DF5E5E80"
ed448.s    : "ED6956E23BCB50F5A611343F548EF5065758EF4B1E3F6AF42E6C15130ABF4B783F6459CFCA69004AFFFDAC8348B4D50338BFB0624E6C101100"

[verify]
secp256k1: T
secp256r1: T
ed25519  : T
ed448    : T

前半の*** RSAと、後半の*** ECDSA, EdDSAで別物です。
前半のRSAは以前記事で書きましたが、 巨大な素数を作って暗号・復号しています。

後半のECDSA, EdDSAは楕円曲線の 秘密鍵、公開鍵、あと"Hello"の署名と検証を行っています。
もし鍵を生成したいなら適当にどうぞ。

作った内容はCommon Lisp版と全く同じです。
C言語版では、鍵生成から署名、検証まで、 全部文字列だけでなんとなかるのが面白かったです。
例えば秘密鍵と公開鍵の生成は次のようになります。

#include "elliptic.h"
#include "random.h"
#include "signature.h"
#include <stdio.h>

int main(void)
{
    char private_key[200];
    char public_key[200];

    init_fixrandom();
    init_elliptic();
    private_string_secp256k1(private_key);
    public_string_secp256k1(private_key, public_key);
    printf("%s\n", private_key);
    printf("%s\n", public_key);

    return 0;
}

実行結果は下記の通り（一例）

E6BD585DCFB42001FD3DDBB0D53D7E1797210DF692567388A0D86D85571E7491
0322703DCC45AD4FBF96D04F662B82B1810B0A1C3E7D233BF570F3CFEC61B840B0

今回は配布物はなんの外部モジュールも必要とせず、 これ単体だけで動作するようになっています。
言い換えるならば必要なものは全部自分で作っています。
なので、副産物としていろいろできたので紹介しておきます。

あと、今回せっかく苦労して作ったので、 楕円曲線DSAのマイ公開鍵を公開しておきます。
それに何の意味があるのかはわかりませんが。

  • secp256k1
    • 03FEEF09658067CFBE3BE8685DDCE8E9C03B4A397ADC4A0255CE0B29FC63BCDC9C
  • secp256r1
    • 03CD92CF7B1C9CE9858383806B8540D72FB022BE577E21DE02B8EAA27371DB7AF2
  • ed25519
    • 75AB16F53A060E7AF9A4B8ECEA3D4DEF058AED2C626FEC96D5505C4A7D922960
  • ed448
    • 99AFC3768EE41B96F208EBAF8627908690DC6A5AC64659F93D0A46C2092B61E84AD14DD03F7B3F146799C29F65682126D517B7E1EA57716E00

これであってるのかな。
全然違ってたら悲しい。

Common Lispで楕円曲線DSAを実装する6（署名）

Common Lisp

前回の続きからです。
Common Lispで楕円曲線DSAを実装する5（鍵生成） - nptclのブログ

9. 署名と検証

DSAの機能である、署名と検証を行います。

9.1. 署名

署名の方法は、ECDSAとEdDSAで違います。
ECDSAは、署名の計算でメッセージ Mを一度しか使わないのに対して、 EdDSAは二度使います。
つまりECDSAはメッセージ Mをストリーミングできるのに、EdDSAはそれができません。
EdDSAはメッセージ Mをバッファリングする必要があると思います。

ちゃんと理由があって、ECDSAでは署名に際に乱数を使っているからです。
ところが「署名のたびに乱数を使うのは危険なのでやめませんか？」という流れがあったたようです。
いまこの実装をする気はありませんが、RFC6979に書いています。

Deterministic Usage of the Digital Signature Algorithm (DSA) and Elliptic Curve Digital Signature Algorithm (ECDSA)
https://datatracker.ietf.org/doc/html/rfc6979

EdDSAだと改良されています。
メッセージ Mのハッシュを余計にとることで、乱数の代用をしたわけです。
ということで、EdDSAはメッセージ Mを二度使います。
たぶんバッファリングが必要になりました。
安全のための犠牲です。
でもまあ検証ならともかく、署名にストリーミングは必要ないかな。

あと定数について注意してほしいことがあります。
署名と検証では、位数 nの \mod nがたくさん出てきます。
今までは有限体の計算ということで、 素数 pの \mod pばかりだったのですが、 今度は nが頻繁に出てくるので注意してください。

9.1.1. secp256k1, secp256r1

署名に必要なのは下記の情報です。

まずは自分の秘密鍵 qとは別に、秘密鍵と公開鍵のペア (k, B)を生成します。
生成する方法は説明済みですが、 kは新しく生成した乱数であり、 B = kGです。
kのことをナンス（nonce）というらしいです。
そうなんすね。

BのX座標を x_Bとしたとき、次の値 rを求めます。

\displaystyle r \equiv x_B \mod n

ここで nは座標に対する位数であり、曲線のパラメーターに載っています。
もし r=0のときはエラーなので、 kの乱数生成からやり直してください。

メッセージ Mのハッシュ値 Hを取得します。
ここで使用する関数はSHA-256です。

\displaystyle H = \text{SHA-256}(M)

得られたハッシュ値 HをBig Endianの整数 eとします。
下記の式より値 sを算出します。

\displaystyle s \equiv k^{-1} (e + r \ q) \mod n

この式は、素数 pではなく位数 nのあまりなので気を付けて下さい。
とくに k^{-1}の求め方は、 そもそも位数 n自体も素数なので、 式自体はいつものやつで問題ありません。

\displaystyle k^{-1} \equiv k^{n-2} \mod n

nも素数だったんですか。
知らなかった。

各変数の内容をまとめます。

  • qは、もともとの秘密鍵
  • kは、今回作成した秘密鍵
  • rは、今回作成した公開鍵のX座標（ただし \mod n）
  • eは、メッセージ MのSHA-256
  • nは、 Gに対する位数であり、曲線パラメータに乗っています

もし計算結果が s=0のときは良くないので、 秘密鍵 kの作成からやり直してください。
以上により、署名 (r, s)を求めることができました。

話題だけ出しておきますが、 secgのPDFには、 sはマイナスの値でも検証に合格してしまう みたいな話題がありました。
マイナスとは -s \mod nのです。
ものによっては、プラスマイナスの両方を算出して、 小さいほうを sにすることもあるそうです。
今は何もしません。

実装について説明します。
署名と検証で使うハッシュ関数は同じものを使うことができます。
データに対してBig Endianの整数を出力するものを定義します。

(defun sign-sha-secp256k1 (message)
  (let ((sha (make-sha256encode)))
    (read-sha256encode sha message)
    (vector-big-integer
      (calc-sha256encode sha))))

署名の実装は、 k選びに失敗したら 再度実行するという処理が必要になります。

(defun sign-loop-secp256k1 (private message)
  (let* ((k (make-private))
         (b (affine (make-public k)))
         (r (modn (point2-x b))))
    (unless (zerop r)
      (let* ((e (sign-sha-secp256k1 message))
             (s (modn (* (inverse-n k) (+ e (* r private))))))
        (unless (zerop s)
          (values r s))))))

(defun sign-secp256k1 (private message)
  (multiple-value-bind (r s) (sign-loop-secp256k1 private message)
    (if r
      (values r s)
      (sign-secp256k1 private message))))

9.1.2. ed25519

署名に必要なのは下記の情報です。

まずは秘密鍵 qから次の値を生成します。

  • ハッシュ値の下位32byteの整数、 uとする。（公開鍵生成に使用）
  • ハッシュ値の上位32byteの整数、 vとする。（署名に使用）

公開鍵で説明した通り、make-public-sign-ed25519関数で生成できます。
取得した uを用いて、公開鍵 Aを計算します。

\displaystyle A = uG
\displaystyle a = \text{encode}(A)

次の手順で整数 wを作成します。

  • SHA-512( v \| M)を計算する。
  • hash値64byteã‚’Little Endianの整数 w'とする。
  • w \equiv w' \mod nとなる wを計算する。

記号 \|は連結を表します。
nは位数であり、曲線のパラメーターに載っています。
いつもは素数 pの剰余を計算しますが、 今は pではなく nなので注意してください。

署名データである rを求めます。

\displaystyle R = w G
\displaystyle r = \text{encode}(R)

Gは起点であり、曲線パラメーターに載っています。
次の手順で整数 kを作成します。

  • SHA-512( r \| a \| M)を計算する。
  • hash値64byteã‚’Little Endianの整数 kとする。

記号 \|は連結を表します。
ここで再び Mが出てきます。
最初の Mと今回の Mがそれぞれ独立していればいいのですが、 そういうわけではないので Mのバッファリングなどが必要になります。

最後に sを求めます。

\displaystyle s \equiv (w + k \cdot u) \mod n

rと sの連結 r \| sが署名のデータです。

実装について説明します。
まずはハッシュ関数ですが、次の2通りのパターンがあります。

  • SHA-512( x \| M)
  • SHA-512( x \| y \| M)

そこで、3つの引数を取るハッシュ関数を作成しました。

(defun sign-sha-ed25519 (x y message)
  (let ((sha (make-sha512encode)))
    (when x (little-endian-sha512encode sha x 32))
    (when y (little-endian-sha512encode sha y 32))
    (read-sha512encode sha message)
    (modn (vector-little-integer
            (calc-sha512encode sha)))))

こちらは検証時にも使います。
署名の実装は下記の通り。

(defun sign-ed25519 (private message)
  (multiple-value-bind (u v) (make-public-sign-ed25519 private)
    (let* ((a (encode (multiple u *elliptic-g*)))
           (w (sign-sha-ed25519 v nil message))
           (r (encode (multiple w *elliptic-g*)))
           (k (sign-sha-ed25519 r a message))
           (s (modn (+ w (* k u)))))
      (values r s))))

9.1.3. ed448

ed25519とほぼ同じですが、ハッシュ関数だけ独特です。
署名と検証で使うハッシュ関数SHAKE-256-ed448を次のように定義します。

  • SHAKE-256-ed448(x)
    • sha3のshake-256で次のデータを読む
      • "SigEd448"、ASCII文字、合計8byte
      • 0x00, 合計1byte [phflag]
      • 0x00, 合計1byte [OLEN(context), context]
      • x, 引数のデータ
    • shake-256のハッシュ値を114byte返却

ハッシュ関数の実装を先行して説明します。

(defun sign-sha-ed448 (x y message &optional (sha1 0) (sha2 #()))
  (let ((sha (make-shake-256-encode)))
    (read-sha3encode sha (map 'vector #'char-code "SigEd448"))
    (byte-sha3encode sha sha1)
    (byte-sha3encode sha (length sha2))
    (read-sha3encode sha sha2)
    (when x (little-endian-sha3encode sha x 57))
    (when y (little-endian-sha3encode sha y 57))
    (read-sha3encode sha message)
    (modn (vector-little-integer
            (result-sha3encode sha 114)))))

それでは署名の説明を行います。
署名に必要なのは下記の情報です。

秘密鍵 qから次の値を生成します。

  • ハッシュ値の下位57byteの整数、 uとする。（公開鍵生成に使用）
  • ハッシュ値の上位57byteの整数、 vとする。（署名に使用）

make-public-sign-ed448関数で生成できます。
取得した uを用いて、公開鍵 Aを計算します。

\displaystyle A = uG
\displaystyle a = \text{encode}(A)

次の手順で整数 wを作成します。

  • SHAKE-256-ed448( v \| M)を計算する。
  • hash値114byteã‚’Little Endianの整数 w'とする。
  • w \equiv w' \mod nとなる wを計算する。

署名データである rを求めます。

\displaystyle R = w G
\displaystyle r = \text{encode}(R)

次の手順で整数 kを作成します。

  • SHAKE-256-ed448( r \| a \| M)を計算する。
  • hash値114byteã‚’Little Endianの整数 kとする。

最後に sを求めます。

\displaystyle s \equiv (w + k \cdot u) \mod n

rと sの連結 r \| sが署名のデータです。

署名の実装は下記の通り。

(defun sign-ed448 (private message)
  (multiple-value-bind (u v) (make-public-sign-ed448 private)
    (let* ((a (encode (multiple u *elliptic-g*)))
           (w (sign-sha-ed448 v nil message))
           (r (encode (multiple w *elliptic-g*)))
           (k (sign-sha-ed448 r a message))
           (s (modn (+ w (* k u)))))
      (values r s))))

9.2. 検証

検証を行います。

9.2.1. secp256k1, secp256r1

検証に必要なのは下記の情報です。

  • 公開鍵 A
  • メッセージ M
  • 署名 r\|s

公開鍵 Aは座標であることとします。
もしencodeされた値のときはdecodeしてください。

まずは署名 r, sが、どちらも1以上位数 n未満であることを確認します。
ゼロか n以上の場合は何かがおかしいので、失敗を返却して中断します。

メッセージ Mのハッシュ値 Hを取得します。
ここで使用する関数はSHA-256です。

\displaystyle H = \text{SHA-256}(M)

得られたハッシュ値 HをBig Endianの整数 eとします。
Mから eを計算するまでの方法は署名の時と同じです。

次の下記の式を求めます。

\displaystyle u_1 = e s^{-1} \mod n
\displaystyle u_2 = r s^{-1} \mod n

ここで注意してほしいのが、 sの逆数 s^{-1}です。
署名のときにも話題にしましたが、 今は素数 pのあまりではなく、位数 nのあまりです。
位数 n自体も素数なので、 s^{-1}は次の式で求めることができます。

\displaystyle s^{-1} \equiv s^{n-2} \mod n

座標 Rを求めます。

\displaystyle R = u_1 G + u_2 A

このとき、 R=Oなら失敗を返却して中断します。
RのX座標を x_Rとしたとき、

\displaystyle v \equiv x_R \mod n

を求めます。

検証の結果は次の通りです。

  • v = rなら検証成功
  • v \neq rなら検証失敗

実装は下記の通り。

(defun verify-secp256k1 (public message r s)
  (and (<= 1 r (1- *elliptic-n*))
       (<= 1 s (1- *elliptic-n*))
       (let* ((e (sign-sha-secp256k1 message))
              (s1 (inverse-n s))
              (u1 (modn (* e s1)))
              (u2 (modn (* r s1)))
              (p (addition
                   (multiple u1 *elliptic-g*)
                   (multiple u2 public))))
         (unless (zerop (point3-z p))
           (let* ((a (affine p))
                  (v (modn (point2-x a))))
             (= v r))))))

9.2.2. ed25519

検証に必要なデータは下記の通り。

  • 公開鍵 A
  • メッセージ M
  • 署名 r\|s

まずは署名データ sの範囲を確認します

  • n \le sなら不正な値なので検証失敗、中断

公開鍵 Aは座標であることとします。
公開鍵 Aのencodeを aとします。

次の手順で整数 kを計算します。

  • SHA-512( r \| a \| M)を計算する。
  • hash値64byteã‚’Little Endianの整数 kとする。

次に値 rをdecodeして、座標 Rを算出します。
得られた値 kと、 R, sにより、 次の式があっているかどうかを確かめます。

\displaystyle s G = R + k A

Gは起点であり、曲線のパラメーターに載っています。

式があっていたら検証は成功です。
あっていなかったら検証は失敗です。

実装は下記の通り。

(defun verify-ed25519 (public message r s)
  (when (< s *elliptic-n*)
    (let ((a (encode public))
          (p (decode r)))
      (when p
        (let* ((k (sign-sha-ed25519 r a message))
               (x (multiple s *elliptic-g*))
               (y (addition p (multiple k public))))
          (equal-point x y))))))

9.2.3. ed448

ed25519とほぼ同じです。
検証に必要なデータは下記の通り。

  • 公開鍵 A
  • メッセージ M
  • 署名 r\|s

署名データ sの範囲を確認します

  • n \le sなら不正な値なので検証失敗、中断

公開鍵 Aは座標であることとします。
公開鍵 Aのencodeを aとします。

次の手順で整数 kを計算します。

  • SHAKE-256-ed448( r \| a \| M)を計算する。
  • hash値114byteã‚’Little Endianの整数 kとする。

次に値 rをdecodeして、座標 Rを算出します。
得られた値 kと、 R, sにより、 次の式があっているかどうかを確かめます。

\displaystyle s G = R + k A

Gは起点であり、曲線のパラメーターに載っています。

式があっていたら検証は成功です。
あっていなかったら検証は失敗です。

実装は下記の通り。

(defun verify-ed448 (public message r s)
  (when (< s *elliptic-n*)
    (let ((a (encode public))
          (p (decode r)))
      (when p
        (let* ((k (sign-sha-ed448 r a message))
               (x (multiple s *elliptic-g*))
               (y (addition p (multiple k public))))
          (equal-point x y))))))

9.3. 署名と検証の確認

それでは、DSAの説明にあった状況をやってみましょう。
次の場合を想定します。

  • Helloが検証できたとき
  • Helloがなりすましで送信したとき
  • Helloが改ざんされたとき

Helloが検証できたときを実行します。
次のようなコードを実行します。

(defun verify-test ()
  (let* ((private (make-private))
         (public (make-public private))
         (msg1 (map 'vector #'char-code "Hello"))
         (msg2 (map 'vector #'char-code "Hello")))
    (multiple-value-bind (r s) (sign private msg1)
      (let ((v (verify public msg2 r s)))
        (format t "veriry: ~X~%" v)))))

(let ((*random-state* (make-random-state t)))
  (with-elliptic-secp256k1 (verify-test))
  (with-elliptic-secp256r1 (verify-test))
  (with-elliptic-ed25519 (verify-test))
  (with-elliptic-ed448 (verify-test)))

実行結果は下記の通り。

veriry: T
veriry: T
veriry: T
veriry: T

Helloがなりすましで送信したときとは、違う公開鍵で検証した場合です。
publicを作成するときに、privateに1を足してみます。

(defun verify-test ()
  (let* ((private (make-private))
         (public (make-public (mod (1+ private) *elliptic-p*)))
         (msg1 (map 'vector #'char-code "Hello"))
         (msg2 (map 'vector #'char-code "Hello")))
    (multiple-value-bind (r s) (sign private msg1)
      (let ((v (verify public msg2 r s)))
        (format t "veriry: ~X~%" v)))))

実行結果は失敗です。

Helloが改ざんされたときとは、メッセージが違っているときです。
msg2の内容を、Halloにしてみましょう。

(defun verify-test ()
  (let* ((private (make-private))
         (public (make-public private))
         (msg1 (map 'vector #'char-code "Hello"))
         (msg2 (map 'vector #'char-code "Hallo")))
    (multiple-value-bind (r s) (sign private msg1)
      (let ((v (verify public msg2 r s)))
        (format t "veriry: ~X~%" v)))))

実行結果は失敗です。

このようにして、署名と検証を行うことができました。

終わりです

長かった。
今度はこれをC言語で作ろうと思いますが、いったん終わり。

Common Lispで楕円曲線DSAを実装する5（鍵生成）

Common Lisp

前回の続きからです。
Common Lispで楕円曲線DSAを実装する4（確認） - nptclのブログ

7. DSAとは何か

あまりにも今さら過ぎますが、DSAの説明をします。
DSAとは、ディジタル署名アルゴリズム(Digital Signature Algorithm)だそうです。
ディジタル署名とは、「そのメッセージは、確かにあの人が書きましたよ」というのを確認するためのものです。

例として「私」が、みんなに向けてHelloというメッセージを 送信することを考えます。
Helloは大変重要なメッセージなので、 悪意がある人たちによって「私」のなりすましがあったり、 あるいはHelloが改ざんされてHalloになったりすることが懸念されます。

事前に私は、秘密鍵と公開鍵を作成します。
秘密鍵は、単なる乱数です。
秘密なので隠しておきましょう。

公開鍵は、起点 Gを秘密鍵やらhashやらを用いてスカラー倍したものです。
公開なので自分のページやメールで、 みんながわかる所に公開しておきましょう。

それでは署名を行います。
署名には下記のものが必要です。

署名によって次のものが得られます。

  • 署名 r \| s

rと sは、それぞれ32byteくらいのデータです。
rと sを連結したものを r\|sと表現し、それを署名データとします。

ではみんなにメッセージを送りましょう。
すでにみんなに周知されているものは下記になります。

  • 公開鍵

今回新たにみんなに送信するものは下記になります。

  • メッセージHello
  • 署名 r\|s

これで私の作業は終わりです。

場面が変わり、メッセージHelloを受け取ったみなさんの立場になります。
みなさんは、本当に私がHelloというメッセージを送信したのか疑っているわけです。
そこで次のデータをそろえます。

  • 公開鍵
  • メッセージHello
  • 署名 r\|s

これらのデータを使い、 署名が正しいかどうかを「はい」か「いいえ」で答えるというのが検証です。
「はい」ならHelloは「私」が書いたものです。
「いいえ」ならHelloは「私」が書いたものではありません。
なりすましか改ざんが行われたということです。

7.1 アリスとボブについて

何なのあいつら

7.2 曲線とエンディアン

ここからはDSAそのものを作成していくことになります。

自分が思ったことなのですが、 曲線ごとに使うエンディアンが違うようです。
作った人が違うからなのでしょう。
次のような感じを受けました。

  • secp256k1, secp256r1
    • Big Endianを使用（ただし未確認）
  • ed25519, ed448
    • Little Endianを使用（こちらは確定）

未確認だの確定だのあいまいですが、 secp256k1, secp256r1は自分が見てそう思った感想レベルです。
ただed25519とed448はRFCに出力例が載っており、 Little Endianとして扱わないとその値になりませんでした。

何にせよ、Endianを扱うための関数を用意する必要があります。
まずはLittle Endianから。

(defun integer-little-vector (v size)
  (let ((a (make-array size :element-type '(unsigned-byte 8))))
    (dotimes (i size)
      (setf (aref a i) (ldb (byte 8 (* i 8)) v)))
    a))

(defun vector-little-integer (v &key (start 0) end)
  (unless end
    (setq end (length v)))
  (let ((r 0) (k 0))
    (loop for i from start below end
          do
          (setq r (logior r (ash (aref v i) (* k 8))))
          (incf k 1))
    r))

上記の2つの関数は、次のような機能を持ちます。

  • integer-little-vector
    • æ•´æ•°ã‚’Little Endianと見なして、8byteの配列を返却
  • vector-little-integer
    • 8byteの配列をLittle Endianの列とみなして、整数を返却。

難しくはないと思います。

* (integer-little-vector #x12345678 4)
-> #(#x78 #x56 #x34 #x12)

* (vector-little-integer #(#x78 #x56 #x34 #x12))
-> #x12345678

Big Endianの方は次の通り。

(defun integer-big-vector (v size)
  (let ((a (make-array size :element-type '(unsigned-byte 8))))
    (dotimes (i size)
      (let ((k (- size i 1)))
        (setf (aref a i) (ldb (byte 8 (* k 8)) v))))
    a))

(defun vector-big-integer (v &key (start 0) end)
  (unless end
    (setq end (length v)))
  (let ((r 0) (k (- end start 1)))
    (loop for i from start below end
          do
          (setq r (logior r (ash (aref v i) (* k 8))))
          (decf k 1))
    r))

8. 鍵を生成

秘密鍵と公開鍵の生成を説明します。
それぞれの形式を覚えておいてください。

  • 秘密鍵は、整数です。
  • 公開鍵は、座標です。

たぶん一般的に鍵というのは、 16進数で記載されたテキストなんじゃないかなと思います。
でもここでは整数か座標として扱います。

8.1. 秘密鍵

最初に秘密鍵を作ります。

秘密鍵はただの乱数です。
しかし細かいところをちゃんと見ていきましょう。

8.1.1. 乱数を使うときの注意

暗号のための乱数は、非常に気を使った方がいいようです。
残念ながら今の私には/dev/urandomではなく /dev/randomを使って下さいくらいしか知識がありません。
しかもこれって確か乱数を得るためのものではなく、 乱数の初期値を設定するためのものだったと思う。

もし暗号に使う乱数を生成する場合、 疑似乱数生成で最良と考えられるMersenne Twisterですら 予測可能な乱数列なのでそのまま使うべきではないと言っています。

Mersenne Twister Home Page
よく聞かれる質問
http://www.math.sci.hiroshima-u.ac.jp/m-mat/MT/faq.html

SHAと組み合わせて下さいだそうです。
まじめにやるならよく考えましょう。
自分もとりあえずSHA-256あたりを絡めてみようと思います。
最近SHAを作ったので、使ってみたくなりました。

ということで、SHAを使った例を示します。

(defun make-private-256bit (&optional (n 4))
  (let ((hash (make-sha256encode)))
    (dotimes (i n)
      (little-endian-sha256encode hash (random (ash 1 256)) 32))
    (vector-little-integer
      (calc-sha256encode hash))))

256bitの乱数を、4回だけSHA-256に突っ込んでいます。
めんどうなら次の命令で代用できます。

(random (ash 1 256))

少しだけでもよくなっているのを期待します。
ed448の場合はSHA-512でも使ってください。

言い忘れていましたが、このためにSHAを自作しました。
好きに使って下さい。

https://github.com/nptcl/fixed/blob/main/sha.lisp

8.1.2. secp256k1, secp256r1

公式サイトのPDFを参考にします。

3.2.1 Elliptic Curve Key Pair Generation Primitive
https://www.secg.org/sec1-v2.pdf
Randomly or pseudorandomly select an integer d in the interval [1,n − 1].

つまり、1～n-1の乱数です。
256bitの乱数を出したら、nで割ってあまりを取りましょう。
0の場合はやり直しましょう。
あとあとでエラーになるっぽい。

★注意
乱数の範囲は、曲線パラメーターの n未満です。
素数 p未満ではないようです。
私は間違えました！

実装例を示します。

(defun modn (x)
  (mod x *elliptic-n*))

(defun make-private-secp256k1 ()
  (let ((x (modn (make-private-256bit))))
    (if (zerop x)
      (make-private-secp256k1)
      x)))

(defun make-private-secp256r1 ()
  (make-private-secp256k1))

(modn x)は、 x \mod nを計算するだけの関数です。
署名以降でたくさん使いますので、今のうちに用意しておきましょう。

8.1.3. ed25519

RFCを参考にします。

5.1.5. Key Generation
https://datatracker.ietf.org/doc/html/rfc8032
The private key is 32 octets (256 bits, corresponding to b) of cryptographically secure random data.

つまりは256bitの乱数そのものです。
nで割る必要はありません。

さらに言うと、どうせ秘密鍵のハッシュが使われるので、0でも問題ありません。

(defun make-private-ed25519 ()
  (make-private-256bit))

8.1.4. ed448

RFCを参考にします。

5.2.5. Key Generation
https://datatracker.ietf.org/doc/html/rfc8032
The private key is 57 octets (456 bits, corresponding to b) of cryptographically secure random data.

SHA-512で乱数を作ってみます。

(defun make-private-ed448 (&optional (n 4))
  (let ((hash (make-sha512encode)))
    (dotimes (i n)
      (little-endian-sha512encode hash (random (ash 1 512)) 64))
    (vector-little-integer
      (calc-sha512encode hash) :end 57)))

面倒なら次の命令で代用できます。

(random (ash 1 456))

あとLittle Endian関数を使っていますが、 乱数を乱数にするだけなのでどっちでもいいです。

8.2. 公開鍵

こちらは種類によって出し方が大きく変わります。
ただどれもこれも基本的には、起点 Gのスカラー倍が公開鍵です。

8.2.1. secp256k1, secp256r1

起点 Gを秘密鍵でスカラー倍したものが公開鍵です。
簡単ですね。

(defun make-public-secp256k1 (private)
  (multiple private *elliptic-g*))

(defun make-public-secp256r1 (private)
  (multiple private *elliptic-g*))

8.2.2. ed25519

EdDSAの場合はいろいろと作業があります。

5.1.5. Key Generation
https://datatracker.ietf.org/doc/html/rfc8032

全部Little Endianとして扱って下さい。
SHA-512の下位32byteとは0～31byteの事です。
上位32byteの32～63byteは署名の時に使います。

実装の話をします。
いきなり公開鍵まで作成してしまうのではなく、 次の2つの値を返却する関数を作成します。

実装は次のようになります。

(defun make-public-sign-ed25519 (private)
  (let ((sha (make-sha512encode)))
    (little-endian-sha512encode sha private 32)
    (let* ((v (calc-sha512encode sha))
           (a (vector-little-integer v :start 0 :end 32))
           (b (vector-little-integer v :start 32 :end 64)))
      (let ((v (ash 1 254)))
        (setq a (logand a (- v 8)))  ;; cofactor
        (setq a (logior a v)))
      (values a b))))

公開鍵だけを生成したいときは、 最初の返却値のみを使います。

(defun make-public-ed25519 (private)
  (multiple
    (make-public-sign-ed25519 private)
    *elliptic-g*))

参考情報ですが、 1～3bitを0にクリアするのはcofactorが8だからだそうです。
ed448はcofactorが4なので、2ビットしかクリアしません。
なんでそういうことするのかは全然知りませんが。

8.2.3. ed448

ed25519とは微妙に異なっていますが、だいたい同じです。

5.2.5. Key Generation
https://datatracker.ietf.org/doc/html/rfc8032

  • 秘密鍵をSHAKE-256によって、114byteのハッシュ値を出力する。
  • ハッシュ値の下位57byteã‚’æ•´æ•° sに変換
  • æ•´æ•° sに対して次のbit値をセット
    • 0 bitã‚’0
    • 1 bitã‚’0
    • 447 bitã‚’1
    • 448～455 bitã‚’0 （最後の1byte=8bitを、全部0にする）
  • 起点 Gのスカラー倍 sを公開鍵とする。
(defun make-public-sign-ed448 (private)
  (let ((sha (make-shake-256-encode)))
    (little-endian-sha3encode sha private 57)
    (let* ((v (result-sha3encode sha 114))
           (a (vector-little-integer v :start 0 :end 57))
           (b (vector-little-integer v :start 57 :end 114)))
      (let ((v (ash 1 447)))
        (setq a (logand a (- v 4)))  ;; cofactor
        (setq a (logior a v)))
      (values a b))))

(defun make-public-ed448 (private)
  (multiple
    (make-public-sign-ed448 private)
    *elliptic-g*))

8.3. 座標のencode

公開鍵ができたので、事前にみんなに周知したくなります。
しかし公開鍵は座標です。
座標をそのまま誰かに送ったら、いや、そんなのわからんし、となると思います。
そこでわかりやすい形式に変換する、encodeを行いましょう。

secp256k1とsecp256r1は、まさにこのような目的で人のためにencodeします。
ed25519とed448は、それに加えて署名の処理でencodeを使います。
曲線によって用途が微妙に違います。
人に見て欲しいだけの場合、encodeは配列を返却することにしました。
一方、その後の署名の処理でも使う場合、encodeは整数を返却します。

  • encodeの返却値は、
    • secp256k1とsecp256r1は、byteの配列
    • ed25519とed449は、整数

encodeの内容を見ていきましょう。
座標を A、同等の意味をもつ値を k（配列か整数）としたとき、

  • 座標 Aから、 kへの変換が、encode
  • kから、座標 Aへの変換が、decode

となります。

座標 Aはアフィン座標にすることで x, yそれぞれ32byteくらいのデータになります。
ただ、この座標はたぶん曲線上に位置しているので、 片方あれば算出することができます。
どちらかの値だけを保存することで短くすることができることを覚えておきましょう。

8.3.1. secp256k1, secp256r1

ECDSAのencodeは、次の場所で定義されています。

2.3.3 Elliptic-Curve-Point-to-Octet-String Conversion
https://www.secg.org/sec1-v2.pdf

このPDFにはencode/decodeという語は出てきません。
また、ECDSAだけの話になりますが、 compression, uncompressionという 2つの表現方法があります。
compressionは、X座標だけを表します（全33byte長）。
uncompressionは、X座標の次にY座標を表します（全65byte長）。
compressionかuncompressionはただの選択肢であり、 どっちでもいいので好きな方を選んでください。
それでは順にみていきましょう。

座標 Aをencodeすることを考えます。
もし A=Oなら、次の1byteで終わりです。

#(0x00)

判定方法は、 Z=0のときで問題ありません。

生成コードは次の通り。

(integer-big-vector #x00 1)

A \neq O`であるときを考えます。
座標 Aをアフィン座標に変換して話を進めます。

最初にuncompressの場合を見ていきます。
次の順番でデータを連結します。

  • 0x04の1byte
  • Big Endianの x値
  • Big Endianの y値

例えばx=0x0100, y=0x0200であり、 例文のために32bit長としたときは次のようになります。

#(04 00 00 01 00 00 00 02 00)

返却値は配列であり、要素0番目が左(04の値)になります。

生成コードは次の通り。

(defun encode-uncompress-secp256k1 (x y)
  (integer-big-vector
    (logior (ash #x04 (* 256 2)) (ash x 256) y)
    (1+ 64)))

続いてcompressを示します。
最初の1byteは、 yの最下位ビットによって変わります。 次の順番でデータを連結します。

  • yの最下位ビットが、
    • 0のときは0x02
    • 1のときは0x03
  • Big Endianの x値

例えばx=0x00000100, y=0x00000200であり、 32bit長のときは次のようになります。

#(02 00 00 01 00)

他の例として、x=0x00000100, y=0x00000201であり、 32bit長のときは次のようになります。

#(03 00 00 01 00)

生成コードは次の通り。

(defun encode-compress-secp256k1 (x y)
  (integer-big-vector
    (logior (ash (if (logbitp 0 y) #x03 #x02) 256) x)
    (1+ 32)))

8.3.2. ed25519

座標 Aをアフィン座標 (x,y)にします。
エンコードされた整数 kは下記のように作成されます。

  • k \leftarrow y
  • xの最下位ビットを kの最上位ビット(255bit)にコピー

整数 kは256bit長(32byte)の整数になります。
基本は yの値だけを保存するのですが、 他に xの符号情報も必要だとのこと。
xの符号は最下位ビットです（最上位ではなく）。

生成コードは次の通り。

(defun encode-ed25519 (v)
  (let* ((a (affine v))
         (x (point2-x a))
         (y (point2-y a)))
    (when (logbitp 0 x)
      (setq y (logior y (ash 1 255))))
    y))

返却値 kは、署名でも使うので整数です。
ところが公開鍵としてみんなに周知したいこともあると思います。
そんなときは、little endianで32byteの配列に直してください。
つまり、例として256bitは長すぎるので32bitで示しますが、

#x12345678 -> $(#x78 #x56 #x34 #x12)

のようになります。

8.3.3. ed448

ed25519とほぼ同じですが、特徴があります。
まずやり方から。

  • k \leftarrow y
  • xの最下位ビットを kの455bitにコピー

生成コードは次の通り。

(defun encode-ed448 (v)
  (let* ((a (affine v))
         (x (point2-x a))
         (y (point2-y a)))
    (when (logbitp 0 x)
      (setq y (logior y (ash 1 455))))
    y))

8.4. 座標のdecode

encodeされたデータをdecodeして座標に戻しましょう。
すでに見てきた通り、曲線によってencodeの返却値の形式が違うので、 decodeの引数の形式も、それに合わせる必要があります。
次のようになります。

  • decodeの引数は、
    • secp256k1とsecp256r1は、byteの配列
    • ed25519とed449は、整数

返却値はどちらも座標です。

8.4.1. secp256k1, secp256r1

引数のbyteの配列 kから、座標を求めます。

次の場所で定義されています。

2.3.4 Octet-String-to-Elliptic-Curve-Point Conversion
https://www.secg.org/sec1-v2.pdf

kの最初の1 byteを見て形式を判定します。
こんな感じで実装しました。

(defun decode-secp256k1 (k)
  (case (aref k 0)
    (#x00 (make-point3 0 0 0))
    (#x02 (decode-compress-secp256k1 k 0))
    (#x03 (decode-compress-secp256k1 k 1))
    (#x04 (decode-uncompress-secp256k1 k))))

もし最初のbyteが0x00なら、座標は無限遠点の Oです。

もし最初のbyteが0x04なら、compressionではない座標データです。
続く32byteを、Big Endianの整数として xの値にします。
続く32byteを、Big Endianの整数として yの値にします。
得られた値から、射影座標 (x, y, 1)で返却します。

配列 kからそのまま値を取り出しましょう。
値がvalidかどうかは判定したほうがいいと思います。

(defun decode-uncompress-secp256k1 (k)
  (let ((p (make-point3
             (vector-big-integer k :start 1 :end 33)
             (vector-big-integer k :start 33 :end 65))))
    (when (valid p)
      p)))

もし最初のbyteが0x02か0x03なら、compressionです。
この値は yの符号を表しています。
次のように符号 y_0を決めてください。

  • 最初の1byteが0x02なら、 y_0=0
  • 最初の1byteが0x03なら、 y_0=1

続く32byteを、Big Endianの整数として xの値にします。
yの値そのものは格納されていませんので、次の手順で算出する必要があります。

次の式より \alphaを求めます。

\displaystyle \alpha \equiv x^{3} + a x + b \mod p

ここで、 aと bは曲線のパラメーターです。

次に \alphaの平方根 \betaを求めます。
平方根の求め方はすでに説明済みです。

最後に符号を合わせます。
算出した平方根 \betaの最下位ビットを \beta_0としたとき、

  • \beta_0 = y_0なら、 y = p - \beta
  • \beta_0 \neq y_0なら、 y = \beta

以上で、 xと yが求まりました。
射影座標 (x, y, 1)がDecodeの返却です。

実装は次の通りです。

(defun decode-compress-secp256k1 (k y0)
  (let* ((x (vector-big-integer k :start 1 :end 33))
         (a (modp (+ (* x x x) (* *elliptic-a* x) *elliptic-b*)))
         (y (square-root-mod-4 a)))
    (when y
      (unless (= (logand y #x01) y0)
        (setq y (- *elliptic-p* y)))
      (make-point3 x y))))

8.4.2. ed25519

引数の整数 kから、座標を求めます。

まずは kから次の値を算出します。

  • kの0～254bitã‚’ yとする。
  • kの255bitã‚’ xの符号 x_0とする。

encodeの逆の処理をしているだけなのでわかると思います。
求めた yのチェックを行います。

  • yが p以上の場合はdecode失敗

次に yと x_0を用いて、 xを求めましょう。
xの算出は、曲線の式を用います。
次の u, vを求めます。

\displaystyle u = y^{2} - 1
\displaystyle v = d \ y^{2} + 1

この値を用いて仮決めの xを求めます。

\displaystyle x \equiv u \ v^{3} (u \ v^{7}) ^{(p-5)/8} \mod p

確認のための値 v \ x^{2}を計算します。

  1. v \ x^{2} \equiv u \mod pなら、 xが平方根
  2. v \ x^{2} \equiv -u \mod pなら、平方根 xは、 x \leftarrow x \cdot 2^{(p-1)/4}
  3. それ以外なら平方根はなし（decode失敗）

符号 x_0で場合分けをします。

  • もし x=0でありかつ x_0=1のときは、decodeは失敗
  • もし xの最下位ビットと x_0が等しくないなら、 x \leftarrow p - xを代入

結果の (x, y)がdecodeの返却値です。

実装は、式の計算があり少々複雑です。

(defun decode-ed25519-x (y)
  (when (< y *elliptic-p*)
    (let* ((yy (* y y))
           (u1 (modp (1- yy)))
           (v1 (1+ (* *elliptic-d* yy)))
           (v2 (* v1 v1))
           (v3 (* v1 v2))
           (v4 (* v2 v2))
           (uv3 (* u1 v3))
           (p (/ (- *elliptic-p* 5) 8))
           (x (mulp uv3 (power-mod (* uv3 v4) p *elliptic-p*)))
           (v1x2 (mulp v1 x x)))
      (cond ((= v1x2 u1) x)
            ((= v1x2 (- *elliptic-p* u1))
             (mulp x (power-mod 2 (/ (1- *elliptic-p*) 4) *elliptic-p*)))))))

(defun decode-ed25519 (k)
  (let* ((y (ldb (byte 255 0) k))
         (x0 (ldb (byte 1 255) k))
         (x (decode-ed25519-x y)))
    (cond ((null x) nil)
          ((and (= x 0) (= x0 1)) nil)
          ((/= (logand x #x01) x0) (make-point4 (- *elliptic-p* x) y))
          (t (make-point4 x y)))))

7.3.4. ed448

流れはed25519と同じです。
kから次の値を算出します。

  • kの0～454bitã‚’ yとする。
  • kの455bitã‚’ xの符号 x_0とする。

求めた yのチェックを行います。

  • yが p以上の場合はdecode失敗

u, vを求めます。

\displaystyle u = y^{2} - 1
\displaystyle v = d \ y^{2} - 1

ed25519とは符号が違うので注意。
仮決めの xを求めます。

\displaystyle x \equiv u^{3} \ v (u^{5} \ v^{3}) ^{(p-3)/4} \mod p

確認のための値 v \ x^{2}を計算します。

  1. v \ x^{2} \equiv u \mod pなら、 xが平方根
  2. それ以外なら平方根はなし（decode失敗）

符号 x_0で場合分けをします。

  • もし x=0でありかつ x_0=1のときは、decodeは失敗
  • もし xの最下位ビットと x_0が等しくないなら、 x \leftarrow p - xを代入

結果の (x, y)がdecodeの返却値です。

(defun decode-ed448-x (y)
  (when (< y *elliptic-p*)
    (let* ((yy (* y y))
           (u1 (modp (1- yy)))
           (u2 (* u1 u1))
           (u3 (* u1 u2))
           (u5 (* u2 u3))
           (v1 (1- (* *elliptic-d* yy)))
           (v3 (* v1 v1 v1))
           (u3v1 (* u3 v1))
           (u5v3 (* u5 v3))
           (p (/ (- *elliptic-p* 3) 4))
           (x (mulp u3v1 (power-mod u5v3 p *elliptic-p*)))
           (v1x2 (mulp v1 x x)))
      (when (= v1x2 u1)
        x))))

(defun decode-ed448 (k)
  (let* ((y (ldb (byte 455 0) k))
         (x0 (ldb (byte 1 455) k))
         (x (decode-ed448-x y)))
    (cond ((null x) nil)
          ((and (= x 0) (= x0 1)) nil)
          ((/= (logand x #x01) x0) (make-point3 (- *elliptic-p* x) y))
          (t (make-point3 x y)))))

続きます

次は署名と検証をやります。
Common Lispで楕円曲線DSAを実装する6（署名） - nptclのブログ

Common Lispで楕円曲線DSAを実装する4（確認）

Common Lisp

前回の続きからです。
Common Lispで楕円曲線DSAを実装する3（乗算など） - nptclのブログ

6. 演算のテスト

ここまで来たらちゃんとテストします。
次の機能が正常に機能することを確認しましょう。

演算の内容は、使用する曲線によって全部異なりますので、 地道に全通りやる必要があります。
下記の4種類を全部やります。

  • secp256k1
  • secp256r1
  • ed25519
  • ed448

6.1. 位数について

その前に、いまいち理解できてない「位数」の話題をします。
乗算のテストの時に使います。

位数とは、どうも2つ意味があるようです。

使用する曲線のことを Eとか言ったりします。
そのとき \# Eとは、 Eの位数といい、 曲線上に存在できる点の数を表します。
曲線の位数 \# Eは、少なくとも実装するときには全く使いません。

一方、曲線パラメーターにある nは、起点 Gに属する位数だそうです。
次のようになります。

\displaystyle n G = O

つまりは加算を n回くらい繰り返すと 元に戻ってきてぐるぐる回るということなのでしょう。
場合によって nは、 Lとか qとか言われたりするようです。
また、次の関係が成り立つそうです。

\displaystyle \# E = n \cdot h

ここで hは曲線パラメーターに載っており、 cofactorとか言われてます。
何者かは知りません。

参考までに、secp256k1の位数 nは次のようになります。

n = 0xfffffffffffffffffffffffffffffffebaaedce6af48a03bbfd25e8cd0364141

これ一体いくつなんだ。

でもなんとなくわかりました。
テストで使うのは次の式です。

\displaystyle n G = O

これでaddition, doubling, multipleが 正しく実装されているかどうかを確認します。

6.2. 曲線上の座標かどうかを確認

この確認はめちゃ便利です。

ある座標が、ちゃんと曲線上にあるのかどうかを確認しましょう。
確認の方法は簡単で、曲線の式にそのまま代入するだけです。
しかしアフィン座標に変換する必要があることだけは注意してください。

例としてsecp256k1の確認をしてみましょう。
確認する座標 Aを次のように表します。

\displaystyle A = (X, Y, Z)

まずはアフィン座標に変換します。

\displaystyle x = X Z^{-1}
\displaystyle y = Y Z^{-1}

アフィン座標を用いて、次の式が成り立っているかを確認します。

\displaystyle y^{2} = x^{3} + a x + b

ここで、 aと bは、曲線のパラメーターです。
確認は簡単ですが、これを4つの曲線で作らなければいけないのが厄介です。
曲線の式は2種類しかありませんので、 とっととやってしまいましょう。

ed25519のみの話ですが、拡張された射影座標の場合は、 ついでに Tも確認しましょう。
どういうことかというと、

\displaystyle A = (X, Y, Z, T)

のとき、曲線上にあるかどうかを確認した後、

\displaystyle xy = T Z^{-1}

も併せて確認します。
加算の際に Tは独立して計算されるので、 確認する意味はあります。

完成したら起点 Gでやってみるといいと思います。
絶対にあっているはずです。

6.2.1. 座標のチェック

確認する内容は大きく2つ。

  • 曲線上に存在すること
  • T値が正しいこと

まずは曲線上であること確認します。
曲線は2つ使っています。

  • Weierstrass曲線
    • secp256k1
    • secp256r1
  • Edwards曲線
    • ed25519
    • ed448

その2種類の曲線の確認を行います。
方法は簡単で、ただ式に代入するだけです。
実装します。

Weierstrass曲線の確認コード

(defun valid-weierstrass (v)
  (let* ((a (affine v))
         (x (point2-x a))
         (y (point2-y a)))
    (zerop
      (modp (- (* y y)
               (+ (* x x x) (* *elliptic-a* x) *elliptic-b*))))))

Edwards曲線の確認コード

(defun valid-edwards (v)
  (let* ((a (affine v))
         (x (point2-x a))
         (y (point2-y a))
         (xx (* x x))
         (yy (* y y)))
    (zerop
      (modp (- (+ (* *elliptic-a* xx) yy)
               (+ 1 (* *elliptic-d* xx yy)))))))

6.2.2. T値のチェック

これは Tを持つ、ed25519でのみ行われる確認です。
Tには次のような関係があります。

\displaystyle x \cdot y = T/Z

単純に確認します。

(defun valid-point4 (v)
  (let* ((z (inverse (point4-z v)))
         (x (mulp (point4-x v) z))
         (y (mulp (point4-y v) z))
         (xy (mulp (point4-xy v) z)))
    (zerop
      (modp (- (* x y) xy)))))

(defun valid-ed25519 (v)
  (and (valid-edwards v)
       (valid-point4 v)))

6.2.3. Common Lispでの実装

確認の内容は曲線によって異なります。
次のものを使用します。

  • secp256k1, secp256r1
    • valid-weierstrass関数を使用
  • ed25519
    • valid-ed25519関数を使用
  • ed448
    • valid-edwards関数を使用

曲線ごとに場合分けをする場合、 additionなどと同じようにspecial変数を使います。

(defvar *elliptic-valid*)

(defun valid (v)
  (funcall *elliptic-valid* v))

束縛する内容は次のようになります。

;; secp256k1, secp256r1
(let ((*elliptic-valid* #'valid-weierstrass))
  ...)
;; ed25519
(let ((*elliptic-valid* #'valid-ed25519))
  ...)
;; ed448
(let ((*elliptic-valid* #'valid-edwards))
  ...)

6.3. 射影座標で同じ場所かどうか確認

2つの射影座標 A, Bが同じ場所かどうかを確認するにはどうしたらいいでしょうか。
Zがあるので、 Xと Yだけを比較しても 等しいかどうかがわからないのです。

それぞれをアフィン座標に変換して比較する必要があります。
ただ、RFCのPythonのコードでは次のようなものがありました。

def point_equal(P, Q):
    # x1 / z1 == x2 / z2  <==>  x1 * z2 == x2 * z1
    if (P[0] * Q[2] - Q[0] * P[2]) % p != 0:
        return False
    if (P[1] * Q[2] - Q[1] * P[2]) % p != 0:
        return False
    return True

つまりは逆数を求めずとも、乗算だけで調べられるようです。
速度が必要な場合は便利なので実装しましょう。

一応日本語で手順を示します。 座標 Aと Bを次のようにします。

\displaystyle A = (X_A, Y_A, Z_A)
\displaystyle B = (X_B, Y_B, Z_B)

このとき、次の判定を順番に行ってください。

  1. X_A Z_B \not\equiv X_B Z_A \mod pのとき、等しくない
  2. Y_A Z_B \not\equiv Y_B Z_A \mod pのとき、等しくない
  3. 上記が当てはまらない場合は、等しい

実装例を示します。

(defun equal-point (p q)
  (let ((pz (point3-z p))
        (qz (point3-z q)))
    (and (zerop (modp (- (* (point3-x p) qz) (* (point3-x q) pz))))
         (zerop (modp (- (* (point3-y p) qz) (* (point3-y q) pz)))))))

アフィン座標に直して比較するのでもいいと思います。
わかりやすいですが、めちゃくちゃ非効率です。

(defun equal-point-affine (x y)
   (let ((a (affine x))
         (b (affine y)))
     (and (eql (point2-x a) (point2-x b))
          (eql (point2-y a) (point2-y b)))))

6.4. 演算の確認

addition, doubling, multipleの確認を行いましょう。
確認する内容は次のとおり。

  • addtion O+Gが Gであることをチェック
  • addtion G+Oが Gであることをチェック
  • doubling Oが Oであることをチェック
  • doubling Oにadditionを実行してチェック
  • addtion G+Gが曲線上にあるかチェック
  • doubling 2Gが曲線上にあるかチェック
  • multiple 2Gが曲線上にあるかチェック
  • 上記3つの座標が等しいかを確認

スカラー倍もできるところまで確認しましょう。

  • multiple 5Gと、 G+G+G+G+Gの確認
  • multiple 5Gと、 2G+2G+Gの確認
  • \text{multiple} (nG) = Oの確認
  • \text{multiple} (nG) + G = Gの確認
  • \text{multiple} (1000 G) + \text{multiple} (n - 1000) = Oの確認

ここまで出来たら、ひとまずは問題ないでしょう。
もし問題があった場合は、署名のあとの検証の時に引っかかると思います。

あと、座標の確認に次のコードを使います。

(defun print-point (v &optional s)
  (let ((s (or s t)))
    (format s "X: ~X~%" (point3-x v))
    (format s "Y: ~X~%" (point3-y v))
    (format s "Z: ~X~%" (point3-z v))
    (format s "  -> Valid: ~X~%" (valid v))))

ただ座標を標準出力に出すだけです。

6.4.1. Oの確認

下記2つは、何事もなくOKでした（結果は省略）。

  • addtion O+Gが Gであることをチェック
  • addtion G+Oが Gであることをチェック

ただし覚えておかなければいけないことがあり、 Oに対してvalidの返却値が異なっています。

  • secp256k1では、(valid O)はnil
  • secp256r1では、(valid O)はnil
  • ed25519では、(valid O)はt
  • ed448では、(valid O)はt

ちゃんと理由があって、アフィン座標に変換すると

  • secp256k1では、 (0,0)
  • secp256r1では、 (0,0)
  • ed25519では、 (0,1)
  • ed448では、 (0,1)

になるので、それが曲線上にあるかどうかを判定しているだけです。
この結果は覚えておいてください。

下記について

  • doubling Oが Oであることをチェック
  • doubling Oにadditionを実行してチェック

secp256k1, secp256r1は、 O = (1,0,0)のようにXに何かが入っていたとき、 doublingを取るとYに何らかの値が入ります。

内容を出力してみます。

X: 0
Y: FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEFFFFF56F
Z: 0
  -> Valid: NIL
X: 0
Y: FFFFFFFF00000001000000000000000000000000FFFFFFFFFFFFFFFFFFFFF93F
Z: 0
  -> Valid: NIL

しかし Z=0なので、これは問題なし。

ed25519, ed448は、そもそも値が変化します。

X: 0
Y: 7FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEC
Z: 7FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEC
  -> Valid: T
X: 0
Y: FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFE
Z: FFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFEFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFFE
  -> Valid: T

これは正しいのでしょうか。
アフィン座標は正しそうですが、加算すればどうなるかやってみます。

確認するコードは次の通り。

(defun check-equal2 ()
  (let* ((g *elliptic-g*)
         (o (doubling (make-point4 0 1 1 0)))
         (r (addition g o)))
    (print-point r)
    (format t "equal: ~A~%" (equal-point g r))))

実行結果は下記の通り。

X: 5A4DB4F35B94FFB029388C7F75B7171A4B31D8254969ECB2558B5823C97547B
Y: 1999999999999999999999999999999999999999999999999999999999999999
Z: 4
  -> Valid: T
equal: T
X: 4F1970C66BED0DED221D15A622BF36DA9E146570470F1767EA6DE324A3D3A46412AE1AF72AB66511433B80E18B00938E2626A82BC70CC05E
Y: 693F46716EB6BC248876203756C9C7624BEA73736CA3984087789C1E05A0C2D73AD3FF1CE67C39C4FDBD132C4ED7C8AD9808795BF230FA14
Z: 1
  -> Valid: T
equal: T

問題なさそうですね。

6.4.1. G+Gの確認

次のような確認のコードを作りました。

(defun check-2g ()
  (let* ((g *elliptic-g*)
         (x (addition g g))
         (y (doubling g))
         (z (multiple 2 g)))
    (format t "Check: ~A, ~A, ~A, ~A, ~A, ~A~%"
            (valid g) (valid x) (valid y) (valid z)
            (equal-point2 x y) (equal-point2 x z))))

(format t "~&secp256k1: ") (with-elliptic-secp256k1 (check-2g))
(format t "~&secp256r1: ") (with-elliptic-secp256r1 (check-2g))
(format t "~&ed25519  : ") (with-elliptic-ed25519 (check-2g))
(format t "~&ed448    : ") (with-elliptic-ed448 (check-2g))

実行結果は次のとおり。

secp256k1: Check: T, T, T, T, T, T
secp256r1: Check: T, T, T, T, T, T
ed25519  : Check: T, T, T, T, T, T
ed448    : Check: T, T, T, T, T, T

6.4.2. 5Gの確認

次のような確認のコードを作りました。

(defun check-5g ()
  (let* ((g *elliptic-g*)
         (x (addition
              (addition
                (addition
                  (addition g g) g) g) g))
         (y (addition
              (addition (doubling g) (doubling g)) g))
         (z (multiple 5 g)))
    (format t "Check: ~A, ~A, ~A, ~A, ~A, ~A~%"
            (valid g) (valid x) (valid y) (valid z)
            (equal-point2 x y) (equal-point2 x z))))

(format t "~&secp256k1: ") (with-elliptic-secp256k1 (check-5g))
(format t "~&secp256r1: ") (with-elliptic-secp256r1 (check-5g))
(format t "~&ed25519  : ") (with-elliptic-ed25519 (check-5g))
(format t "~&ed448    : ") (with-elliptic-ed448 (check-5g))

実行結果は次のとおり。

secp256k1: Check: T, T, T, T, T, T
secp256r1: Check: T, T, T, T, T, T
ed25519  : Check: T, T, T, T, T, T
ed448    : Check: T, T, T, T, T, T

6.4.1. nGの確認

ここは少し詳しくやります。

まずは単純に nGを求めてみます。

(defun check-equal ()
  (let* ((g *elliptic-g*)
         (n *elliptic-n*)
         (v (multiple n g)))
    (print-point v)))

(with-elliptic-secp256k1 (check-equal))
(with-elliptic-secp256r1 (check-equal))
(with-elliptic-ed25519 (check-equal))
(with-elliptic-ed448 (check-equal))

実行結果は下記の通り。

X: 0
Y: C0763AD39DC594D1509227D46650EEC24211044CC52DB049C500EF738489D9C6
Z: 0
  -> Valid: NIL
X: 0
Y: 5E2FFD6BB8888CC31BE98BAB3E36FBB12FCE891128CBFA8E9ADF5EBFF836DE1D
Z: 0
  -> Valid: NIL
X: 0
Y: 7B28334BD670BFDB3E6C2C75B4D19ADA29D561FEB787C0ECA004B27857423FE5
Z: 7B28334BD670BFDB3E6C2C75B4D19ADA29D561FEB787C0ECA004B27857423FE5
  -> Valid: T
X: 0
Y: F75B3B69D1456D113231CD7C343FB1B95DD34B30558BDF5426A2BA922049929C833BEE409B5EC2163906BC7A835A51212C3CBB1DDD8FBF12
Z: F75B3B69D1456D113231CD7C343FB1B95DD34B30558BDF5426A2BA922049929C833BEE409B5EC2163906BC7A835A51212C3CBB1DDD8FBF12
  -> Valid: T

ちゃんと、すべてのパターンにおいて Oが返却されています。

では、この Oに Gを足してみます。

(defun check-equal ()
  (let* ((g *elliptic-g*)
         (n *elliptic-n*)
         (v (multiple n g))
         (r (addition v g)))
    (format t "equal: ~A~%" (equal-point g r))))

(with-elliptic-secp256k1 (check-equal))
(with-elliptic-secp256r1 (check-equal))
(with-elliptic-ed25519 (check-equal))
(with-elliptic-ed448 (check-equal))

実行結果は下記の通り。

equal: T
equal: T
equal: T
equal: T

問題なしです。

では最後に1000で分けたものを実行します。

(defun check-equal ()
  (let* ((g *elliptic-g*)
         (n *elliptic-n*)
         (v1 (multiple 1000 g))
         (v2 (multiple (- n 1000) g))
         (r (addition v1 v2)))
    (print-point r)))

(with-elliptic-secp256k1 (check-equal))
(with-elliptic-secp256r1 (check-equal))
(with-elliptic-ed25519 (check-equal))
(with-elliptic-ed448 (check-equal))

結果は出しませんが、ちゃんと Oでした。

これで、addition, doubling, multipleは問題なさそうです。

続きます

次は秘密鍵と公開鍵です。
Common Lispで楕円曲線DSAを実装する5（鍵生成） - nptclのブログ

Common Lispで楕円曲線DSAを実装する3（乗算など）

Common Lisp

前回の続きからです。
Common Lispで楕円曲線DSAを実装する2（加算） - nptclのブログ

5. その他の演算

次の演算をやります。

  • ä¹—ç®—
  • 逆数
    • アフィン座標への変換
  • 平方根

あとでいっぱい使います。

5.1. ä¹—ç®—

乗算を行います。
乗算は本当にいろんなところで使います。

乗算とは、座標のスカラー倍の事です。
スカラーとは単なる整数ですが、乗算の対象となるのは座標です。
よく

\displaystyle P = n G

みたいにかかれたりします。
座標 Pは、座標 Gをスカラー値である n倍した場所です。
例えば n=5の場合は、

\displaystyle P = 5 G = G + G + G + G + G

となるわけです。
この演算の難しい点としては、スカラー値 nがものすごく大きいことにあります。

RSAをやったことがある人なら、

\displaystyle P \equiv G^n \mod p

みたいな演算に頭を悩ませたことがあるかもしれません。
事情は全く同じであり、巨大なスカラー値 nをビットごとに判定して doublingとaddtionを繰り返します。

5.1.1. 乗算の実装方法

RFCのサンプルコードを見ましょう。
これはed25519のものです。

def point_mul(s, P):
    Q = (0, 1, 1, 0)  # Neutral element
    while s > 0:
        if s & 1:
            Q = point_add(Q, P)
        P = point_add(P, P)
        s >>= 1
    return Q

このコードで注意しなければいけないのは次の2点です

  • point_add(P, P)は、正しく計算できない場合がある
  • Qの初期値は、曲線の種類に合った値にする

まずpoint_add(P, P)の計算が正しく行われないときは、 doublingを別途用意してください。
すでに説明済みなので問題ないでしょう。

Qの初期値は、必ず単位元の Oにする必要があります。
つまりは次のようになります。

  • secp256k1, secp256r1のとき、 O=(0,0,0)
  • ed25519のとき、 O=(0,1,1,0)
  • ed448のとき、 O=(0,1,1)

これを実現するために、special変数を新たに用意します。

(defvar *elliptic-o*)

曲線によって適切な値を割り当ててください。
letによる例を示します。

;;  secp256k1, secp256r1
(let ((*elliptic-o* (make-point3 0 0 0)))
  ...)
;;  ed25519
(let ((*elliptic-o* (make-point4 0 1 1 0)))
  ...)
;;  ed448
(let ((*elliptic-o* (make-point4 0 1 1)))
  ...)

曲線によって初期値を変えるのが面倒な場合は Oにこだわらず、 全く別の値にする方法もあるとは思います。
例えば Qにnilを突っ込んでおいて、 point_add(Q, P)を次のようにします。

  • point_add(Q, P)
    • Qがnilの場合は Pを返却
    • そうでなければ Q+Pを返却

しかし、そのように実装した場合は、 引数のスカラーが 0のときをちゃんと想定してください。
s=0のとき、返却は Oです。
どのように実装したところで Oが必要になります。

5.1.2. 乗算の実装

Common Lispの実装を示します。

(defun multiple (s p &optional (q *elliptic-o*))
  (if (< 0 s)
    (multiple
      (ash s -1)
      (doubling p)
      (if (logbitp 0 s)
        (addition q p)
        q))
    q))

再起呼び出しの作りになっています。
最初の呼び出しでは、qが*elliptic-o*ですので、 これが Oになります。

乗算のテストは別の章でやります。

5.2. 逆数

逆数の求め方を説明します。
逆数は射影座標からアフィン座標に変換するときに用います。

例えば次の射影座標があったとします。

\displaystyle G = (X, Y, Z)

これをアフィン座標に変換するときは次の計算を行います。

\displaystyle (x, y) = (X/Z, Y/Z)

つまりは Zの逆数 Z^{-1}を求めたのち、

\displaystyle (x, y) = (X Z^{-1}, Y Z^{-1})

という乗算を計算すればいいわけです。

それでは Zの逆数を求めてみましょう。
逆数の計算は、フェルマーの小定理より次の式が導かれるとのこと。

\displaystyle Z^{-1} \equiv Z^{p-2} \mod p

ここで pは素数である必要があります。
たいていの場合、 pは曲線パラメーターの素数 pになります。

逆数といっても、有限体での出来事なので結果は整数です。
小数になるわけではありません。

まとめると、とても簡単な結論になります。
Zの逆数を求めたいなら、 Zの p-2乗を計算しろ。

5.2.2. 逆数の実装

まずは、べき乗の計算が必要です。
RSAなんかではいっつも計算しているやつです。
次の式、

\displaystyle x^y \mod n

を求める関数をpower-modという名前にします。

(defun power-mod (x y n &optional (r 1))
  (if (< 0 y)
    (power-mod
      (mod (* x x) n)
      (ash y -1)
      n
      (if (logbitp 0 y)
        (mod (* r x) n)
        r))
    r))

再起呼び出しにより実装しています。
アルゴリズムとしては乗算と同じなので、 形がとてもよく似ています。

逆数を求めるコードは非常に簡単です。
その関数をinverseとします。

(defun inverse (x)
  (power-mod x (- *elliptic-p* 2) *elliptic-p*))

5.2.3. 逆数のテスト

いくつか例を出してみます。

(let ((*random-state* (make-random-state t)))
  (with-elliptic-secp256k1
    (dotimes (i 4)
      (let* ((x (random *elliptic-p*))
             (y (inverse x))
             (z (inverse y)))
        (format t "     x: ~X~%" x)
        (format t "inv1 x: ~X~%" y)
        (format t "inv2 x: ~X~%" z)
        (format t " equal: ~A~%" (= x z))))))

実行例は下記の通り

     x: AF9064D272615B0B8F8D0207B623F34418F6B95F1929CE565625C45C311C4B10
inv1 x: 1CFE207D9EFBE7CF0B4D8826C51EB95627B1413E109B4599495CCDA19BE93D6
inv2 x: AF9064D272615B0B8F8D0207B623F34418F6B95F1929CE565625C45C311C4B10
 equal: T
     x: D559DE25AEC052C9E0A53B78317D4C966C8790408B4607E44E3F97BDB733CC87
inv1 x: 89D7E9896239C480706CD1304E2E324631808624C673AC6E21F9BB77CEDB6B61
inv2 x: D559DE25AEC052C9E0A53B78317D4C966C8790408B4607E44E3F97BDB733CC87
 equal: T
     x: B002513CEF087FC21A9F8C1CDECAEE691702A438293421AFF426D7744CD47566
inv1 x: FB9656792D4EF972D46144D1AB2489EC77D0C5CB9E7E6E3A7722401AD87A273F
inv2 x: B002513CEF087FC21A9F8C1CDECAEE691702A438293421AFF426D7744CD47566
 equal: T
     x: 735E87838D7F97172DE9DF17B81553C94BFA85235ED74D854F49D1899E32E775
inv1 x: 7DD810A7D8893A094B94E466D1417A60F2D1C8EB32A7BEF80DDE0FF0AEB00C8E
inv2 x: 735E87838D7F97172DE9DF17B81553C94BFA85235ED74D854F49D1899E32E775
 equal: T

5.2.4. アフィン座標への変換

逆数を使って、座標変換をしましょう。

ここでは射影座標からアフィン座標への変換をします。
座標 Aが次のようにあらわされるとします。

\displaystyle A = (X, Y, Z)

あるいは次のようでも問題ありません。

\displaystyle A = (X, Y, Z, T)

変換したあとのアフィン座標を (x, y)とした場合、次のような関係があります。

\displaystyle (x, y) = (X/Z, Y/Z)

つまりは、 Zの逆数 Z^{-1}を計算したのち、乗算すればいいのです。
計算は次のようになります。

\displaystyle x = X Z^{-1}
\displaystyle y = Y Z^{-1}

実装は次のようになります。

(defun affine (v)
  (let ((z (inverse (point3-z v))))
    (make-point2
      (mulp (point3-x v) z)
      (mulp (point3-y v) z))))

5.3. 平方根

平方根は、例えば座標の xだけがわかっているとき、 yを求めるときに使います。
具体的には座標のDecodeのときに使います。

では aの平方根を求めましょう。
方法は全部EdDSAのRFCに載っています。
すばらしい。

RFCの情報では、次の2通りの状況で平方根を求めることができます。

ではさっそく素数 pを割って確認します。
まずは 4で割ります。

;; secp256k1
* (with-elliptic-secp256k1 (rem *elliptic-p* 4))
3
;; secp256r1
* (with-elliptic-secp256r1 (rem *elliptic-p* 4))
3
;; ed25519
* (with-elliptic-ed25519 (rem *elliptic-p* 4))
1
;; ed448
* (with-elliptic-ed448 (rem *elliptic-p* 4))
3

次に 8で割ります。

;; secp256k1
* (with-elliptic-secp256k1 (rem *elliptic-p* 8))
7
;; secp256r1
* (with-elliptic-secp256r1 (rem *elliptic-p* 8))
7
;; ed25519
* (with-elliptic-ed25519 (rem *elliptic-p* 8))
5
;; ed448
* (with-elliptic-ed448 (rem *elliptic-p* 8))
7

以上により、次のような結果が得られました。

  • p \equiv 3 \mod 4
    • secp256k1
    • secp256r1
    • ed448
  • p \equiv 5 \mod 8
    • ed25519

それでは順番にやっていこうと思います。
次の計算を説明します。

5.3.1. p \equiv 3 \mod 4のときの平方根

計算方法はここに記載されています。

こちらはsecp256k1, secp256r1でも使えます。
RFCで見るところは、ed25519ではなくed448の方ですので注意。

あと注意してほしいのですが、 ここに書かれている方法はed448では使いません。
RFCでは、結局別の方法で実装します。
たぶん参考として載せてくれてるんだと思います。

それでは気を取り直して続きを行きます。
まずは平方根の候補として xを次の式で求めます。

\displaystyle x \equiv a^{(p+1)/4} \mod p

算出したら x^{2}を計算してください。
次の場合について。

  1. x^{2} \equiv a \mod pなら、 aの平方根は x
  2. それ以外なら aの平方根はなし

実装するとこんな感じです。

(defun square-root-mod-4 (a)
  (let* ((x (power-mod a (/ (+ *elliptic-p* 1) 4) *elliptic-p*))
         (x2 (mulp x x)))
    (if (= x2 a)
      x)))

平方根がない場合は、nilを返却します。

5.3.2. p \equiv 5 \mod 8のときの平方根

最初に言っておくと、これは使いません。
ed25519では、別の方法で算出するように説明しています。
興味がないならすっ飛ばしてください。

計算方法はここに記載されています。

まずは平方根の候補として xを次の式で求めます。

\displaystyle x \equiv a^{(p+3)/8} \mod p

算出したら x^{2}を計算してください。
その値によって、次の3つの場合があります。

  1. x^{2} \equiv a \mod pなら、 aの平方根は x
  2. x^{2} \equiv -a \mod pなら、 aの平方根は 2^{(p-1)/4} \cdot x
  3. それ以外なら aの平方根はなし

なお -a \mod pは、 p-a \mod pと同じです。

実装するとこんな感じです。

(defun square-root-mod-8 (a)
  (let* ((x (power-mod a (/ (+ *elliptic-p* 3) 8) *elliptic-p*))
         (x2 (mulp x x)))
    (cond ((= x2 a) x)
          ((= x2 (- *elliptic-p* a))
           (mulp x (power-mod 2 (/ (- *elliptic-p* 1) 4) *elliptic-p*))))))

平方根がない場合は、nilを返却します。

5.3.3. 実装のテスト

乱数で適当な値 xを求め、 x^{2}を求めてから 実装したコードで元の xを求めてみます。

まずはsquare-root-mod-4から。

(let ((*random-state* (make-random-state t)))
  (with-elliptic-secp256k1
    (let* ((x (random *elliptic-p*))
           (y (modp (- *elliptic-p* x)))
           (a (mulp x x))
           (z (square-root-mod-4 a)))
      (format t "~X~%" x)
      (format t "~X~%" y)
      (format t "~X~%" z)
      (format t "equal: ~A~%" (or (= x z) (= y z))))))

乱数の値 xがマイナスであることも考慮してください。
マイナスの時は、二乗した結果プラスになるので、 その平方根では最初の値と一致しません。
上記の判定では、正と負の両方の場合を考慮しています。

正の場合の実行例を示します。

3475A6298CF6C10E3586EE8A12BCD0F584F1579AFE4D53351DE19DE9D5300572
CB8A59D673093EF1CA791175ED432F0A7B0EA86501B2ACCAE21E62152ACFF6BD
3475A6298CF6C10E3586EE8A12BCD0F584F1579AFE4D53351DE19DE9D5300572
equal: T

負の場合の実行例を示します。

617559E463B7B0B19743D53F53546966723BFBBD96BB07544C21102307F3B116
9E8AA61B9C484F4E68BC2AC0ACAB96998DC404426944F8ABB3DEEFDBF80C4B19
9E8AA61B9C484F4E68BC2AC0ACAB96998DC404426944F8ABB3DEEFDBF80C4B19
equal: T

square-root-mod-8も同じです。

(let ((*random-state* (make-random-state t)))
  (with-elliptic-ed25519
    (let* ((x (random *elliptic-p*))
           (y (modp (- *elliptic-p* x)))
           (a (mulp x x))
           (z (square-root-mod-8 a)))
      (format t "~X~%" x)
      (format t "~X~%" y)
      (format t "~X~%" z)
      (format t "equal: ~A~%" (or (= x z) (= y z))))))

正の場合の実行例を示します。

73EC17B0C2C37B013F0793E1C0BDC116118231F773262F18FC045981F57BFF2A
C13E84F3D3C84FEC0F86C1E3F423EE9EE7DCE088CD9D0E703FBA67E0A8400C3
73EC17B0C2C37B013F0793E1C0BDC116118231F773262F18FC045981F57BFF2A
equal: T

負の場合の実行例を示します。

518EB5A3821819110B19C0925EC335DFEDE8F70F8040A8715B83AFA42D8EC2B3
2E714A5C7DE7E6EEF4E63F6DA13CCA20121708F07FBF578EA47C505BD2713D3A
2E714A5C7DE7E6EEF4E63F6DA13CCA20121708F07FBF578EA47C505BD2713D3A
equal: T

5.3.4. 符号について

今回座標を求めるために平方根を算出するのですが、 符号だけは正しく算出できません。
そこで符号の情報だけ別で保存しておき、あとで復帰するという方法が行われます。

例えばアフィン座標 (x,y)の xだけを保存するとします。
その場合でも、 yの符号だけは別途保存してください。
この場合、符号情報は yの最下位ビットです（最上位ではない）。

あらかじめ保存しておいた符号情報を y_0とします。
y_0には 0か 1が入っています。
いま、何らかの値 aの平方根を用いて、仮決めの yを得たとします。

(setq y (square-root-mod-4 a))

この状態では、 yの符号があっていない可能性があります。
そこで、仮決めの yの最下位ビットと y_0が等しいか確認します。
等しくない場合は符号を反転させます。

(when (/= (logand y #x01) y0)
  (setq y (- *elliptic-p* y)))

これで、符号 y_0を含めて、正しい yを求めることができました。

5.3.5. 平方根の求め方について

「平方剰余」というらしいです。
もしまじめにやりたいなら平方剰余で調べてください。

求め方は次のサイトが参考になります。

平方根の出し方はsecg.orgのPDFには載ってませんでした。
探し方が悪かっただけであるのかもしれません。
偶然EdDSAのRFCに載ってるやつで代用できたのでよかったです。

続きます

次はここまでのテストです。
Common Lispで楕円曲線DSAを実装する4（確認） - nptclのブログ

Common Lispで楕円曲線DSAを実装する2（加算）

Common Lisp

前回の続きからです。
Common Lispで楕円曲線DSAを実装する1 - nptclのブログ

4. 加算

楕円曲線暗号の加算を行います。
なんでいきなり加算なんだと不思議に思う方には申し訳ありません。
自分では説明できません。
わかんなくても、やっていけば後で使います。

点 Pと点 Qの加算 P+Qを次のように表します。

\displaystyle P = (x_1, y_1)
\displaystyle Q = (x_2, y_2)
\displaystyle P + Q = (x_3, y_3)

つまりこれを計算します。

\displaystyle (x_3, y_3) = (x_1, y_1) + (x_2, y_2)

単なるベクトルの演算ではないので注意してください。
加算の演算内容は使う曲線によって変わります。
上記4つはそれぞれ違った形の曲線ですのでひとつずつ見ていきましょう。

一応言っておきますが、全部整数で計算します。
浮動小数は出てきません。

4.1. 加算の計算方法

便利な世の中になったもので、 ed25519とed448については、 RFCに詳しく記載されていました。
もう何も考えずに丸ごと見てしまいましょう。

RFC8032
Edwards-Curve Digital Signature Algorithm (EdDSA)
https://datatracker.ietf.org/doc/html/rfc8032

この「5.1.4. Point Addition」に 加算のすべてが書かれています。
式を使うことがないといったのは、 すでにRFCに手順が載っているからです。

ただ覚えておかなければいけないことがあります。
5.1.4.を見るといきなり座標変換の話から始まっています。
「extended homogeneous coordinates」だそうです。
検索すると「拡張された同次座標系」と出てきますが、 他の人たちのサイトを見る限りだと「同次座標系」ではなく 「射影座標 (projective coordinates)」と呼んでいました。
おそらくは「拡張された射影座標」でよいのかと思います。
これは何者かというと、今まで座標に xと yを用いていましたが、 それを (X, Y, Z, T)の4つにしましょうという話です。
単なる射影座標だと (X, Y, Z)の3つだったと思うので、 拡張されて4つになったんでしょうね。

今までの (x, y)はアフィン座標というらしいです。
対応関係は下記の通り。

  • x = X/Z
  • y = Y/Z
  • x \cdot y = T/Z

Z=1のときアフィン座標そのものであり、 Z=0のときは無限遠点 Oになるんだそうです。
アフィン座標 (x, y)から射影座標へは (x, y, 1, x \cdot y)になります。
逆に (X, Y, Z, T)をアフィン座標にするときは、 Zが 0ではないとき (X/Z, Y/Z)です。

何のためにこんなことするのかというと、演算の最適化のためです。
めちゃ時間がかかる乗算や除算の処理を減らして、 可能な限り加算に回そうという工夫らしいです。
でも今はそんなことより、RFCで公開されている手順に従って 楽に加算を実現させたいという思いしかありません。

加算は次のようになると書いています。

A = (Y1-X1)*(Y2-X2)
B = (Y1+X1)*(Y2+X2)
C = T1*2*d*T2
D = Z1*2*Z2
E = B-A
F = D-C
G = D+C
H = B+A
X3 = E*F
Y3 = G*H
T3 = E*H
Z3 = F*G

こちらはed25519の加算です。
こんなのをもこもこ計算するだけで、

\displaystyle (X_3, Y_3, Z_3, T_3) = (X_1, Y_1, Z_1, T_1) + (X_2, Y_2, Z_2, T_2)

が得られるということになります。
面白いですね。
今説明したのが「adding」または「addition」と呼ばれるものです。
そのほかにも点 (x, y)を2倍する「doubling」も紹介されています。
doublingはスカラー乗算で使います。

5.2.4.ではed448の加算が紹介されています。
ed448や、secp256k1, secp256r1は、 普通の射影座標 (X, Y, Z)であることに注意してください。

さらに6. Ed25519 Python Illustrationでは、 ed25519の実装例を見ることができます。
もはや欲しいものは全部手に入りました。
加算以外にも考えなければいけないことはあるのですが、 このRFCに全部載っています。
本当に素晴らしい。

ただし、EdDSAだけです。
自分が調べた限りだと、ECDSAについてはここまで詳しく書かれた RFCはありませんでした。
本当にないのかな。
ECDSAバージョンは自分で作るしかなさそうです。
ありがたいことに手順を公開しているサイトがあります。

Explicit-Formulas Database
https://hyperelliptic.org/EFD

このページの「Twisted Edwards curves: a*x2+y2=1+d*x2*y2」が、 いま見てきた情報です。
下記のページに全く同じ計算が載っています。
https://hyperelliptic.org/EFD/g1p/auto-twisted-extended-1.html

ECDSAの場合は「Short Weierstrass curves: y2=x3+a*x+b」を選んでください。
その中に場合分けされた「Projective coordinates」が いくつも紹介されているので、 それに従って作成できるはずです。
あとで説明します。

まずは情報があるものから作りましょう。
ed25519から行きます。

4.2. 加算の実装

まずはサンプルを見ましょう。
Pythonによる実装が次のようになります。

def point_add(P, Q):
    A, B = (P[1]-P[0]) * (Q[1]-Q[0]) % p, (P[1]+P[0]) * (Q[1]+Q[0]) % p;
    C, D = 2 * P[3] * Q[3] * d % p, 2 * P[2] * Q[2] % p;
    E, F, G, H = B-A, D-C, D+C, B+A;
    return (E*F, G*H, F*G, E*H);

これを一列に並び変えます。

def point_add(P, Q):
    A = (P[1]-P[0]) * (Q[1]-Q[0]) % p;
    B = (P[1]+P[0]) * (Q[1]+Q[0]) % p;
    C = 2 * P[3] * Q[3] * d % p;
    D = 2 * P[2] * Q[2] % p;
    E = B-A;
    F = D-C;
    G = D+C;
    H = B+A;
    return (E*F, G*H, F*G, E*H);

あとはこれをC言語やらに地道に変換していけばいいわけです。

これらの演算は、有限体と呼ばれる状況下で行われます。
Pythonは超高級な言語なので、 掛け算とかを連続で行っても適切にメモリ確保してくれますが、 ちょうどいいタイミングになったら 素数 pで割り算してあまりにしてやりましょう。
C言語だと、どのタイミングであまりを求めればいいのか、 結構真剣に考えることになります。

以降では、使用する演算を列挙していこうと思います。
結局は全部下記のページのものです。

Explicit-Formulas Database
https://www.hyperelliptic.org/EFD/index.html

列挙する内容は次の通り。

  • secp256k1
  • secp256r1
  • ed25519
  • ed448

それぞれに対しての次の演算方法です。

  • addition
  • doubling

dounlingは、additionで代用できる可能性があります。
最後に実装の注意点をまとめます。

あと、座標系に注意してください。
射影座標 (X, Y, Z)は次のものが該当します。

  • secp256k1
  • secp256r1
  • ed448

拡張された射影座標 (X, Y, Z, T)は次のものが該当します。

  • ed25519

4.3. addition

自分が使用した実装を列挙します。
最後に注意点をまとめますのでご確認ください。

4.3.1. secp256k1, secp256r1

使用する座標は、射影座標 (X, Y, Z)です。

Projective coordinates for short Weierstrass curves
https://hyperelliptic.org/EFD/g1p/auto-shortw-projective.html#addition-add-2007-bl

U1 = X1*Z2
U2 = X2*Z1
S1 = Y1*Z2
S2 = Y2*Z1
ZZ = Z1*Z2
T = U1+U2
TT = T^2
M = S1+S2
R = TT-U1*U2+a*ZZ^2
F = ZZ*M
L = M*F
LL = L^2
G = (T+L)^2-TT-LL
W = 2*R^2-G
X3 = 2*F*W
Y3 = R*(G-2*W)-2*LL
Z3 = 4*F*F^2

4.3.2. ed25519

使用する座標は、拡張された射影座標 (X, Y, Z, T)です。

Extended coordinates with a=-1 for twisted Edwards curves
https://www.hyperelliptic.org/EFD/g1p/auto-twisted-extended-1.html#addition-add-2008-hwcd-3

RFC8032
https://datatracker.ietf.org/doc/html/rfc8032

A = (Y1-X1)*(Y2-X2)
B = (Y1+X1)*(Y2+X2)
C = T1*2*d*T2
D = Z1*2*Z2
E = B-A
F = D-C
G = D+C
H = B+A
X3 = E*F
Y3 = G*H
T3 = E*H
Z3 = F*G

T3とZ3の順番に注意。

4.3.3. ed448

使用する座標は、射影座標 (X, Y, Z)です。

Projective coordinates for Edwards curves
https://www.hyperelliptic.org/EFD/g1p/auto-edwards-projective.html#addition-add-2007-bl

RFC8032
https://datatracker.ietf.org/doc/html/rfc8032

A = Z1*Z2
B = A^2
C = X1*X2
D = Y1*Y2
E = d*C*D
F = B-E
G = B+E
H = (X1+Y1)*(X2+Y2)
X3 = A*F*(H-C-D)
Y3 = A*G*(D-C)
Z3 = F*G

4.4. doubling

自分が使用した実装を列挙します。
最後に注意点をまとめますのでご確認ください。

4.4.1. secp256k1, secp256r1

Projective coordinates for short Weierstrass curves
https://hyperelliptic.org/EFD/g1p/auto-shortw-projective.html#doubling-dbl-2007-bl

XX = X1^2
ZZ = Z1^2
w = a*ZZ+3*XX
s = 2*Y1*Z1
ss = s^2
sss = s*ss
R = Y1*s
RR = R^2
B = (X1+R)^2-XX-RR
h = w^2-2*B
X3 = h*s
Y3 = w*(B-h)-2*RR
Z3 = sss

4.4.2. ed25519

Extended coordinates with a=-1 for twisted Edwards curves
https://www.hyperelliptic.org/EFD/g1p/auto-twisted-extended-1.html#doubling-dbl-2008-hwcd

RFC8032
https://datatracker.ietf.org/doc/html/rfc8032

A = X1^2
B = Y1^2
C = 2*Z1^2
H = A+B
E = H-(X1+Y1)^2
G = A-B
F = C+G
X3 = E*F
Y3 = G*H
T3 = E*H
Z3 = F*G

T3とZ3の順番に注意。

4.4.3. ed448

Projective coordinates for Edwards curves
https://www.hyperelliptic.org/EFD/g1p/auto-edwards-projective.html#doubling-dbl-2007-bl

RFC8032
https://datatracker.ietf.org/doc/html/rfc8032

B = (X1+Y1)^2
C = X1^2
D = Y1^2
E = C+D
H = Z1^2
J = E-2*H
X3 = (B-E)*J
Y3 = E*(C-D)
Z3 = E*J

4.4. 実装の注意点

いろいろと注意しなければいけないことがあります。
もし自分で実装するなら、この章の内容を覚えておいた方がいいと思います。

話題は次の通り。

  • doublingã‚’additionで代用するときの注意
  • 無限遠点 Oの加算
  • 加算の実装例

4.4.1. doublingをadditionで代用するときの注意

doublingはaddtionで代用できるかもしれません。
doublingとはつまり、

\displaystyle P = Q + Q

の事なので、

P = doubling(Q) = addtion(Q, Q)

で計算できる場合があります。
なんで断言しないかというと、できない場合もあるからです。
RFCに載っている、ed25519とed448は代用できます。
わざわざdoublingも載せてくれていますが、 面倒ならadditionだけで大丈夫です。

問題はsecp256k1とsecp256r1なのです。
自分は、下記のサイトから好きなのを選んで実装していました。

Projective coordinates for short Weierstrass curves
https://hyperelliptic.org/EFD/g1p/auto-shortw-projective.html

ところが、下記のもの

https://hyperelliptic.org/EFD/g1p/auto-shortw-projective.html#addition-add-1998-cmo-2

は、doublingの代用ができませんでした。
つまり次の計算はできません。

addition(P, P)

結果は X, Y, Zが全部 0になります。

元の論文を見たわけではないので、それが良いのか悪いのか全然わかりません。
このような場合は、加算をするまえに二つの座標を確認する必要がありそうです。
もし問題ありそうな場合、doublingに切り替えるような実装になると思います。
自分には難しすぎるので、違う式で算出することにしました。

additionのコードはたくさんありますが、 実際に使ってみない限りわからないこともあると思います。
自分で実装する場合は、いろいろ作って試してみてください。

4.4.2. 無限遠点 Oの加算

無限遠点 Oの加算は、ちゃんと考えなければいけません。
そして、できるようにしなければなりません。

まずは、いったい何が無限遠点 Oなのかをまとめます。
おそらくは次のような判定になると思います。

  • secp256k1, secp256r1
    • Z=0なら O
  • ed25519
    • アフィン座標 (0,1)なら O、例えば射影座標 (0, 1, 1, 0)
  • ed448
    • アフィン座標 (0,1)なら O、例えば射影座標 (0, 1, 1)

ed25519とed448は具体的な座標があり、 これを中立元(neutral point)というらしいです。
もしかしたら中立元と Oは同一ではないのかもしれませんが、 完全に Oとして機能しています。
この辺、勉強不足でよくわかりません。

次の2つの性質を持っています。

  • P = O+P = P+O
  • nG = O

ということなので、ed25519とed448なら無限遠点 Oの加算は、 中立元がうまいことやってくれるので問題ありません。
しかしsecp256k1とsecp256r1はそうではありません。
ちゃんと考えないとできません。
普通にそのまま加算すると、変な値になります。

secp256k1とsecp256r1の場合は次のように場合分けをします。

  • 次の手順で P + Qを実行する
    • Pの Zが 0なら、 Qを返却
    • Qの Zが 0なら、 Pを返却
    • それ以外なら、additionのコードで計算する

上記のように Zによる場合分けの実装は、 どうも必ずやらなければいけないようです。
次に実装の例を示します。

4.4.3. 加算の実装例

場合分けが必要な、secp256k1のadditionについて示します。

(defun addition-weierstrass (p1 p2)
  (let* ((x1 (point3-x p1))
         (y1 (point3-y p1))
         (z1 (point3-z p1))
         (x2 (point3-x p2))
         (y2 (point3-y p2))
         (z2 (point3-z p2))
         (u1 (* x1 z2))
         (u2 (* x2 z1))
         (s1 (* y1 z2))
         (s2 (* y2 z1))
         (zz (* z1 z2))
         (t1 (+ u1 u2))
         (t2 (* t1 t1))
         (m (+ s1 s2))
         (r (+ (- t2 (* u1 u2)) (* *elliptic-a* zz zz)))
         (f (* zz m))
         (k1 (* m f))
         (k2 (* k1 k1))
         (g1 (+ t1 k1))
         (g2 (- (* g1 g1) t2 k2))
         (w (- (* 2 r r) g2)))
    (make-point3
      (modp (* 2 f w))
      (modp (- (* r (- g2 (* 2 w))) (* 2 k2)))
      (modp (* 4 f f f)))))

(defun addition-secp256k1 (p1 p2)
  (cond ((zerop (point3-z p1)) p2)
        ((zerop (point3-z p2)) p1)
        (t (addition-weierstrass p1 p2))))

addition-weierstrassだけだと、 無限遠点 Oの加算の時に失敗します。
そこで Zの値を見て場合分けをすることで対処します。
これをしないと、加算や乗算の最中に偶然にも Oにあたってしまった場合、 正しく算出できません。
具体的な症状としては、後々の話になりますが、公開鍵を生成するときに失敗します。

ひとつだと寂しいので、ed25519のadditionを示します。

(defun addition-ed25519 (p1 p2)
  (let* ((x1 (point4-x p1))
         (y1 (point4-y p1))
         (z1 (point4-z p1))
         (t1 (point4-xy p1))
         (x2 (point4-x p2))
         (y2 (point4-y p2))
         (z2 (point4-z p2))
         (t2 (point4-xy p2))
         (a (* (- y1 x1) (- y2 x2)))
         (b (* (+ y1 x1) (+ y2 x2)))
         (c (* t1 2 *elliptic-d* t2))
         (d (* z1 2 z2))
         (e (- b a))
         (f (- d c))
         (g (+ d c))
         (h (+ b a)))
    (make-point4
      (modp (* e f))
      (modp (* g h))
      (modp (* f g))
      (modp (* e h)))))

続きます

次は乗算やらいろんな演算をします。
Common Lispで楕円曲線DSAを実装する3（乗算など） - nptclのブログ