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
