- k を 2^n で割った余りでフィルタリングが可能。詳しくはこちら
- k=3m+2 と書ける場合は除外可能。
- 2m+1→3(2m+1)+1=6m+4→3m+2 となるので、自分より小さい k が大きいg(k) を持つから。g(2m+1)=g(3m+2)+2
- k=6m+4 と書ける場合は除外可能。
- g(6m+4)=g(3m+2)+1=g(2m+1)-1
- N/2までのg の最大値 = G(N/2) ( at k=k0 ) を既知とする。このときN/2からNまで偶数を除外して調べてもよい。g の最大値がG(N/2)と同じか小さかった場合、代わりに G(N/2)+1 を最大値とする。k=2*k0 でこの値となる。
これらのフィルタリングによって、例えば 3*2^24 で割った余りでフィルタすると全体の 7% (3507928/50331648) だけ調べれば良くなる。調べるべき余りの数のテーブルを作成して使う。3*2^24 の場合 unsigned int で 14Mb。
- N/4までのg の最大値 = G(N/4) を既知とする。k を2^n で割った余りで分類し、初めて初期値の1/4 以下になるステップ数 g0 とその時の飛びテーブルを用意しておく。1/4 以下にならない場合は計算できる最大ステップ数までにしておく。g(k) は g0+G(N/4) 以上にはならないので、これが現在の最大値より小さい場合はそれ以上計算しない。N/2 から N まで調べる場合、g の最大値の初期値として G(N/2) を与えておくと、だいたい半分くらいカットできる。なお 2^n は上記フィルタリングの値と同じにして、余りのテーブルと一緒にまとめて構造体に用意しておく。これによりメモリアクセスする番地が連続的になるのでキャッシュ性能が上がると思われる。カットしなかった場合でも、24〜36 ステップ程度進むので損にはならない(オーバーフローが起こらないようにする。24 bit の場合、初期値が 2^49 までOK)。
- 繰り返し部分で、x を 2^BITS で割った余りで分類し x → (x>>BITS)*a[x&(1>>BITS-1)]+b[x&(1>>BITS-1)] とする場合の飛びテーブルは連続アクセスとならないので、キャッシュに全部収まる程度が最速となる。手元のマシンでは BITS=10 程度(28Kb)。
これらを使って計算した場合、N=2^34 〜1.7*10^10 が 69秒(うち7秒が初期化)で終了。
このスピードだと世界記録の2^58 オーダーは50年、てことになる。並列計算すれば同じ土俵に立つことになるか。
time ao 2^42 1438 1349 g(3743559068799)=1550 21246.980u 21251.740s 13:35:25.62 86.8% 0+0k 0+0io 120pf+0w
大昔に数年かかったという４兆は半日で終了。
追記
3^n で割った余りでもフィルタできそうだ。9m+4 とか。18で割って分類すると、今6個だけやってるのが5個になる。ただしnを1増やしてもカットできる数は最大1しか増えない。
追記
ここにいろいろ高速化について書いてある。だいたい書いてあることは既にやったな。
http://www.ericr.nl/wondrous/techpage.html
g(k)の計算は他の量に比べて数桁遅いんだそうで。
