1024ビット素数生成の難しさ

 (glitchcomet.com)
4 ポイント 投稿者 GN⁺ 2024-05-05 | 1件のコメント | WhatsAppで共有
  • 2048ビットRSA鍵に必要な2つの約1024ビット素数をRustで自作生成し、外部依存なしで乱数生成から大整数演算まで実装した実験
  • 単純な trial division は16ビットでは約40msで十分だったが、64ビットでも最適化後に6.4秒かかり、1024ビットへ拡張するのは難しかった
  • Fermatテストは高速だが pseudoprime を除外できない場合があるため、最終判定には k=10 のMiller-Rabinテストを使用
  • 標準整数型の限界を超えるために自作の BigInt を作り、bool配列からbyte配列、u64チャンク構造へ変更することで、1024ビット素数の生成時間は32分台から60〜90秒程度まで短縮
  • 最終実装は u64チャンクBigInt、高速な除算、小さな素数によるtrial division、候補値を +2 で増やす処理、16スレッド並列実行を組み合わせ、平均約40msで1024ビット素数を見つけたが、検証済みの暗号ライブラリではない

RSA用の1024ビット素数を自作する

  • 目標は RSA鍵生成に使える素数を自分で生成することだった
    • 2048ビットRSA鍵は2つの素数の積で作られるため、それぞれ約1024ビットの大きさの素数が必要になる
    • 課題は自然と 1024ビット素数生成に絞られた
  • 実験には3つの制約を設けた
    • コードはゼロから書き、外部依存を使わない
    • 外部ハードウェアやクラウドは使わず、AMD Ryzen 7 CPUと16GB RAMを搭載したノートPCを使う
    • 「妥当な時間」内に素数を生成する
  • 言語には最近学んでいた Rust を選んだ
    • 低レベルの概念を扱うのに十分近く、コード片を理解するには十分高レベルだと考えた

16ビットと64ビットで見えたtrial divisionの限界

  • 基本的な流れは、Nビット乱数を繰り返し生成し、素数判定を通過したら終了する方式
  • 乱数はRustの rand crate ではなく、Linuxの /dev/urandom を直接読んで作った
    • /dev/urandom はLinuxカーネルの CSPRNG にアクセスする疑似デバイスファイル
    • カーネルはユーザー環境からエントロピーを収集し、ChaCha20ベースの決定的ストリーム暗号を定期的にシードする
  • 16ビット乱数は先頭ビットと末尾ビットを 1 に設定した
    • 末尾ビット 1 は奇数を保証するため
    • 先頭ビット 1 は必要なビット範囲全体を使うための仕組み
  • 16ビットでは、3 から sqrt(num) まで割ってみる trial division だけでも約40msで素数を見つけられた
    • 実行例は Prime found: 44809、全体の時間は約0.038秒だった
  • 64ビットへ拡張すると、単純なtrial divisionは約30秒かかった
    • その後、6k±1 形式の候補だけを調べ、小さな素数リストで先に割ってみる方式に改善した
    • 改善後の64ビット素数生成時間は約6.414秒だった
  • 64ビットでも6秒かかるため、この方式では 1024ビット素数生成には到達しにくいという限界が明確になった

確率的素数判定への移行

  • 決定的アルゴリズムのAPR-CLとECPPを調べたが、数学的に複雑で、手が届く説明も少なく、実装対象にするのは難しかった
  • OpenSSLのソースコードとNISTの推奨を確認した後、RSAを含む実利用の場面で 確率的素数判定が広く使われていることを確認した
  • 以後のアルゴリズムは、数が「素数であることを証明する」のではなく、一定の精度で probable prime と判定する方式に変わった

  • Fermatテスト

    • Fermatの小定理は、p が素数で a が p で割り切れない場合、a^(p-1) = 1 mod p が成り立つという関係を使う
    • 単純な累乗は u128 でオーバーフローが発生するため、モジュラ累乗を実装した
    • pow() は指数を u32 として受け取り、u128 をさらに大きな指数に上げるとオーバーフローする可能性がある
    • 乗算自体も u128 の範囲を超える可能性があるため、一時的に64ビット数を u128 内に保存する方式で進めた
    • Fermatテストは高速だが、Fermat pseudoprime のため合成数を素数と誤判定する可能性がある
    • こうした合成数はまれでも十分な数があるため、Fermatテストだけでは信頼しにくいと考えた

  • Miller-Rabinテスト

    • Miller-RabinはFermatテストと同じ原理に基づくが、より強力な確率的素数判定アルゴリズムとして使われる
    • 実装では n-1 = 2^s × d の形に2の累乗を分離した後、複数の条件を検査する
    • a^d = 1 mod n
    • または、ある 0 <= r < s について a^(2^r × d) = n - 1 mod n
    • 128ビット実験ではFermatテストと同程度に、約0.042秒で素数を見つけた
    • Miller-Rabinの最悪時の誤り上限は 4^-k、大きな n では平均的に 8^-k 程度
    • k=10 のときの平均誤り確率の計算は 0.000000000931323% だった
    • これはコインを30回連続で投げてすべて表が出る確率である 2^-30 と同じだと比較した
    • 実際の暗号用途では、ランダムなbase選択や敵対的条件にさらに注意する必要がある

BigIntを自作する

  • Rustの標準整数型だけでは64ビットを超える十分大きな数を扱いにくく、任意精度整数（BigInt） の実装が必要になった
  • 外部のbigint crateを使わないという制約のため、BigIntも自作した

  • 試行1: 数字の桁配列

    • 最初は大きな数を10進数の桁配列として保存する方式を試した
    • 加算と乗算は筆算のように実装できたが、除算の実装で行き詰まり断念した

  • 試行2: boolベースの二進配列

    • 2つ目の方式は、数を0と1の配列として保存する構造だった
    • BigInt は [bool; 2048] 配列を使用する
    • 1024ビット同士を掛けると最大2048ビットの領域が必要になるため、2048ビットを確保した
    • 加算と減算はfull adder方式で実装した
    • 乗算は二進数の特性を利用して shift-and-add 方式で処理した
    • 除算は二進long divisionで実装した
    • この実装で初めて1024ビット素数を見つけることに成功したが、実行時間は約32分44.90秒だった
    • 技術的には目標を達成したが、「妥当な時間」という制約には合わなかった

  • 試行3: byteチャンク

    • bool配列の各 bool が1ビットではなく1バイトを占めることを確認した
    • [bool; 2048] は2048ビットではなく2048バイトを使う
    • その後、2048ビットを256バイト配列に保存する方式へ変更した
    • 加算、減算、乗算は大きな変更なしで動作し、除算はbyteチャンクをビット列のように扱うよう調整した
    • この方式で1024ビット素数の生成時間は 4分43秒まで短縮された

  • 試行4: u64チャンク

    • byteチャンク方式は、実質的には高い基数の桁を使うdigitベースのBigIntだった
    • 次の段階では、2048ビットを u64 チャンク32個に保存した
    • 各チャンクは1つの「桁」のように動作する
    • 2つの u64 チャンクを掛けた結果を保持するために u128 を使用した
    • この構造では、1024ビット数を10進数309桁ではなく u64 チャンク16個で表現できる
    • 1024ビット素数の生成時間は 60〜90秒まで改善した

ボトルネックの最適化

  • 簡単なベンチマークでは、binary実装とu64チャンク実装の差がはっきり出た
    • a + b と a - b: 5537.35ns → 123.57ns
    • a * b: 1292283.14ns → 842.32ns
    • a / b と a % b: 733446.76ns → 44440.12ns
    • a < b と a > b: 2506.02ns → 58.91ns
  • その後の最適化は主に、除算、乗算、Miller-Rabin内部の演算、候補生成ロジックに集中した

  • 除算

    • 最大のボトルネックは 除算だった
    • u64チャンク構造でも、従来の除算は依然として1ビットずつlong divisionを行っていた
    • Handbook of Applied Cryptography の598ページのアルゴリズムを参考に、radixベースのlong divisionを実装した
    • dividendの先頭3つの「桁」とdivisorの先頭2つの「桁」から、現在のquotientの「桁」を推定する方式
    • この実装は除算1回あたり約40,000nsを節約した
    • divisorが単一の u64 チャンクなら、u128 を使ってより直接的なlong divisionを行うよう特別扱いした
    • Miller-Rabinではこのようなケースが頻繁に現れる

  • 乗算

    • 乗算は中間結果保存用のBigIntをなくすようにループを組み替え、約2倍速くなった
    • 使用中のチャンク数を計算し、0でないチャンクだけに対してループを回すように変更した
    • BigIntはほとんどの場合1024ビット以下の数を保存するため、2048ビット領域の半分が空いていることが多い
    • KaratsubaやFFTベースの乗算も検討したが、自作するには複雑で、現在の乗算が十分高速になったと判断した

  • Miller-Rabin内部の最適化

    • Miller-Rabin実装では、コストの高い演算を減らすことに集中した
    • x = mod_exp(x, 2, n) の代わりに x = (x * x) % n を直接実行した
    • 最初の mod_exp() は単純化したインライン版に置き換え、関数呼び出しのオーバーヘッドを減らした
    • 偶数判定に num.is_even() を追加し、% 2 の計算を避けた
    • d / 2 は d >>= 1 に変更した
    • += 1、-= 1 は increase() と decrease() で特別扱いした
    • 特に is_even() と d >>= 1 はそれぞれ約70,000nsの効果を出した
    • 最終ベンチマークでは、u64チャンク最適化版が大きく高速化した
    • a * b: 842.32ns → 295.04ns
    • a / b と a % b: 44440.12ns → 831.77ns
    • a / 2: 75121.58ns → 60.89ns
    • a % 2 == 0: 78400.87ns → 21.65ns
    • a - 1: 103.15ns → 67.54ns

最終的な1024ビット素数生成器

  • 最終関数はまず /dev/urandom から1024ビット乱数を読む
    • 最上位ビットを立てて1024ビットの大きさを保証する
    • 最下位ビットを立てて奇数を保証する
  • その後、新しい乱数を毎回読み直すのではなく、候補値に 2 を足して次の奇数候補へ移動する
    • increase_by_2() はほとんどの場合、u64 チャンク1つの加算だけを行う
  • Miller-Rabinの前に、小さな素数リストで先にtrial divisionを行う
    • 最終コードでは最初の1000個の小さな素数を使用した
    • 小さな素数は単一の u64 チャンクに収まるため、高速な単一チャンク除算の特別処理を活用できる
  • この問題は共有メモリやスレッド間同期を必要としない embarrassingly parallel な形として扱える
    • 16個のCPUスレッドがそれぞれ素数を探し、最初に結果を送ってきたスレッドの値を使う
  • 最終実行例では約0.086秒のelapsed timeを記録した
    • CPU使用率は690%と表示された
  • 100回実行の平均は 0.04109 ± 0.00307 秒だった
    • 平均して約 40ms で1024ビット素数を見つける
    • 個々の prime_1024bit() 呼び出しはランダム性のため、約8msから約800msまで変動することがある
    • 並列実行で最速の結果を選び、ばらつきを緩和する

コードと限界

  • 全コードとリポジトリは github で公開されている
  • 議論リンクは hackernews と reddit にある
  • この実装は実際に暗号学的に安全だとは見なしにくく、目的も暗号用ライブラリの作成というより 学習と実装実験に近い

関連記事

1件のコメント

 
GN⁺ 2024-05-05
Hacker Newsのコメント

  • 大きな素数を見つける作業をプルーフ・オブ・ワーク関数の一部として使っていた暗号通貨がいくつかあり、8年ほど前には非常に高速な素数判定実装だけでもかなり稼げた
    しばらく riecoin マイニングソフトウェアの作者兼メンテナーをしていたが、理由はよく分からず、ただ素数が好きだったからだと思う
    この記事は高速な素数判定の最優先最適化であるMontgomery乗算を落としている: https://en.m.wikipedia.org/wiki/Montgomery_modular_multiplic...
    これは実用的な高速モジュラー累乗実装の基盤になる
    当時は学界にいて、今は Nvidia にいると認識している Niall Emmart が、本当にとてつもなく高速な GPU 多倍長整数ライブラリCGBNを公開した: https://github.com/NVlabs/CGBN
    いまだに自分の知る限り最速のバッチ・モジュラー累乗実装で、少しオタクっぽく感嘆するなら息をのむほどだ
    いつか、このおかげで小さな暗号通貨の生産を5年ほど支配していた話を書きたい。あと Python には pow(x, y, m) の3引数形式で x^y % m を計算する、かなり良いモジュラー累乗が入っている
    これを使えば、自分で実装したいときに Fermat あるいは Miller-Rabin の素数判定をとても簡単に作れて、かなり楽しい。自作したくないなら gmp ライブラリの mpz_probab_prime() も良い。もちろん gmp のほうが速いが、大きな素数で遊ぶときの2行の Fermat テストの楽しさにはなかなか勝てない

    • Niall は高速多重スカラー乗算に関する ZPrize 受賞提出物の1つにも関わっていた
      バッチ・モジュラー累乗と密接だが、素数モジュロではなく楕円曲線上で動く点が違う。CGBN の作業を引き継いだ結果だと見ている
      昨年、Stanford の暗号学ランチセミナーで良い発表をしており、スライドと録画もオンラインにある
      https://cbr.stanford.edu/seminarTalks/slides_20230526_niall_...
      https://www.youtube.com/watch?v=KAWlySN7Hm8
    • そうした暗号通貨がなぜそんな特注のプルーフ・オブ・ワーク関数を使っていたのか気になる
      暗号ではとにかく素数を使う、という漠然とした考えだけがあって、いつなぜ使うのか分かっていなかったのか、それとももっと深い理由があったのか知りたい
    • pow(x,e,mod) のおかげで Perl から Python に乗り換えた

  • 与えられた最大数の範囲があれば、Miller-Rabinを事実上決定的にするのは簡単だ
    その範囲内のすべての疑似素数をまとめて弾けることが証明された底を選べばよい
    リストも長くならない。Miller-Rabin は本当に強力だ

    • 1024ビット数の範囲ではその底が何なのか気になる
      オンラインでは答えを見つけられなかった
    • それに、単に素数を探しているだけなら、素数らしく見える候補を選んで決定的テストで確認すればよい

  • インラインアセンブリ1行で、多倍長整数の筆算式乗算は簡単になる: https://github.com/jcalvinowens/toy-rsa/blob/master/bfi.c#L4...
    過去に戻って C 言語で1つだけ変えられるなら、拡張乗算の概念を入れたい。Rust にもないのが惜しい。ハードウェア支援はどこにでもある。Cortex M0 は除算すらないが拡張乗算はある
    ずっと前に書いた、とても醜いおもちゃの RSA 実装から出てきたコードだ: https://github.com/jcalvinowens/toy-rsa
    Fermat テストだけでも持ちこたえられた理由は、素数が実際には素数でないとアルゴリズムが動かないからだった。Fermat テストは高速で、暗号化/復号を1回行えば、どちらかが Fermat の嘘つきである極めて小さい可能性も排除できる
    ただし、素数ではない P/Q 値でもメッセージを正常に暗号化・復号できる RSA 鍵ペアが存在しないと証明できるのかは分からない。実装としてはもちろん定石ではないが、答えを見つけたことはない

    • 興味深いことに、C には今や多倍長整数がある
      C23 で _BitInt(N) 型が追加され、たとえば 128 バイト型として _BitInt(1024) を使える
      ただしコンパイラ対応は限定的だ。Clang で N を 128 より大きく許可するには -fexperimental-max-bitint-width=N フラグを指定できる。N が 128 より大きい _BitInt(N) を除算するとコンパイラがそのままクラッシュするが、+、-、* は期待どおり動く
    • Zig ではこれは比較的簡単だ
      @mulWithOverflow 組み込み関数があり、結果とともにオーバーフロービットを返し、整数型は (u|i)65535 まである
      何をしたいかによって、オーバーフローを検知してからより大きい型に昇格させることも、先に昇格させてから必要に応じて切り詰めることもできる
      さらに別演算子 *| で飽和乗算、*% でラップアラウンド乗算をサポートしており、そうした意味論が必要なときに使える。それ以外のオーバーフローは、安全性チェック対象の未定義動作として Debug と ReleaseSafe ビルドモードでは panic になる
    • p と q が互いに素なCarmichael 数なら、RSA は依然としてメッセージを正常に暗号化・復号できる
      ただし p*q はより小さい素因数を持つことになり、因数分解しやすくなるので安全性は下がる
    • たいていの C コンパイラと Rust では、より大きい型にキャストしてから掛ければ、まさに欲しい機械語命令が生成されると理解している
    • Philip Zimmermann の1994年の元祖Pretty Good Privacy(PGP) は、既知のすべての16ビット素数で割るふるいだけを使っており、その表はエラトステネスのふるいで作っていた。その後に Fermat テストを適用していた

  • これにどれくらい時間がかかったのか気になる。
    学部の研究プロジェクトで大きな整数の乗算をやったが、ほぼ2学期かかった。Karatsuba、Toom-Cook、複素 FFT、いくつかの NTT、Schonhage-Strassen を実装した。
    素数はほとんど数学の魔法に近い。興味のある人には Silverman の A Friendly Introduction to Number Theory はすばらしい数学の入門書だ。
    ちなみにページのリンクは 40250519 ではなく 4025051 になっている

  • 良い記事。私も最近 [0] の初期バージョン向けに自分で 多倍長整数コード を少し書いたが、数学論文の高レベルな説明を実際の演算に落とし込む作業がどれほどもどかしいかを思い出した。
    ただし小さな異論がある。
    u64 の全範囲を使うなら、数は 2^64-1 進法ではなく 2^64 進法 だ。各ワードは 0 から 2^64-1 までの範囲を持ち、10進数の各桁が 0 から 9 までであるのと同じだ。
    [0] https://github.com/LegionMammal978/bigfoot-sim

  • 最後の最適化のように、失敗したとき新しい乱数を生成せず数を 2 ずつ増やすと、セキュリティが少し損なわれる。
    素数は一様に分布していないため、大きな 素数ギャップ のすぐ後ろにある素数へと偏る。

    • 調べているときにその話を読んだ。
      実行速度と素数のランダム性のあいだのトレードオフで、16 個のスレッドがそれぞれ乱数から開始して素数を見つけるまで競争するなら、十分なランダム性が追加されると判断して速度を選んだ。
      速度よりランダム性を重視するなら、+=2 を rng() 呼び出しに変えるのは簡単な変更だ。

  • 良い記事で、よく書けている。
    筆者は base-255 ではなく base-256 を意味していたのだと思う。

  • 1～2KB の数がいくつかある程度なら L1 キャッシュには十分収まるし、仮にそうでなくてもアクセス時間が約 3ns の L2 キャッシュがメガバイト以上ある。
    記事では L1 キャッシュミスのせいで RAM の読み書きを待っていたのではないかと述べているが、その後この点には戻っていない。
    また、これは素数生成だけを扱っているので RSA の落とし穴のほとんどを避けており、urandom は安全なはずだ。コードが正しく動作するなら、大きく問題になる点は多くない。
    RSA には避けるべき 弱い素数 に関する問題がいくつかあるが、ここで実際に問題になるほど一般的なのかは分からない。

  • 何十年も前の大学1年のプロジェクトを思い出す。
    プロジェクトパートナーであり友人でもあって、のちに卒業生代表になった人がアイデアを出し、中核となる数学を実装して、4096ビット RSA 暗号化 を作るというものだった。
    最終実装で素数生成がどれほど遅かったかを覚えている。PA-RISC ワークステーションで生成に約 20 分かかった。
    数学オタクだったその友人はプロジェクトが終わったあともコードの最適化を続け、素数判定や多倍長整数演算の実装に関する論文を読んでいたのを覚えている。
    たとえば筆算式の乗算で、ある数が 0 ならその乗算を飛ばして結果を 0 にするようにしたところ、ものすごい改善があった。

    • 遅いハードウェアでは 楕円曲線鍵 を生成するほうがはるかに良い。
      さもなければ長く待つか、将来まで持つセキュリティを犠牲にすることになる。

  • 下位ビットを 1 に設定するのは分かる。偶数は決して素数ではないからだ。もちろん 2 は例外だ。
    でも、なぜ上位ビットも 1 に設定するのかは分からない。素数や暗号に詳しいわけではないが、不必要に エントロピー 1 ビット を捨てているように見える。何か見落としているのだろうか？

    • 上位ビットが常に設定されていて、そのビットまで含めて素数をエンコードすれば、素数は常に同じバイト数でエンコードされる。
      可変長のバイトエンコーディングは、仕様が非常に明確で十分にテストされていない限り、異なるソフトウェア間でデータをやり取りする際に問題を起こしかねない。
      サーバー公開鍵に先行 0 がある場合の RSA ベース DHE で起きる問題を見ればよい。
    • 2 桁の数字を生成するのと同じことだ。
      最初の桁が 0 なら 2 桁の数字ではない。
    • 最上位ビットを 1 に設定すると 1 ビットのエントロピーは失うが、素数が十分に大きいことを保証できる。
      さらに言えば RSA では 2 つの素数を掛け合わせる。片方が 1024 ビットなら、記憶が正しければもう片方は約 200 ビットでも鍵に必要なエントロピーのビット数に達しうる。
      だから両方の素数を 1024ビット にしておけば、いくらか余裕もできる。
    • エントロピー 1 ビットを捨てるのは確かだが、それでも 1022 ビット残っている。
      誰かが 1024 ビット素数を求めているときに、1020 ビット素数でも大丈夫かと悩むよりはこちらのほうが安全そうだ。ふつう 00042 を 5 桁の数字とは見なさないのと同じだ。
      技術的には、正確にどこで使うかによって最適な選択は変わりうるが、この記事のやり方のほうがより安全なデフォルトに見える。
    • 確かに、50 ビット程度で終わらない素数を作るために 1 ビットのエントロピーを失うのは、十分に妥当な折衷に思える。