■[理系ネタ][J] 3n+1 問題：まだまだ高速化

http://ll.jus.or.jp/2006/blog/doukaku2

記号

• g(k):k からスタートして1に行くステップ数
• G(n): g(k)
• h(n): g(k)=G(n)となる k

定理：2 h(n) ＞ n

証明：2 h(n) ≦ n とすると、2 h(n) から 1+G(n)ステップで1に行くので、g(2h(n)) ＞G(n) で最大値の定義に反する。

高速化 1

n/2 から n までだけ調べれば良い。

高速化 2

k からスタートして i ステップ目で k1 (k1＜n) になったとする。h(k1) + i ＜ h(k-1) なら記録更新の望みがないので計算を打ち切って次の k に行っても良い。これを使う場合、1からnへ昇順でいく必要がある。これは思ったほど効かなかったんでやめ。

高速化 3

k を 2^a で割った余りで場合分けして、何ステップか一気に進む。a=10 あたりが最適っぽい。テーブルがキャッシュに収まるからか。

定理２

８で割って４余る数と、その次の数は４ステップ後に同じ数になる。

これを使うと計算時間を7/8 に出来る（微妙）。

追記：TBしていただいたこちらの記事:http://www.typemiss.net/blog/kounoike/20060716-87

いやすごい。1/3くらいに減らせるのか。

定理３

h(n) は 3で割って2余る数ではない。

証明

k=3m+2 とする。2m+1 (＜3m+2)からスタートすると 6m+4→3m+2 と到達するので g(2m+1)=g(3m+2)+2。自分より小さいkが大きい g を持つのでダメ。

高速化４

3で割って2余る数は飛ばす。時間が2/3になると思ったがあまり効かない。なぜだ。

7/17　追記：6で割って余りが2,4,5 だと飛ばす、てのもできる。for ループを手動アンロールして %6 計算と if 文を無くしたら速くなった。

結果

以上の改造パーツを #ifdef で囲ってゴテゴテつけたコードでやってみた。

• 10^7: g(8400511) = 686 (0.470u 0.490s 0:01.04 92.3%)
• 10^8: g(63728127) = 950 (4.210u 4.310s 0:08.73 97.5%)
• 10^9: g(670617279) = 987 (45.980u 46.140s 1:32.47 99.6%)
• 10^10: g(9780657630) = 1133 (553.950u 554.270s 18:31.78 99.6%)
• 10^11 (これは64bit 越えて間違ってるかも): g(75128138247) = 1229 (22125.850u 0.360s 6:09:25.00 99.8%)

10^10 も10分以内。計算時間は T(n) = 2.5e-9*n*log(n) とフィットできる。このページにある四兆を数年、てのと比較すると T(四兆)=3.4 日。とりあえずムーアの法則に感謝ってところ。

（追記）10^11 以降は64bit でも足りなくなるらしい。多倍長だと一気に重くなりそう。

世界記録は n=2^58 らしいけど、これを代入すると T = 三万年とか出た。まだまだ高速化が甘いらしい。

7/24 追記

64bit 溢れそうになったら96bit 多倍長に移行するコード＋6で割った余りでスキップ＋2^16で割った余りでスキップ、でやってみた。

time a.out 10000000000
114.820u 115.040s 3:50.25 99.8% 0+0k 0+0io 101pf+0w

10の10乗で2分を切った。ここまでは桁溢れなし。

time a.out 100000000000
g(75128138247)=1229
1305.420u 1305.860s 43:37.88 99.7%      0+0k 0+0io 101pf+0w

10の11乗は 64bit over 頻出。しかし全体に対する割合はまだ小さい。次は 2^37 あたりから倍々で行くか。

2^37: g(133561134663)=1235
1970.970u 1971.630s 1:06:32.99 98.7%    0+0k 0+0io 95pf+0w

2^38: g(202485402111)=1308
4726.010u 4727.130s 2:47:21.74 94.1%    0+0k 0+0io 112pf+0w

2^41: g(2081751768559)=1438
44328.780u 44335.610s 29:48:31.58 82.6% 0+0k 0+0io 114pf+0w

• ソース、Ver 0.3

• ソース、Ver 0.5 112bit
• スキップテーブル初期化
• 飛びテーブル初期化

Permalink | コメント(0) | トラックバック(0)