3n+1 問題:まだまだ高速化

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

記号

  • g(k):k からスタートして1に行くステップ数
  • G(n): \Large \frac{max}{\small 1\leq k \leq 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 あたりが最適っぽい。テーブルがキャッシュに収まるからか。

定理2

8で割って4余る数と、その次の数は4ステップ後に同じ数になる。
これを使うと計算時間を7/8 に出来る(微妙)。
追記:TBしていただいたこちらの記事:http://www.typemiss.net/blog/kounoike/20060716-87
いやすごい。1/3くらいに減らせるのか。

定理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 を持つのでダメ。

高速化4

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