JGeek Log

虫記事
[カミキリ星撮表]
科学記事 SF About 物置 Amazon.co.jp
ディアスポラ数理研 /白熱光メモ 読み中
 | 

2006-07-12

[][] 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 SempronLINUX + 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);
}

[][] 3n+1 問題: 並列化

長くなったんで別エントリで。単体CPUでは弾さんのところを見ると最速に達したっぽいんで、複数CPUでの並列化を考えてみる。といっても単純で、k のループを単純に分割する。

    MPI_Init など並列計算オマジナイ
    nprocs = プロセッサ数
    myrank = 自分のプロセッサ番号

    for(k=myrank;k<n;k+=nprocs)
    {
   ...
    }
    通信して各プロセッサが持っている g の最大値のうちで最大のものを表示

小さいkのテーブルの計算が無駄に重複するが、n が大きいと無視できると思われる。各プロセッサの計算時間が不均一になるが、多分±数10%に収まると思う。いま手元に並列機がないんで試せないけど。

 | 
Contact: Mitsuhiro Itakura/板倉充洋
ita.mitsu spam @ gmail, sausage, spam, egg, and com
最近のコメント

#. DATE NAME

1. 05/21 ita
2. 05/21 しげしげ
3. 05/21 ita
4. 05/21 yama
5. 05/21 umajin
6. 04/21 ita
7. 05/14 しげしげ
8. 05/06 ita
9. 05/06 しげしげ
最近のTB

#. DATE  NAME

CALENDAR
<< 2006/07 >>
1
2 3 4 5 6 7 8
9 10 11 12 13 14 15
16 17 18 19 20 21 22
23 24 25 26 27 28 29
30 31