素数発見難度の評価
at NewPGen.exe(Paul Jobling) / PRP.exe(George Woltman) / Proth.exe(Yves Gallot)
文中に出てくるln(n)は自然対数関数を表します。
素数が大きくなるほど発見が難しくなる要素は以下の3つになります
1.数が大きくなるほど、素数の出現頻度が下がる
素数の出現頻度は桁数に反比例します
2.数が大きくなるほど、1回の計算時間が増える
計算（乗算や除算）にかかる時間はFFTを使用した場合、
桁数に対して、nln(n)ln(ln(n))に比例します。[n=桁数]
3.数が大きくなるほど、素数判定に必要な計算回数が増える
Proth判定などの素数判定に要するステップ数は桁数に比例します。
1,2,3を総合すると、桁数nの素数を見つける難度d(n)は基本的に
d(n)=cn3ln(n)ln(ln(n)) [1]
と表されます。(cは定数)
大雑把にいえば、素数発見難度は桁数のほぼ3乗に比例する、ということです。
もちろん、これは特殊な形式(メルセンヌ素数なども含む)をした素数に限っての話
です。一般の数に対応する素数判定はこれよりはるかに時間を要します。
以下、NewPGen.exe, PRP.exe, Proth.exeを組み合わせて使用した場合の
k*2n+1(または-1)の形式をした素数発見難度を考察してみました。
素数候補のふるい落としについての考察
素数探索で最初に行うのは、探索範囲内の数のうち、小さな素因数を持つもの、
例えば2,3,5,7,11,…で割れるようなものを除去していくことです。特に大きな
数に対しては探索範囲内全ての数にPRPテストやProthテストを実施するより、
その前に小さな素数で割リ切れるか、試し割を実施してPRPテスト対象数を減らし
ます。これによって探索効率が遥かによくなります。
ふるいというのは、試し割を1個1個の候補に対して行うのでは無く、探索範囲
の全候補に対してある素数pで割れる候補を同時に消去していくステップです。
ふるいによって、どの程度素数候補が絞れるのか、考察してみます。
以下、NewPGen.exeで、k*2n+1(または-1)のモードでふるいをかけた場合について
の場合です。kを固定する場合については、
(1) kの値によって、素数頻度の差が非常に大きいうえに、その頻度を正確に予測す
ることが比較的困難なこと。
(2) kを固定するモードではふるい速度が大幅に遅くなること
の理由で触れていません。ただ、大雑把な推定では、k固定モードではn固定モード
と比較し、次のステップであるPRPテストの候補対象が平均して10-15%程度多くな
ると見積もっています。
まず、底が2の場合、通常はkが奇数の場合だけを対象とするので、偶数であるk
は除去されます。残りは、奇数である全体の1/2のk値となります。
次に残った奇数について、最小の奇素数3で割った場合、余りは0か1か2のいずれ
かで、これは周期的に現われます。よって、残りのうち余りが0となる1/3の候補が
消えます。言い換えれば、(3-1)/3の候補が残ります。
最初から考えれば、全体の1/2*{(3-1)/3}=1/3の候補数となります。
残りの候補について、次の素数5で割った場合、余りは0,1,2,3,4のいずれかで、
それぞれが1/5の割合で出現しますから、残りのさらに1/5の候補が消えて行くこと
になります。言い換えれば、(5-1)/5の候補が残ります。
最初から考えれば、全体の1/2*{(3-1)/3}*{(5-1)/5}=4/15の候補数となります。
以下同様にあるpまでふるいを実施した場合、残る候補は以下のように考えられます。
(kmax-kmin)[{(2-1)/2}*{(3-1)/3}*{(5-1)/5}*{(7-1)/7}*…*{(p-1)/p}] [2]
[2]式をpの関数とおいて、プロットすると、以下の近似式が得られます。
R(p) = (kmax-kmin) / {1.7722ln(p)+0.1226} [3]
式を得るため、UBASICを用いて107までの素数で計算を行った結果を用いていま
す。この結果は、累乗の逆数無限和であるゼータ関数
ζ(s)= 1/1s+1/2s+1/3s+…
は以下の無限積でも表現できること、
= 2s/(2s-1)*3s/(3s-1)*5s/(5s-1)*7s/(7s-1)*…
また、自然数の逆数和
1+1/2+1/3+1/4+1/5+…+1/x
=ln(x)+γ [γ=0.577216…]
と近似できることからも頷けるかと思います。
[3]式における定数値は有限個の素数項の積がどこまでの有限逆数和に相当するかの
補正になります。
結局、ふるいをかけた時の残りの候補数は篩がけで用いた上限素数pの対数値（言い
換えれば桁数）にほぼ反比例する、ということです。
上記の[3]式が実際にどの程度近似されているか、以下の設定でNewPGen.exeを実際
に走らせて確認してみました。
k*bn+1 mode
base=2, n=10000, k min=1000000, k max=1999999
これらからも、上記の推定が十分な近似を示していることが分かります。
[この例ではR(p)は過小評価となっていますが、過大評価となる例もあります。]
NewPGen.exeでの最適なふるいの終了時点について
NewPGen.exeによるふるいは具体的にはPRP1回のテストに要する時間より短い
間隔で候補が減少している間はふるいを継続し、ほぼ同じ頻度になればそこでふるい
を停止する、ということになります。
ふるいを実際に実行すれば終了時点はおおよそ分かりますし、候補が消える頻度で
自動的にふるいが終了するように設定するモードもあるのですが、ここでは最適終了
時点とそれに要する時間を予測するための考察を行います。
まず、ふるいによって残る候補の数は上記[3]式によって与えられていますから、
これを微分することで、ふるいによって候補が減っていく速度を示す式が得られます。
[3]式R(p)の導関数をR'(p)とおくと、
R'(p)=-1.7722(kmax-kmin) / p{1.7722ln(p)+0.1227}2 [4]
[4]式がふるいで用いる素数pの時点での候補が減っていく割合を示します。
大まかに言えば、候補数が減少する速度はふるいで用いる素数pの大きさに対して、
p*(pの対数)2に反比例するということです。実際、私が計測したNewPGenの処理速
度などをもとに定数を設定するなど、[4]式を修正して、以下の式が得られます。
S = CV(kmax-kmin) / p{ln(p)+0.07}2 [5]
S : ふるいがけ素数pにおける１秒当たりに消える候補数
C : NewPGen処理速度によって決まる定数
C=2.6*107 (p<232の場合) 2.1*107 (p>232の場合)
V : パソコンの処理速度[単位GHz]
kmax,kmin : ふるいの範囲であるk値の最大値と最小値
p : ふるいで用いる素数p
ここでの注意点
ふるいの速度は上の式に挙げた以外にも、
1. ふるい開始の直後は前処理の影響と消える候補が極端に多いため、速度は遅くなり
ます。
2.ふるい落とし素数の大きさが232(約43億)を超えると速度が約20%低下します。
3.2.以外にもp値によって小さな速度変化が生じる場合があるようです。
4.kmax-kminの値、つまりふるい落としを実施する範囲の大きさによってpの増加速度
が変化します。範囲が広ければ若干遅くなり、10000000と50000では約10%の違いが
あります。[5]式ではkmax-kmin=10000000の場合をもとに定数を決定しています。
5.それ以外にも残っている候補数が多ければ、処理が増えるため遅くなり、少なくなる
と速くなる傾向があると考えられます。
上記のうち2.についてのみ[5]式で考慮しています。それ以外の要素は比較的小さいとみ
なして無視していますが、おおよその目安を知りたいのであれば、十分と考えています。
１個の候補が消えるのに必要な時間に関しては、[5]式の逆数で与えられます。
Tsie = p{ln(p)+0.07}2 / CV(kmax-kmin) [6]
Tsie : ふるいがけ素数pにおける１候補が消えるのに必要な時間 単位:秒
各変数および定数は[5]式と同じ
大きな素数を探す場合のふるいがけでは[6]式の方が使いやすいでしょう。
PRPテストにおける時間
ふるいの終了タイミングはPRPテストの速度にも依存します。そこでPRPテストに要
する時間を考えます。
先の素数探索難度で紹介した3つのうち2と3の要素によって決まるため、1回のPRP
テストに要する時間はターゲットとなる素数の桁数(べき指数)nとすると、
Tprp =A*n2*ln(n)*ln(ln(n)) [7]
Aは定数。
と表されます。実際はFFTサイズの変わり目などで不連続に処理速度が変わるのですが、
ここではそれは無視します。
PRP.exeを実際に走らせた結果から、定数は以下のように設定しました。
Tprp = Dn2ln(n)ln(ln(n)) / V [8]
Tprp : PRP.exe1テストに要する時間 単位 秒
D : PRP.exeの処理速度によって決まる定数 5.1*10-10
n : k*2n+1のnの値
V : パソコンの処理速度 単位GHz
以上の結果から、ふるいを終了させるタイミングは
Tsie = p{ln(p)+0.07}2 / CV(kmax-kmin) [6]
Tprp = Dn2ln(n)ln(ln(n)) / V [8]
この2つの式の値が等しくなるところになることがわかります。
ふるいに必要な時間はこの関係式から得た、pの値を単位時間あたりのpの増加速度
で割れば求められます。私の測定結果では
3.9*107 /sec./GHz (p<232)
3.1*107 /sec./GHz (p>232)
となっています。(NewPGen.exe V2.80より)
素数を探す際に必要なリソース（CPU時間）の見積り
実際に素数を探索する場合、どの程度の範囲を探索すればよいのでしょうか。まず、
素数発見難易度の1つ目の要素で紹介したように、素数出現頻度は桁数に反比例して少
なくなります。よって、ターゲットとなる桁数に比例して探索範囲を広げる必要があり
ます。k*2n+1の形式で、ターゲットとなる桁数と対応する2nの指数、および素数が1
つ存在する平均区間の対応は以下のようになります。
表1 素数の頻度
素数の存在は個々でみればランダムなので、各桁数に対応するk値の上記範囲を全て
探査したとしても、素数が必ず見つかるわけではありません。基本的に上記kの平均
間隔区間を全部探査して素数が見つかる確率は
1-1/e (eは自然対数の底=2.718281828…)
=約63.2%です。
よって、素数を確実に発見したい場合は、上記区間より広い範囲を探索することが必
要になります。この範囲の3倍を探索すると、約95%の確率で少なくとも一つは発見
できます。
以下の条件で、素数探索に要する時間を見積もって見ます
a.素数探索ソフトはNewPGen.exe , PRP.exe および Proth.exeを組み合わせて使用する
b.ふるい終了タイミングは対象となる数の大きさでの素数平均出現区間をふるいにかけ
たとして計算される、最適終了時とする。
c.素数のタイプはk*2n+1とし、NewPGen.exeで素数候補をふるいがけする場合はnを
固定したモードで実施する。
d.この推定は私のP3-750M(600M)Hzノートパソコンで実施した結果をベースとした。
計算式は上記の[6]および[8]を元にしている。
e.CPUの種類や条件により、ふるい速度やPRP速度が、計算量を単純に外挿できなく
なる場合があるが、この推定ではそれらは一切無視している。
まず、ふるいについて。
再度[6],[8]式を比較します。
Tsie = p{ln(p)+0.07}2 / CV(kmax-kmin) [6]
Tprp = Dn2ln(n)ln(ln(n)) / V [8]
両式の変数のうち、Vは共通ですから無視できます。
[6]式のうち、kmax-kminの値は探索する素数の桁数によって決まります。b.の条件より、
kmax-kmin = 0.723n と表されます。つまり、対象となる桁数によってのみ、ふるい上限
pの値が決定されます。実際にnに対応するpを直接計算するのは難しいので、nを決定
して、p値が合うようにいくつか代入して得られた数値表を以下に示しました。ふるいを
実施する範囲は特に指定していませんが、表1で示した範囲よりも広い範囲で実施するの
がよいでしょう。
表2 NewPGen.exeによるふるいに必要な時間
ここで、ふるい上限p値は単位1012(=trillion)です。
NewPGen.exeではp値の上限が260(約1150000*1012)となっています。
ふるいに要する時間（PRPに必要な時間と共通の単位です）
GHzs : 1GHzのマシンで必要な計算量 単位 秒
GHzh : 1GHzのマシンで必要な計算量 単位 時間
GHzd : 1GHzのマシンで必要な計算量 単位 日
GHzy : 1GHzのマシンで必要な計算量 単位 年
素数1個当りに残る候補数
表1で示した、素数となるkの平均間隔の区間において、ふるい終了時に何個の候補が
残るかを示します。
表より100万桁素数を見つけようとする場合は、ふるいの段階で1GHzマシンならば8
年以上かける必要があることになります。なお、NewPGen.exeにはふるいを複数のマシ
ンで分割して実施する機能があり、マシンの台数を集めることで、処理時間の短縮を図
れます。
PRPテストに必要な時間
ふるいの次に、残った候補に対してPRPテストを実施します。
[8]式に値を代入して得られた素数判定に必要な時間を表にしました。数によっては実際
の速度とやや異なる場合があります。
表3 PRP.exeによるPRPテストに必要な時間
50%とあるのは素数発見確率が50%に到達するのに必要な計算量で、素数が平均1個
存在する範囲のおよそ69％を探索した時点です。
90%とあるのは素数発見確率が90%に到達するのに必要な計算量で、素数が平均1個
存在する範囲のおよそ2.3倍を探索した時点です。
最終的には、このテストで素数らしいと判定された数をProth.exeなどで素数であるこ
とを確認するのですが、全体の時間から見ればわずかなので、ここでは考慮しません。
表2と3で示した時間の合計が素数探索に必要な実質的なリソース量となります。100
万桁素数の場合、発見確率50％まで到達するのにふるいの時間と合わせ、1GHzマシン
でおよそ200年要することがわかります。逆に言えば、200GHz分のCPUパワーを集
めれば1年で発見確率50％に到達できる、ということになります。
個人で大きな素数を発見しようとする場合、時間さえかければ20万桁、30万桁以上で
も発見可能ではないかとみています。続ける根気次第です。もちろん、上記で示したよ
うに最適化することが大切です。
以上
2004年05月02日 T.Nohara