h(4294967295) = 2610744987 where g = 1051

itaさんの追記を読んでいたら、私も行き着くところまで行きたくなった。

キミならどう書く 2.0 - ROUND 2 - ? Lightweight Language Ring

このとき
h(n) = k, 1 ≦ k ≦ n ∧ g(k) = max (g(1),g(2),…,g(n))
について h(100) を求めよ．

JGeek Log - 3n+1 問題

time a.out 1000000000 (テーブルサイズ = n * 0.03)
g(670617279) = 987

とはいっても、行ったのはunsigned 32bit intの限界までなのだけど。

以下、FreeBSD 6.1, Dual Xeon 2.8MHz による結果。

% /usr/bin/time -l ./a.out
h(4294967295) = 2610744987 where g(2610744987) = 1051
     2629.46 real      2611.42 user         2.90 sys
    131892  maximum resident set size
         4  average shared memory size
      7574  average unshared data size
       128  average unshared stack size
     32860  page reclaims
         0  page faults
         0  swaps
         0  block input operations
         0  block output operations
         0  messages sent
         0  messages received
         0  signals received
        54  voluntary context switches
    215379  involuntary context switches

見ての通り、一時間かかっていない。また、n = 10億までの結果はitaさんのものと一致する。

% /usr/bin/time ./a.out 1000000
h(1000000) = 837799 where g(837799) = 525
        0.16 real         0.16 user         0.00 sys
% /usr/bin/time ./a.out 10000000
h(10000000) = 8400511 where g(8400511) = 686
        1.77 real         1.73 user         0.03 sys
% /usr/bin/time ./a.out 100000000
h(100000000) = 63728127 where g(63728127) = 950
       16.35 real        16.06 user         0.23 sys
% /usr/bin/time ./a.out 1000000000
h(1000000000) = 670617279 where g(670617279) = 987
      357.63 real       355.45 user         0.24 sys

itaさんはもう少し先に行ったみたいだけど、「渡し」も「自身」がなかったのでとりあえず無難な線で一区切り。

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

実はMacBook Proはもっと優秀で、n=1億でこんな結果になった。

% /usr/bin/time ./a.out 100000000
h(100000000) = 63728127 where g(63728127) = 950
       14.67 real        14.29 user         0.28 sys

しかもjkondoのポッドキャストを聞きながら、である。Intel Core Duoって優秀かも。

コードは以下のとおり。見ての通り、cacheはunsigned shortで実装している。それで充分であることはg(n)を観察すればわかるので。その結果、200MB弱のメモリーでitaさんよりよい成績を上げている。あと、実行時間の長いプログラムのための工夫として、65565回ごとに途中結果を表示するようにしている。この部分をはしょればもう少し早くなるかも知れない。

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

#ifndef CACHE_SIZE
#define CACHE_SIZE 64*1024*1024 /* large enough */
#endif

#define MAX_U16        65535
typedef unsigned short U16;
typedef long long      I64;

U16 *cache;
int cachesize;
#define lookup(n)  ((n) < cachesize ? cache[(n)] : 0)
#define store(n,g) if ((n) < cachesize && g <= MAX_U16) cache[(n)] = (g)

I64 g(I64 n){
    I64 result = lookup(n);
    if (result) return result;
    I64 i = n;
    I64 l;
    result = 1;
    while (i > 1){
        i = i & 1 ? i*3+1 : i >> 1;
        l = lookup(i);
        if (l){
            result += l;
            break;
        }else{
            result++;
        }
    }
    if (i < 1)       fprintf(stderr, "overflow! g(%qd) = %qd\n", n, i);
    store(n, result);
    return result;
}

I64 *h(I64 n){
    I64 i, nmax = 0, gmax = 0, gnext = 0;
    static I64 result[2];
    for (i = 1; i <= n; i++){
        gnext = g(i);
#ifdef VERBOSE
        fprintf(stderr, "g(%qd) = %qd\r", i, gnext);
#endif
        if (gnext > gmax){
            nmax = i; gmax = gnext;
        }
        if (i & 0xffff) continue;
        printf("h(%qd) = %qd where g(%qd) = %qd\r", i, nmax, nmax, gmax);
    }
    result[0] = nmax; result[1] = gmax;
    return result;
}

#define min(x,y) ((x) < (y) ? (x) : (y))

void init_cache(I64 n){
    cachesize = min(n, CACHE_SIZE);
    cache     = (U16 *)calloc(cachesize, sizeof(U16));
    if (cache == NULL) exit(-1);
}

int main(int argc, char **argv){
    I64 n = 0xffffffff; /* look @ this! */
    if (argc > 1){
        n = atoi(argv[1]);
    }
    init_cache(n);
    I64 *hg = h(n);
    printf("h(%qd) = %qd where g(%qd) = %qd\n", n, hg[0], hg[0], hg[1]);
    return 0;
}

Enjoy!

Dan the Number Cruncher

追記:

逆まわし、すなわち最大の数からチェックした方がcacheがより有効に使われることに気がついた。

/usr/bin/time -l ./a.out 100000000
h(100000000) = 63728127 where g(63728127) = 950
      184.75 real       183.95 user         0.26 sys

上に以下のpatch。

--- gogo.c.orig 2006-07-14 13:33:49.000000000 +0900
+++ gogo.c      2006-07-14 15:12:55.000000000 +0900
@@ -38,7 +38,7 @@
 I64 *h(I64 n){
     I64 i, nmax = 0, gmax = 0, gnext = 0;
     static I64 result[2];
-    for (i = 1; i <= n; i++){
+    for (i = n; i > 1 ; i--){
        gnext = g(i);
 #ifdef VERBOSE
        fprintf(stderr, "g(%qd) = %qd\r", i, gnext);
@@ -47,7 +47,8 @@
            nmax = i; gmax = gnext;
        }
        if (i & 0xffff) continue;
-       printf("h(%qd) = %qd where g(%qd) = %qd\r", i, nmax, nmax, gmax);
+       printf("nmax = %qd, gmax = %qd, n = %qd\r", nmax, gmax, i);
+       fflush(stdout);
     }
     result[0] = nmax; result[1] = gmax;
     return result;
@@ -64,7 +65,7 @@
 int main(int argc, char **argv){
     I64 n = 0xffffffff; /* look @ this! */
     if (argc > 1){
-       n = atoi(argv[1]);
+       n = atoll(argv[1]);
     }
     init_cache(n);
     I64 *hg = h(n);