高速な素因数分解アルゴリズム(Lenstra の楕円曲線法)を実装する

はじめに

春休みに入ってからというもの,腰痛が悪化する一方です.

TL; DR

Lenstra の楕円曲線法とは?

素因数分解アルゴリズムの一種で,現時点で見つかっているアルゴリズムの中では3番目に高速だと言われています.Hendrik Lenstra によって考案されたらしいです.

この記事では,Lenstra の楕円曲線法(以下 ECM と略称します)について,アルゴリズムや実行結果,原理などを簡単にまとめていこうと思います.

前提

くどくならない程度に前提となる数学的基礎を紹介しておきます.興味のない人は読み飛ばした方がいいかもしれません.

方程式  y^2 = x^3 + ax + b *1 の体  K における解全体がなす集合,すなわち
 E' = \{ (x, y) \in K \times K ; y^2 = x^3 + ax + b\}
に,無限遠点としての役割を持つ点*2  O を追加して得られる集合,
 E(K) = E' \cup O
は,Abel 群をなします.「点同士の足し算」として加法が突っ込まれた群です.この群を  K 上の楕円曲線と呼びます.足し算の構成法や Abel 群をなすことの証明に関しては,インターネットや書籍に優れた説明があるので,そちらを参照していただければと思います.

さらにまた,同じ点同士を繰り返し足すことでスカラー倍算も導入できます.すなわち,ある点 P \in E(K)に対し,
 nP = \underbrace{P + P + \cdots + P}_{n \: times}
という計算も定まります.

実数体上で考えた楕円曲線のグラフの一例を以下に示しておきます.

f:id:villach:20190317205727p:plain

ECM の手順

ECM は以下のような手順で構成されています.以下, n を実際に素因数分解したい合成数とします.以下の手順は  n の自明でない約数を求めるものなので,完全に素因数分解するには素数判定をはさみつつ再帰的に手順を実行する必要があります.

  1. 楕円曲線  E(\mathbb{Z}/n\mathbb{Z}) を適当に決め,その上の点  P をランダムに一つ定める.
  2. 適当に  k を決める.
  3.  P k 倍点である  kP を計算する. kP を計算する過程で無限遠 O に飛んでしまったら,そこから  n の非自明な約数が求まる(原理は後述).もし  kP が普通に計算できたら,1. か 2. に戻ってパラメータを選び直し,再試行する.

要は  kP を計算しているだけです. kP の計算過程で  n の非自明な約数が求まるまで,神様に祈りつつ試行を繰り返します.楕円曲線ガチャみたいなものですね.

また, k としては極めて大きな数が使用されるので,この計算にあたっては繰り返し2乗法*3などを用います.

実行結果

大きな数の取り扱いが大変だったので,Python 3 を用いて実装しました.ソースコードはチラ裏レベルのひどい代物なので,実行結果だけの記述とさせてください.

以下が実行結果の一例です.ちょっと出力がシンプルすぎる気もしますね.

f:id:villach:20190317211714p:plain

実際,29130502385925287 = 96600277 * 301557131 です.このくらいの数なら総当たりでも計算可能かもしれませんが,以下の場合はどうでしょう.

f:id:villach:20190317212029p:plain

このくらいの桁数になると,総当たりで計算するのは(少なくとも PC 使用かつ短時間では)不可能でしょう.

計算時間については,手元の PC でだいたい以下のような結果が得られています.入力に使った2つの合成数は,大きな素数2つの積で書かれる,いわば「タチの悪い」ものです.

  • 25227563887000340779(20桁)……約2秒
  • 119798122324797418809504766589(30桁)……約44秒

もっとも,ECM がかなり運ゲー色の強いアルゴリズムであること, k の選択についてあまり最適化をしていないことを考えると,もっと速くなる余地は十分あると思います.逆に,実行状況によってはもっと遅くなることも考えられます.

原理

 n が非自明な素因数  p を持っているとし,位数(要素数 p の体  \mathbb{F}_p 上の楕円曲線  E(\mathbb{F}_p) について考えます.

いま,ステップ 2. で選んだ  k が,幸運にも  E(\mathbb{F}_p) の位数の倍数となっていたとしましょう.このとき, E(\mathbb{F}_p) 上の任意の点  P は, k 倍したら必ず無限遠 O になります*4

 kP無限遠点であるということは, kP を求める式で登場する分母(仮に  d_k とします)が  \mathbb{F}_p 上で  0 となることに他なりません.いいかえると, p d_k を割り切ります.よって, p \mathrm{gcd} (d_k, n) を割り切り,わずかな例外を除いて  n の自明でない約数が得られます.

また,「 k E(\mathbb{F}_p) の位数の倍数である」という条件を満たすためにどうすればよいかを考えると, k をどう選べぶべきかもおぼろげにわかってきます.すなわち, k は約数をたくさん持つように選べばよいのです.これについては流儀が複数あるようです.

 

上述した内容からわかる通り,ECM E(\mathbb{F}_p) の構造に依拠しています.実は,ECM の前身(?)として,Pollard の  p - 1 法が存在します. p-1 法では, (\mathbb{Z}/p\mathbb{Z})^{*} の構造を利用していました.これは  p によって完全に固定されてしまうため, n の因数が「都合よく」なければ,その時点で諦めて試合終了でした.ECM は,代わりに楕円曲線  E(\mathbb{F}_p) を用いることで,仮に試行に失敗しても, E を取り換えて再チャレンジできるようにしています.これが ECM の偉大な進歩なんだそうです.

ECM では, E(\mathbb{F}_p) の位数が重要になってきますが,これについては Hasse の定理 なる定理が存在して, |\#E(\mathbb{F}_p) - p + 1| \le 2\sqrt{p} となることがわかっています.つまり,位数が  p + 1 を中心として  4\sqrt{p} の範囲に散らばっているということです.さらにまた,この「散らばり具合」に関しても 佐藤-Tate 予想という予想があって,そこまで偏らずいい感じに散らばっているらしいと信じられています.したがって, E(\mathbb{F}_p) の位数に望みを託す賭けも,そこまで分の悪い勝負ではないということが分かります.

今後の方針

まず,単純に実行速度を追求したいと思います.具体的には,パラメータの選択手法を洗練したり, kP の計算方法を見直したり*5といった手段が考えられます.

また,背景にある数学的理論をより深く知り,見通しの明確化を図るのも一つの目標です.

ECM の実装や理論面の学習は半分くらい高専の授業の枠内で行っているので,授業で得る知見や教員からのアドバイスも活かしていきたいですね.

参考文献

*1:Weierstrass 標準形と呼ばれます.

*2:厳密には射影幾何学を使って考えます.

*3:競プロ界隈の人にはおなじみだと思います.

*4:この辺は認識が曖昧なのですが,Lagrange の定理からしたがうそうです

*5:繰り返し2乗法の他にも,移動窓法や符号つきm進展開法といった手法が存在するようです.

女子小学生でもわかる剰余環の定義

はじめに

この記事きじんで,数学すうがく興味きょうみをもったおんなは,ぼくへ連絡れんらくしてください.個人こじんレッスンをしてあげます.

 

大人の方へ:

厳密性の話はしないでもらえるとありがたいです.

加算・減算・乗算の3つの演算を行える集合を,と呼びます.たとえば,整数全体は環です(整数同士の和差積は全て整数になります.一方,整数同士の商は整数になるとは限りません.たとえば,1/2 は整数ではありません).

同値関係

以下では R を環とします.
R の要素2つ (a, b) を引数として取り,TRUE/FALSE のどちらかを返す関数 f(a, b) を考えます.たとえば,R が整数集合なら,a < b のとき TRUE を返し,a ≧ b のとき FALSE を返す,といった関数などを挙げられますね.
こうした関数を,特別に2項関係と呼びます.また,2項関係は,関数の形で書く代わりに,適当な記号(たとえば ~ とか)を取り,a ~ b と書いたりします.a < b,a > b,a != b なんかは全て2項関係です.

ここで,以下の3条件を満たす2項関係 ~ を特に同値関係と呼びます.

  1. R に含まれる全ての要素 x について,x ~ x
  2. x ~ y かつ y ~ z なら,x ~ z
  3. x ~ y なら y ~ x

たとえば,整数集合での関係 = (等しい)は,明らかに同値関係です.整数 n, m, l について,n = n ですし,n = m かつ m = l なら n = l ですし,n = m なら m = n です.

また,ある整数 k を決め打ちして,「x と y は,それぞれ k で割った余りが等しい」ことを x ≡ y (mod k) と書きますが,これも同値関係です.

剰余環

さて,同値関係を導入すると,「同値な要素は全部ひとまとめにして考える」ということができますね.

たとえば,先ほどの「k で割った余り」について考えます.k = 3 とすると,整数は全て

  • 3 で割った余りが 0 になるやつ
  • 3 で割った余りが 1 になるやつ
  • 3 で割った余りが 2 になるやつ

に分類されます.すなわち,{ …, -6, -3, 0, 3, 6, … },{ …, -5, -2, 1, 4, 7, … },{ …, -4, -1, 2, 5, 8, … } の3分類ができますね.

ここで,それぞれの分類のうち,代表者を一つずつ選出します.別になんでもいいんですが,簡単さのために 0, 1, 2 を取り出しましょう.つまり,

  • 3 で割った余りが 0 になるやつの代表 → 0
  • 3 で割った余りが 1 になるやつの代表 → 1
  • 3 で割った余りが 2 になるやつの代表 → 2

みたいに考えます.

さて,こうすることで,代表者のみを集めた集合を作れます.{0, 1, 2} ですね.この集合を剰余環と呼びます.要するに国会みたいな感じですね.いまは例として 3 で割った余りのことばかり考えていましたが,同値関係として別のものを取れば,いろんな剰余環を作れます.

うれしさ

剰余環では,同値な要素を「ひとまとめ」にして,代表者だけを集めて考えることができます.
これによるうれしさは,間接民主主義を導入するうれしさと一部共通すると思います.古代ギリシャのポリスみたいに,みんなが顔見知りであるほど人数が少なければ,全員の意志を考える直接民主主義が可能ですね.しかし,現代国家のように,人数が爆発しているところでは,代表者を選出し,直接的にはそうした代表者の意志を考える間接民主主義が必要になってきます.

ましてや数学では,集合の要素が無限に存在することがたくさんあります.ここで,手に負えるサイズまで集合を「圧縮」し,しかもその集合をきちんと有意義に扱える(剰余環も環なので,演算ができます!),というのは大変うれしいことです.

特に現代では,有限な演算しかできないコンピュータにたくさんの処理をさせることが多いので,剰余環の効用はかなりありがたいです.先ほど考えた「3で割った余りについての剰余環」とかはかなり使われています.もっとも,こいつは環じゃなくて体になっちゃうんですが.環のままじゃイカン!ってことなんでしょうか.

第26回コンピュータフェスティバル競技部門に参加しました

宇部高専で開催された第26回コンピュータフェスティバルの競技部門に参加し,なんと優勝を勝ち取ることができたので,嬉しさが残ってるうちに記事をしたためておこうと思います.

参加記事は誰かが部のブログに書いてくれると思うので,主に技術的なことを残しておきます.

基本的な気持ち

最良優先探索(priority_queue を使った探索)で5ターン分先読みし,各状態について評価し,もっとも評価値の高かったものを採用するようにしました.

探索結果は保持せず,ターンが進むごとに一から探索しています.

評価関数

以下の5つを評価用のメトリクスとして採用しました.

  1. タイル点による利得
  2. 自チームで領域を形成することによる利得
  3. 相手チームの領域を破壊することによる利得
  4. 自チーム両エージェント間の距離
  5. 移動先での自由度

1. 2. 3. は特に説明しなくても大丈夫だと思います.付け加えるとすれば,領域点よりタイル点の方に多少加重しました.これは理論的な根拠があるわけではなく,テストプレイを積み重ねての実感に依拠しています.領域点を重視すると,いきおい「大風呂敷」を広げがちになり,安定した得点を得られなくなります.

4. についてですが,エージェントがくっついてもいいことはない(とりうる行動が制限される,相手が遠くで策動しているのを防げない,などの弱点ばかりが目立つ)ため導入しました.5. は,自チームのタイルがあらかた置いてあったり,隅っこだったりするところにあまり行かないようにするため設定しました.
これらのメトリクスは,マップのタイル得点の状況に依存しない数値なので,重みを固定した場合,全体として高得点なマップにおいては 1. 2. 3. との不均等が生じてしまいます.そこで,マップのタイル得点の平均値を計算し,それに応じて 4. 5. の重みを動かすようにしました.

また,体感として,タイルによって得点の上下が激しいマップほど領域点が大事になる気がしたので,タイル得点の標準偏差を求め,それがしきい値を超えた場合に 1. の重みを増やすようにしています.

対戦前の時点で,ぼくとソルバで対戦した場合は9割方ソルバが勝つ,くらいの強さまで持っていけていたはず……です.

対戦の結果

1戦目が一番大変でした.与えられたマップが線対称性を前提しておらず,読み込みの時点でエージェント座標をバグらせてしまいました.また,この混乱により,最初の2ターンはロスしました.最終ターンでなんとか逆転して僅差の勝利です.

以降は,上記の不具合も修正し,ソルバも順調に動いたので,基本的にうまく戦えたと思います.たまに,人間からすれば不合理に思える手(-99点のタイルを占領しにいく,序盤10ターンの間ひたすら前進するなど)もありましたが,ソルバを信頼するように努めていました(こう書くと美しいですが,実際には人力モードへの切り替えを導入していなかっただけです).また,コリジョンが連発することもあり,特に決勝では10ターンくらい睨み合っていましたが,この辺も柔軟に対応できるようシステムを作っていなかったので,我を通しつづけた感じです.

所感

この辺の考察を昨年9月の段階でやっておきたかったですね.

基本的なアルゴリズム自体はプロコン本選のものと一緒なんですが,当時は実装力が弱かった(今が強いというわけでもありませんが)ので,かなりひどいコードでした.冬休みの期間にその方面のコードを読み込み,いい感じのプラクティスを吸収できたのが活きてきたんだと思います.

はじめての制御工学:第14講

内容

ループ整形法について.

フィードバック制御系の開ループ伝達関数 L(s)が持つ周波数特性を好ましいものにするため,ループ整形法なる手法を用いることがある.以下に詳説する.

具体的に必要な特性は,たとえば以下のようになる.

  • 定常特性:低周波数帯域でのゲインが小さくなる.
  • 即応性:ゲイン交差周波数が十分高くなる.
  • 減衰性:位相余裕 PMが十分大きくなる.
  • ロール・オフ特性:高周波数帯域でのゲイン変化が急になる.

これらを満たすために,位相遅れコントローラと位相進みコントローラの2つを利用する.伝達関数はそれぞれ C(s) = \frac{s + \omega_1}{s + \omega_2}, \omega_1 \lt \omega_2 C(s) = \frac{\omega_3}{\omega_4} \frac{s + \omega_4}{s + \omega_3}, \omega_3 \lt \omega_4である.一般には,これらを複数つなぎ合わせてコントローラを構成する.

ロール・オフ特性について.通常,フィードバック制御系には観測ノイズ n(t)が制御対象からの出力に混入する.観測ノイズの影響を少なくするには,観測ノイズが存在する周波数帯域での開ループ伝達関数のゲインを急激に落とす必要がある.

感想

ようやく一通り終わりました.なせばなるものですね.

とはいえ,一周目なこともあり,詳細な計算を追っかけたり,付録に回されている導出過程を読んだりがほとんどできていません.このあたりは二周目以降でフォローしていきたいと思います.

はじめての制御工学:第13講

内容

Nyquist の安定判別法について.

系の開ループ特性を考えたとき,伝達関数が安定となるようパラメータを選べても応答の振動が激しくて困ったことになる場合がある.実用上十分なほど安定することを安定余裕があると呼び,そうでない場合(めっちゃ振動する場合など)は安定余裕がないとか小さいとか呼ぶ.以下,系が安定余裕を持っているかどうか判別する方法について考える.

フィードバック制御系の分母多項式 N_p(s) N_c(s) + D_p(s) D_c(s)で表され,この根は閉ループ極と呼ばれる.ここで, N_p(s), D_p(s)はそれぞれ P(s)の分子及び分母であり, N_c(s), D_c(s)はそれぞれ C(s)の分子及び分母である.以前触れた通り,閉ループ極の実部が全て負であれば系は内部安定となる.

さて,いまフィードバック制御系の4つの伝達関数の分母に出てくる 1 + P(s)C(s) N_p(s), D_p(s), N_c(s), D_c(s)で書いてみると,分母は D_p(s) D_c(s)となる.ここでその根を開ループ極と呼ぶことにする.

このとき,系の設計にあたっては開ループ極のうち不安定なものの個数 Zが知りたいのであって,閉ループ極のうち不安定なものの個数 Pは既知であることが多い.Nyquist の安定判別法は, Pがわかっている状態で Zを与える手法である.

(導出はサボったが)Nyquist の安定判別法は以下の手順からなる.

  1. 開ループ伝達関数のベクトル軌跡 L(i \omega) = P(i \omega) C(i \omega)を描く.
  2. 描いたベクトル軌跡と実軸対称な軌跡を描く.これは \omega: - \infty \to 0として軌跡を描くことにほかならない.こうしてできた軌跡と先のベクトル軌跡を合わせて Nyquist 軌跡と呼ぶ.
  3. Nyquist 軌跡が複素平面上の点 -1を時計回りに回る回数を 1,反時計回りに回る回数を -1としてカウントしていき,その合計を Nとする.
  4.  N = Z - Pである.したがって, Z = N + P 0なら系は内部安定である.

さて,実際には P = 0となるようコントローラを設計することが多い(安定な制御対象に安定なコントローラを付ける場合).このとき,簡略化された Nyquist の安定判別法を利用できる.具体的には,Nyquist 軌跡が点 -1を常に左手に見つつ原点へ収束するなら系は内部安定であり,そうでなければ系は不安定となる.

最後に,安定余裕の測定について考える.ここでは簡略化された Nyquist の安定判別法について限定する.Nyquist 軌跡が点 -1を十分な距離を保ちつつ左手に見ていれば系は余裕を持って安定となる.したがって,ベクトル軌跡と点 -1との距離を反映する指標が安定余裕の判定に役立つことがわかる.

ここで,ゲイン交差周波数 \omega_{gc}と位相交差周波数 \omega_{pc}なる値を導入する.それぞれ, |L(i \omega_{gc}) = 1|, L(i \omega_{pc}) = - \pi)を満たすような角周波数である.このとき,位相余裕 PM = \angle L(i \omega_{gc}) + \piが大きいほど系は安定余裕を大きく持ち,ゲイン余裕 GM = \frac{1}{|L(i \omega_{pc})|}が大きいほど,系は L(s)の増大による安定性の喪失を来しづらくなる.ゲイン交差周波数と位相交差周波数はボード線図をじっと睨むことで得られる.

感想

大変遅くなりました.そして長くなりました.

最後の方になるとダラダラしてしまってダメですね.

ところで,制御工学をサボってるうちに試験週間に突入してしまいました.もっとも試験勉強もサボってるんですが.本格的に試験が始まる前に制御工学を終わらせてしまいたいですね.あと1講!

はじめての制御工学:第12講

内容

ボード線図と周波数伝達関数について.

高次の伝達関数について,部分分数分解をほどこして個別に周波数特性を求めてから合成してもともとの周波数特性を求められる.

2次遅れ系ではゲイン K = 1でもゲイン線図が 0デシベルを超えることがある.これを共振と呼ぶ.

ゲインが -3デシベルに達する周波数をバンド幅と呼び,系の入力追従特性の指標となる.

ステップ応答について考える.ステップ応答は(Fourier 変換を見ることで明らかに)無限の周波数成分を含んでいる.高周波数帯域についてはゲインが負となるため,ステップ応答はすぐさま1に収束することなく(ステップ信号と同一でなく),適当な時間を経たのちに収束する.

系の伝達関数 G(s)がわかっている場合, G(i \omega)なる関数を周波数伝達関数と呼ぶ.周波数伝達関数の大きさはゲインと関係し,偏角は位相差を表す.したがって,周波数伝達関数を見ることで系の周波数特性を分析できる.

周波数伝達関数 \omegaを変数とする複素平面上のベクトルと見なせるから, \omega 0から \inftyまで変化させたときにベクトルがどう動くかを考えると周波数特性の分析に役立つ.ベクトルの軌跡を(そのまんま)ベクトル軌跡と呼ぶ.ベクトル軌跡もボード線図と同じく周波数特性を表す図だが,一つの図に大きさと偏角を同時に描きこめるのが利点である.

感想

大遅刻ですね.ごめんなさい.2日も空いてしまいました.

追記

 G(s) s = i \omegaを代入しているのは,要するに Laplace 変換のかわりに Fourier 変換をしているのに等しいです.通常 Fourier 変換は積分範囲を \mathbb{R}全体で取りますが,入力信号が t \lt 0 0となることから正実数全体での積分である Laplace 変換に機械的に代入することが正当化されます.

はじめての制御工学:第11講

内容

周波数特性について.

 A\sin(\omega t + \theta)なる正弦波を入力した際の応答を周波数応答と呼ぶ.

周波数応答の特性(振幅,位相など)はボード線図を書くことで把握できる.また,応答の振幅が Bとなるとき, 20 (\log A - \log B)をゲイン(利得)と呼び,デシベルの単位で取り扱う.

感想

残すところあと3講のみですね.