3n+1 問題

http://ll.jus.or.jp/2006/blog/doukaku2
あ、これ去年某大学の並列計算の講義でネタにしました。簡単な整数演算を大量にやって何か面白いことが出来る話の一つとして、2CH のトリップや/etc/passwd のクラックなどと共に紹介しました。
参考リンク：http://www.math.grinnell.edu/~chamberl/conference/proceedings.html
高速化の方法として考えられるのは、計算済みの g(n) を保存して、後の計算でそれにヒットしたら計算を省略する、てのだけど、必要なメモリが予測できないのでトリッキーかな。

とりあえず書いてみた

実行結果。AMD Sempron のLINUX + GCC。10000 までは速過ぎて測定不能、100000で0.01sec、1000000 でint のオーバーフローで落ちる。使用メモリは数メガ。

追記

一定の値を越えると 64 bit 演算にスイッチするよう改造してみた。gprof で調べると 64 bit 演算に費している時間は n=10^8 の場合でも全体の3% 程度なんで、全部64bit でやるよりも速くなっていると思われる。 追記：全部64bitでやったほうが速かった。
n=10^8 あたりからメモリ不足でテーブル引きができなくなって遅くなる。

time a.out 100
g(97) = 119
0.000u 0.000s 0:00.00 0.0%      0+0k 0+0io 84pf+0w

time a.out 1000
g(871) = 179
0.000u 0.000s 0:00.00 0.0%      0+0k 0+0io 84pf+0w

time a.out 10000
g(6171) = 262
0.000u 0.000s 0:00.00 0.0%      0+0k 0+0io 84pf+0w

time a.out 100000
g(77031) = 351
0.010u 0.010s 0:00.02 100.0%    0+0k 0+0io 84pf+0w

time a.out 1000000
セグメントエラー

64 bit 版 (バグがあった。以下修正済み)

time a.out 1000000
g(837799) = 525
0.110u 0.150s 0:00.27 96.2%     0+0k 0+0io 90pf+0w

time a.out 10000000
g(8400511) = 686
1.000u 1.080s 0:02.08 100.0%    0+0k 0+0io 97pf+0w

time a.out 100000000 (テーブルサイズ = n * 0.1)
g(63728127) = 950
26.410u 26.470s 0:52.98 99.8%   0+0k 0+0io 101pf+0w

time a.out 1000000000 (テーブルサイズ = n * 0.03)
g(670617279) = 987
434.710u 435.190s 14:31.53 99.8%        0+0k 0+0io 97pf+0w

これは自身なし
n=10000000000
g(9780657630) = 1133
4222.760u 4223.620s 2:23:12.06 98.3%    0+0k 0+0io 95pf+0w

10^10 以降は n が int 上限を越えるので全部 64 bit でやったほうがいいかも。あとテーブルは役に立たなくなるんで素直に全部計算したほうが一番速いかも。たぶん n log(n) の計算時間になる。
n=4000000 の時、テーブルのサイズをいろいろ変えて測ると以下のようになった。テーブルサイズがnと同じ時が最適っぽい。

サイズ	CPU時間
0.1 n	0.960
0.2 n	0.690
0.4 n	0.480
0.6 n	0.390
0.8 n	0.350
1.0 n	0.340
1.2 n	0.380
1.4 n	0.390
1.6 n	0.440

ソースはこちら。64 bit で行く部分は #ifdef GO_64 で囲ってある。

#include <stdio.h>
#include <stdlib.h>

#define GO_64 0x40000000

#ifdef GO_64
void handle_64(unsigned int *x, int *g)
{
    int g0 = 0;
    unsigned long long x64 = *x;

    while (x64 >= GO_64)
    {
	unsigned long long m = x64>>1;
	if (x64 & 1)
	{
	    x64 = 3*m+2;
	    g0 += 2;
	}
	else
	{
	    x64 = m;
	    g0++;
	}

	if (x64 > 0x4000000000000000LL)
	{
	    printf("Error 64 bit is not enough!\n");
	    exit(1);
	}
    }
    *x = (unsigned int) x64;
    (*g) += g0;
}
#endif

int  main(int argc, char **argv)
{
    int i,j,k;
    #define H_BUFSIZE (1<<16)
    int x_log[H_BUFSIZE];
    int g_log[H_BUFSIZE];

    int gmax = 1;
    int kmax = 1;

    int n = atoi(argv[1]);
    int *gtable;
    int gtable_size = n;

    if (argc==3)
    {
        double r = (double)atof(argv[2]);
	gtable_size = (int)(n * r + 1);
    }

    #define GTABLE_MAX 30000000
    if (gtable_size > GTABLE_MAX) gtable_size = GTABLE_MAX;
    gtable = (int *)malloc(gtable_size * sizeof(int));
    if (gtable == NULL)
    {
	printf("Not enough memory\n");
	exit(1);
    }
    printf("Buffer size ratio %f\n",1.0*gtable_size/n);

    for(i=0;i<gtable_size;i++)
	gtable[i] = 0;

    gtable[0] = 1;

    for(i=0;i<n;i++)
    {
	int g=0;
	int n_log = 0;
	unsigned int x = i+1;

        while( (x > gtable_size) || (gtable[x-1] == 0))
	{

	    if (x <= gtable_size)
	    {
		x_log[n_log] = x;
		g_log[n_log] = g;
		n_log++;
	    }

            if (x & 1) x = x*3+1;
	    else x >>=1;
	    g++;

#ifdef GO_64
	    if (x > GO_64) handle_64(&x, &g);
#endif

	}

	g += gtable[x-1];

	for(j=0;j<n_log;j++)
	{
	    int g0 = g_log[j];
	    x = x_log[j];
	    gtable[x-1] = g-g0;
	}

	if (g>gmax)
	{
	    gmax = g;
	    kmax = i+1;
	}

    }//i

    printf("g(%d) = %d\n",kmax, gmax);
}