巨大な素数の効率的な作り方

C++の勉強がてらにRSAのプログラミングをやってみたのですが鍵生成の時間が一定しません。
試しに512bitの素数を50回作ってみると、平均79秒、標準偏差71秒、最長343秒、最短9秒でした。
ここまで標準偏差が大きいと当然50回の平均を取ったところで、一定の値は取りません。
現在はオイラーの公式n^2+n+41で素数を出していますが、もっと効率良く素数を出す方法はないのでしょうか？

No.3 ベストアンサー 回答者： fujeira 回答日時： 2003/03/17 17:59

（１）

大きな素数を見つけるには、Miller-Rabin法を使うのが良いです。確率的素数判定法の１つです。厳密に素数であることを判定できるわけではありませんが、現実的に十分な信頼度で、高速に判定することができます。
なお、n^2+n+41で生成した数の中から素数を探すと、素数に偏りができてしまいます。512bitの一様な乱数を発生させて、素数かどうか判定するようにしてください。

（２）
整数の四則演算にも使えるプログラミングテクニックはいろいろあります。ビットの立っている位置を探すなら、以下の方法も使えるでしょう。

int mrb[256]; //1が立っている最も下位のビット位置をあらかじめテーブルに入れておく。
unsigned char u;

for( ;u;u&=u-1) { //最も下位の1をクリアする。
　mrb[u] ← 1が立っているビット位置を順番に取り出せる。
}

この回答への補足

回答ありがとございます。

実はNo.2の回答を頂いた後、色々と調べていたところ、モンゴメリ乗算という方法を見つけたのでそれのコーディングをしていたのですが、それがようやく今日できました。
x*y (mod a)という計算がaが常に一定で、繰り返し行われるとき効率良く計算できるという方法です。
これをa^b (mod c)の計算に組みこみ、さらにフェルマーの判定式の前に、10個ほどの小さな素数で割り切れるかどうか確認するようにしたところ、50回素数を生成した結果が平均7秒、標準偏差6秒、最長28秒、最短1秒になりました。
二秒程度で鍵ができあがるときもあれば、一分以上かかる場合もあるといったような状態ですが、一応使い物にはなるレベルにはなってきました。

> Miller-Rabin法を使うのが良いです。
最初はMiller-Rabin法でやっていたのですが、あまりにa^b (mod c)の計算速度が遅かったので簡略化のためフェルマーの小定理での判定法にしてしまっていました。
改めて調べて見ると、カーマイケル数の問題があるようなのでMiller-Rabin法に戻そうと思います。

> n^2+n+41で生成した数の中から素数を探すと、素数に偏りができてしまいます。
そこまで考えていたわけではありませんが、何となく使える素数が少なくなりそうと思い、現在は2n+1型の素数生成式しています。
オイラーの公式と比べて平均時間が5割以上、標準偏差が2割程度増えてしまいましたけど。
50回の平均の結果が標準偏差と同程度という状況を何とかしたかったんですけど、2n+1型で取ってくるしか無いのでしょうか？

> int mrb[256]; //1が立っている最も下位のビット位置をあらかじめテーブルに入れておく
値を格納するのに使っているのがlong型なので、この方法ですと32bitのCPU・OSの場合2^32個の配列が必要になるので、ちょっと無理ではないかと。
ビット演算で1バイトずつ取り出すと、結局余計な計算が増えますし。
現在は、
if (x[i / 32] & 1 << i % 32)
result = (result * temp) % c;
といった感じで取り出しています。

補足日時：2003/03/17 20:10

No.2 回答者： minimum 回答日時： 2003/03/06 23:49

以前調べたことがあるので、その資料を漁りつつ…

a^b (mod c) 効率的な方法ですが
どこまで効率化されているか解らないので
小さな所から(知ってるところまで)説明させて頂きます。

まず、
mod(X*Y,c)=mod(mod(X,c)*mod(Y,c),c)
が成り立つので掛け算をするたびに
cで割っておくことができます。

次に
a^b (mod c)　ですが、
まずbを２進数で表しておきます。
次にa^1,a^2,a^4,a^8,a^16,a^32,a^64,…を
順に求めていきます。
そしてbのビットが立っている所のみを
掛け合わせるとa^bが出来上がります。

この回答への補足

アドバイスありがとうございます。

教えていただいた方法はよくPowerModという名前で関数化されている方法ですよね。
残念ながらすでにその方法にしています。
一応ソースを載せておきます。

temp = a % c;
result = 1;
for(i = 0; i < <bのビット数>; i++){
　if(<bのi番目のビットが立っているか調べる>)
　　result = (result * temp) % c;
　temp = (temp * temp) % c;
}

1つの変数では512bitもの数は扱えないので、C++のオーバーロードを使い、配列を1変数として加減乗除をできるようにしているのですが、ネックになっているのは除算・剰余算のようです。
加減乗については細かな調整はできたとしても、元々考えられるアルゴリズムのパターンもたかが知れてますし、根本的なスピードアップは難しいと思っています。

除算・剰余算の現在のアルゴリズムは以下のようになっています。

1.除数が非除数より小さいのなら、除数を左へビットシフトし、非除数とビット数をそろえる
　除数が非除数より大きいのなら、終了
2.ビットシフト後の除数が非除数より小さいのなら、非除数から除数を引き、商の対応する位置にビットを立てる
3.除数を右へ1つだけビットシフト
4.2～3を除数のビット数が元に戻るまで繰り返す
5.残った除数が剰余に対応

この方法では下手をすると500回以上も2～3をループすることになるので、相当時間がかかっていると思います。
ここを何とかできれば相当時間を縮められるはずなのですが。

補足日時：2003/03/07 20:42

No.1 回答者： hjoshua 回答日時： 2003/03/05 01:29

n^2 あたりの計算を効率化したらどうでしょうか。

ただの巨大な数のかけ算をすると、O(n^2)かかってしまいます。平方根の計算方法がありますね。それの逆をやれば多分O(n log n)位になって安定するのではないでしょうか。

この回答への補足

アドバイスありがとうございます。

質問文では詳しく書いていませんでしたが、素数は次のようにして求めています。

1. n^2+n+41の結果が欲しい桁数になるようなnの範囲を求める
2. 1で求めた値の範囲で乱数nを発生させる
3. n^2+n+41を計算する
4. 3の結果が素数か評価し、素数でないなら2へ戻る

最終的にはこうして得られた素数を2つあわせて、RSA用の鍵とします。


主に時間がかかっているのは4の素数の評価のためにa^b (mod c)を計算しているところであって、n^2の計算は微々たるものなんです。(aは1桁の小さな数ですが、b,cは3で出した素数を元にした巨大な数です)
もちろんa^b (mod c)もできるだけ効率化しようとしており、組み上げた当初の半分くらいの時間でできるようになったのですが、私の頭ではそろそろ限界です。
効率的な剰余の計算方法があればいいのですが。

今回どうにかしたいと思っているところは、2～4のループ回数です。
このループ回数をある程度一定の範囲に収まるようにしないことには、たとえ計算速度自体を上げれたとしても、下手すると何分も素数が見つからずに延々ループを繰り返してしまうので・・・
なにか高確率で素数を出すいい方法はないのでしょうか？

補足日時：2003/03/05 18:17