前回の続きからです。
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$ が、 $p \equiv 3 \mod 4$ のもの
素数 $p$ が、 $p \equiv 5 \mod 8$ のもの

ではさっそく素数 $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

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

- secp256k1
- secp256r1
- ed448
- ed25519

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

$p \equiv 3 \mod 4$ のときの平方根
$p \equiv 5 \mod 8$ のときの平方根

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

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

5.2. Ed448ph and Ed448
5.2.1. Modular Arithmetic
https://datatracker.ietf.org/doc/html/rfc8032

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

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

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

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

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

$x^{2} \equiv a \mod p$ なら、 $a$ の平方根は $x$
それ以外なら $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では、別の方法で算出するように説明しています。
興味がないならすっ飛ばしてください。

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

5.1. Ed25519ph, Ed25519ctx, and Ed25519
5.1.1. Modular Arithmetic
https://datatracker.ietf.org/doc/html/rfc8032

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

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

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

$x^{2} \equiv a \mod p$ なら、 $a$ の平方根は $x$
$x^{2} \equiv -a \mod p$ なら、 $a$ の平方根は $2^{(p-1)/4} \cdot x$
それ以外なら $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. 平方根の求め方について

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

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

biginteger（任意精度整数）でelliptic curve/secp256k1（楕円曲線）
平方根（square root）
https://qiita.com/tnakagawa/items/c0396c294eaf67f55fdc
平方剰余
素数または素数べき乗の法
https://ja.wikipedia.org/wiki/%E5%B9%B3%E6%96%B9%E5%89%B0%E4%BD%99
Quadratic residue
Prime or prime power modulus
https://en.wikipedia.org/wiki/Quadratic_residue

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

続きます

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

nptclのブログ

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

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

5. その他の演算

5.1. 乗算

5.1.1. 乗算の実装方法

5.1.2. 乗算の実装

5.2. 逆数

5.2.2. 逆数の実装

5.2.3. 逆数のテスト

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

5.3. 平方根

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

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

5.3.3. 実装のテスト

5.3.4. 符号について

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

続きます

5. その他の演算

5.1. 乗算

5.1.1. 乗算の実装方法

5.1.2. 乗算の実装

5.2. 逆数

5.2.2. 逆数の実装

5.2.3. 逆数のテスト

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

5.3. 平方根

5.3.1. のときの平方根

5.3.2. のときの平方根

5.3.3. 実装のテスト

5.3.4. 符号について

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

続きます

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

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