クラスを用いて有限体上の楕円曲線を構成してみる

はじめに

　このブログはmotercontrolman様主催のMATLAB/Simulink Advent Calendar 2019の15日目に寄稿した記事です。普段は生体信号解析についての研究でMATLABを日々使っていたのですが、来年度からゲームプログラマとして某社に勤務することになったので、MATLABともこれでお別れという事になり、折角なので今回のAdvent Calenderに参加する事にしました。

　今回の記事では楕円曲線という非常に面白い性質を持った代数曲線をMATLABに実装されているクラスを使って実装してみました。楕円曲線は暗号分野に用いられる（有名な例だとビットコインで楕円曲線暗号が使われています）他、数論の分野でもフェルマーの最終定理に使われていたりとかもしています。記事に書く内容を色々と試行錯誤していた関係もあって解説が足らない所がありますので、この記事を見る方は以下の3点を予め把握していると理解が早いかなと思います。

　当初は楕円曲線暗号およびCM法と呼ばれる手法の実装まで手を出す予定でしたが、楕円曲線暗号（および楕円関数論）の想像を遥かに超える沼の深さにそのままズボズボと入っていき、気が付いたら後数日前というところまで十分に理解できずに終わってしまいました……（遅刻は本当にごめんなさい！）。なので、妥協案としてこの記事では恐らく被りがないと思われる、MATLABを用いたクラス設計の紹介をメインとし、おまけとして楕円曲線について勉強してなんとか感覚的に理解できたところをおまけとして紹介することにしました。

　という訳で以下から本題に入っていきます！　プログラムも数学も初心者ですが何か興味深い点が見つかれば幸いです！

目次

楕円曲線の簡単な紹介

　最初に実装するものの簡単な紹介をしたいと思います。

　楕円曲線というのは

　　　　 y^{2} = x^{3} + ax + b

より表される代数曲線のことを言い、座標 (x, y)がどちらも有理数である点 P, Q \in \mathbb{Q} \times \mathbb{Q}について以下の加法則 R = P+Qが成り立ちます。

P=(P_x, P_y), Q=(Q_x, Q_y), R=(R_x, R_y)とすると、

[tex: R_x = grad2 - P_x - Q_x] R_y = grad*(P_x-R_x)-P_y

となる。ただし、

[P \neq Qの場合] 　　　　 grad = \frac{R_y-P_y}{R_x-P_x}

[ P=Qの場合] 　　　　 grad = \frac{3{P_x}^2+a}{2P_x}

　更に、有限体 \mathbb{F}_p上の楕円曲線というのは、楕円曲線のパラメータ a, b、楕円曲線上の点の座標 x, yが全て pで割ったあまりで表現される数で表される楕円曲線です。 　

有限体上の楕円曲線の実装

　という訳で早速、有限体上の楕円曲線の実装を紹介したいと思います。今回は群構造上の加法などの複雑な処理を"上手く"隠蔽するために、楕円曲線のクラスとその楕円曲線上の点を表すクラスの2つを書きました。後者の"点"を表すクラスが前者のクラスのインスタンスを、自分の所属する楕円曲線として参照して、加法を適用する時に適切な場所を移動するといったロジックで動いています。

楕円曲線のクラス

gist0ea559ee06d5f55abdf12e2f69b1f9ee

コンストラクタ

　とりあえず、有限体上の楕円曲線が

　　　　 y^{2} = x^{3} + ax + b \; (\text{mod}\; p)

より表されるので、それに対応するパラメータとして a, b, pをプロパティとして入れました。

あと、これ以降に紹介するメソッドもそうなのですが、64bit以上の巨大数の計算はSymbolic Math Toolboxを利用しています。最近知ったツールボックスですが、特に数式とかの記述に凄く便利で独自のsym型として数値を書いてやると、 \sqrt{2}や \frac{1}{2}などをそのままの形で表示してくれたり、導関数・原始関数の計算をしてくれたり、一応の数学好きとしては色々と感動しまくる機能ばかりでした！　なので、コンピュータ上での高度な数式の表現など興味がある方は使ってみると色々と道が拓けるかも？

y^{2}の計算

　たまにエラーチェックなどを目的として、楕円曲線の式の左辺を計算する場合があるので、それを計算しておきます。

　

楕円曲線上の点のクラス

gista44b4c68e866e17e2ffa9c9fe751d9a3

無限遠点の値

　最初に無限遠点の値がINFINITE_POINT=-9999とConstantを付けた上で定義されていますが、これは実はお気持ち程度のもので、点が無限遠点であるかどうかは条件判別的に判断しています。

コンストラクタEp(E, x, y)

　それぞれ、

に対応しています。また、中を見てみれば分かりますが、Ep(E), Ep(E, x), Ep(E, x, y)の3種類の入力に対応しており、それによって振る舞いが大きく異なります。

[Ep(E)の場合]

　自動的に無限遠点となります。具体的にはisInfinitePointというパラメータをtrueにします。

[Ep(E, x)の場合]

　楕円曲線Eの座標 xに対応する2つの座標の内、 yが正のものを設定します。また、入力されたx座標から楕円曲線上の有理点が作れない場合はエラーを返します。

[Ep(E, x, y)の場合]

　そのまま出力しますが、もしも楕円曲線上の点として認められなかった場合はエラーを返します。

群構造上の加算

　"+"演算子をオーバーロードする為にメソッド名を『plus』にしています。これによって、資料などで見る形そのままで群構造上の足し算を実行できます。主な仕様としては――

  • 一方が無限遠点だった場合は、もう片方の点をそのまま出力
  • 入力された2つの点の y座標が符号だけ逆の関係になっていた場合は、無限遠点を出力

　となり、そのどちらにも当てはまらなかったら、通常の群構造上の足し算を実行します。

プロパティ

  • elipticCurve

　楕円曲線の参照を持っている変数です。（あるいはポインタ変数のようなもの）

  • isInfinitePoint

　この点が無限遠点ならtrueになります。

計算結果

　試しに、

　　　　 y^{2} = x^{3} + x + 1 \; (\text{mod}\; 13)

　をベースにして色々とやってみました。なので最初に

>> E = ElipticCurve(1, 1, 13);
>> [E.a E.b E.p]
 
ans =
 
[ 1, 1, 13]

　より楕円曲線を作りました。シンボリック数だとプロパティとして値を表示した時にすぐに値が見れないのでこのように個別に値を呼び出してやる必要があります。

　次に無限遠点 Oを作ってみます。

>> O = Ep(E)

O = 

  Ep のプロパティ:

                  x: -9.99900000000000e+003
                  y: -9.99900000000000e+003
       elipticCurve: [1×1 ElipticCurve]
    isInfinitePoint: 1

　これも良さそうですね！　続けて、 P = (0, 1)と Q = (12, -5)を2通りのやり方で作ってみます。

>> P = Ep(E, 0);
>> [P.x, P.y]
 
ans =
 
[ 0, 1]
 

>> Q = Ep(E, 12, -5);
>> [Q.x, Q.y]
 
ans =
 
[ 12, -5]

　無事に正しい値を作れました！ 　後は、この3点を使って群構造上の加算が上手くいくかを調べたいと思います！

>> R = P+O;
>> [R.x, R.y]
 
ans =
 
[ 0, 1]

　単位元との計算は上手く行きました。しかし、加法となると――

>> P+P
エラー: Ep (line 46)
入力された座標は楕円曲線上の点ではありません

エラー:  +  (line 103)
            R = Ep(P.elipticCurve, x3, y3);

　とエラーが出てきてしまいます。どうやらP+Pの結果で出てきた座標が（少なくとも）有理点と見なされなかったようです。そこで、エラーが起きる前にブレイクポイントを掛けて、計算された座標の値を調べてみると――

>> P+P
45                  if mod(sym(y)^2, E.p) ~= E.CalculateY2(x)
K>> [x y]
 
ans =
 
[ 1/4, 95/8]

パッと見で合ってないように思えますが、 mattyuu.hatenadiary.com

mattyuuの数学ネタ集様のこの記事によると、有限体にも除法は定義ができるので、それに則って整数に直して本当にあってないかどうか調べてみます。すると、

x = 1/4=10 \in \mathbb{F}_p

y = 95/8=7 \in \mathbb{F}_p

となるので、以下のように有理点になるかどうか確かめると

>> E.CalculateY2(10)
 
ans =
 
10

>> mod(7^2, 13)

ans =

    10.0000000000000e+000

　と、これも実際には楕円曲線上の有理点となることが確認されました！　ただ、分数の形のままで同じ様な計算をすると、異なる値となってしまうため、加法のやり方をもうちょっと工夫するか、エラーチェックのやり方を工夫するかして直したいところなのですが、現状中の人がまだシンボリック数の取り扱い方が分からないのでこの辺りについては今後の課題としたいと思います……。

最後に

　最後はちょっとうまく行きませんでしたが、とりあえずクラスを利用して楕円曲線を比較的扱いやすい形で取り扱えるという点については理解できたのかな？　と個人的には感じております。

　あまりの沼にハマってしまったせいで夢中で勉強しすぎた＋計画性なさすぎたせいで申し訳ないことに結構な遅刻をかましてしまいましたが、楕円曲線および楕円曲線暗号に関わる理論は勉強してて凄く楽しく、実装に関しても自分の実力を知る機会になれたのでとても良い経験ができました。

　最後に今回のアドベントカレンダーの企画を立ち上げてくれたmotercontrolman様はじめ、この記事を最後までご覧になってくれた方、本当にありがとうございました！

　もしよろしければ、この後に続くおまけに（つたない文章ではありますが）楕円曲線の由来について色々とまとめたものがありますので、更に時間のある方はご覧になってください！　あと、githubの方に楕円曲線の構築法の一つとして知られるCM法に手を出した痕跡があるので、物好きな人はそちらもご覧になってください……

github.com

楕円曲線について（おまけ）

　楕円曲線を簡単に説明すると、以下に記述されるワイエルシュトラスの標準形（を変形したもの）

　　　　 y^{2} = x^{3} + ax + b

によって記述される代数曲線です。一見すると、そこまで複雑な式には見えませんが、これが楕円曲線暗号を強固な方式へと足らしめる基盤となっております。ただ、この式が強い暗号方式を導く理由については、ガウス、ヤコビ、アーベルを初めとした名だたる数学者たちが築き上げた楕円関数論という分野から考えるしかありません。

　楕円関数論についてもこの記事を書くための調査でやっと触れた程度で、まだあまり理解しているとは言えないのですが、ここでは楕円関数論の初歩的な部分を軽く触れつつ、ワイエルシュトラスの標準形が強固な暗号方式を築き上げる理由について書いていこうと思います。

楕円関数の発見まで

　知っている方も多いかと思われますが、工学系の方々ならば何度もお世話になっているであろう三角関数――つまり、 \sin, \cos, \tanの3つは円の弧の長さを求める関数の逆関数からも得られます。

　つまり、円の弧の長さを求める関数を弧長積分と呼び、これは微分幾何で曲線の長さを求める為の公式を使って

l(u) = \int_{0}^{u} {\frac{1}{\sqrt{1-x^{2}}}dx}

と表され、この l(u)の逆関数というのが u = \sin{l}ということです。この時、 \cos{l}も円の弧長積分の導関数の逆数から定義することができます。

　では、弧長積分の対象となるのが円ではなく楕円だった場合はどうなるのか？　これについては答えのみ言いますが、確かに円に対しての弧長積分と同様に sinのような周期関数は得られました。しかし、その周期の性質というのがそれらと次元が違っていて、円の場合だと1次元の数直線上の端っこでループするものが得られたのに、楕円の場合だと複素平面上に敷き詰められた平行四辺形のタイルの境界でループするようになったのです！　そして、このような"2重の周期性"を持つ関数（2重周期関数）は、楕円の弧長積分の逆関数から求められることに由来している為、楕円関数と呼ばれるようになりました。

　ここで楕円関数が2重周期関数であるのを発見するまでの過程まで書きたかったのですが、時間的な関係によって断念しました……（一言で言うならば、楕円関数についての3倍角,5倍角の公式を因数分解すると虚数根が見つかった事から始まっております）。もし、気になる方がいれば、繭野孝和著の『楕円関数論への入門』にこの経緯も分かりやすく書かれていたのでこちらをご覧になってください。

楕円関数上の点の振る舞い

　ここで楕円関数についての顕著な振る舞いとして"2重周期性"が出てきたので、楕円曲線が強固な振る舞いをする根拠を記してみます。

　まずは、代表的な暗号方式であるRSA暗号で基本となる

　　（12/16：式に間違いが合ったので修正しました）

　　 de=1 \text{mod} (p-1)(q-1)

　についてですが（ただし、 d, e, p, qは全て整数で、 p, qは互いｎ異なる素数とします）、これについては点 eが1次元の直線上をループしながら一定の歩数 d分歩いていると言ったイメージです。一方で、楕円関数は2重周期性より、2次元の平面上をループしながら一定方向・同じ歩数分歩いていると解釈しなければなりません。

　RSAも楕円曲線もどちらも「ある位置からループする世界を何回か歩いてまた別の位置に付いた後に何歩歩いたか忘れた場合、その歩数を計算で求めるのは難しい」（離散対数問題の例えです）という点を暗号の安全性の根拠としていますが、やはり"ループする世界の広さ"などによってその難易度は変化します。そして、楕円曲線の場合はそのループする世界の広さが2次元で与えられますので、1次元のRSAと比べると格段に離散対数問題の解を得るのが難しそうなのが直感的に分かると思います！（そして、これが楕円曲線暗号でよく言われる「RSAと比べて格段に小さい鍵長で同じ強度のセキュリティが築ける」点に由来しています）

ワイエルシュトラス標準形１：楕円積分の表示を変える

　楕円曲線の強さを説明したところで、暗号でお馴染みのワイエルシュトラス標準形がここからどのようにして導かれるかを説明したいと思います。

　まず、一般的な話から始めると楕円に対する弧長積分 Fは4次関数 Q(\omega)を用いて

　　　　 F(\omega) = \int_{0}^{\omega} {\frac{1}{\sqrt{Q(\omega)}}d\omega}

と書くことができます。ここで、4次関数 Q(\omega)を

　　　　 Q(\omega) = \alpha_{0}({\omega - \alpha_1})({\omega - \alpha_2})({\omega - \alpha_3})({\omega - \alpha_4})

4つの根に対して因数分解した形で表し、その上で楕円積分の形を分かりやすく書くことを目指すと、変数変換を2回を経て最終的に以下の形に変形できることが知られております。

　　　　 z = \int_{\infty}^{u} {\frac{1}{\sqrt{4u^{3} - g_2u - g_3}}du}

　楕円積分の形となっていますが、実はこれもワイエルシュトラスの標準形と呼びます。 g_2と g_3は変数変換によって生まれた各根の慣れの果てと思ってください（笑）これについても繭野孝和著の『楕円関数論への入門』に掲載されていますので、計算過程とか気になった方はそちらをどうぞ（と言っても、2回目の変数変換辺りから訳わかんなくなりましたが……）

ワイエルシュトラス標準形２： \mathscr{P}関数の微分方程式

　さて、一応はワイエルシュトラス標準形を導くことができましたが、このままだとまだ暗号関係での取り扱いが難しいです。そこで、先程求めたワイエルシュトラス標準形の逆関数――つまり楕円関数として、楕円曲線暗号の陰の主役とも言えるワイエルシュトラスの \mathscr{P}関数を登場させます。

　このヤケに難しそうな名前の関数ですが、三角関数で言う cosにあたる導関数の \mathscr{P}'について以下の方程式が成り立つのです……

　　　　 \mathscr{P}'^{2}(z) = 4\mathscr{P}^{3}(z) - g_2\mathscr{P}(z) - g_3

　こ、これは！！！

来ました！　この章の冒頭の方で示したワイエルシュトラスの標準形

　　　　 y^{2} = x^{3} + ax + b

と形式がかなり一致します！　なので、どうにかして \mathscr{P}をx座標、 \mathscr{P}'をy座標とした実平面の上に、この微分方程式と同型な代数曲線を構成してやれば楕円曲線が描ける事が分かります！

　という訳でこれについてもこちらで説明するべきですが――、これもちょっと時間の都合上で間に合わなかったので、説明を非常に分かりやすく説明してくださっているtujimotterのノートブログ様の記事に説明を譲りたいと思います……（他力本願で申し訳ない！）