$ time ./collatz_tab 10000000000
check table: 10226 / 32768 (ratio:0.312073)
h(10000000000) = 9780657630 where g(9780657630:246f8fdde) = 1133
real 5m59.225s
user 5m58.718s
sys 0m0.140s
この結果を出したコードはこちらに。このエントリでは「一捻り」の前までです。
そうそう。「世界記録は2^58」とか色んな話を読みますが,Collatz予想の検証には今回のお題の「f^k(n)=1となるk」を求める場合と「f^k(n)
さて。まずはDanさんのコードの効率化から。ポイントは2点。「キャッシュミスをなくす」ことと,「剰余テーブルを使って数ステップまとめて飛ばす」ことです。
一つ目の改良だけでも結構効きます。コードではlookupのマクロの辺りだけの話です。
g(i)を計算する過程でcache[j]の値が0だった場合(j
cachesize
これに対して改良版では,このときにg(j)の計算を先に行います。その最後でcache[j]に値が入るので,後になってg(j)を求めるときにはキャッシュから見るのですぐに終わります。
この改良で,g(i)のループを下から回すか上から回すかが関係なくなると思います。
二つ目。n=2^k*i+jという形で書くと,jの値によって,数ステップの操作をまとめて考えることができます。k=3の場合についてみてみましょう。
8i + 0 -> 4i + 0 -> 2i + 0 -> i + 0
8i + 1 -> 24i + 4 -> 12i + 2 -> 6i + 1 -> 18i + 4 -> 9i + 2
8i + 2 -> 4i + 1 -> 12i + 4 -> 6i + 2 -> 3i + 1
8i + 3 -> 24i + 10 -> 12i + 5 -> 36i + 16 -> 18i + 8 -> 9i + 4
8i + 4 -> 4i + 2 -> 2i + 1 -> 6i + 4 -> 3i + 2
8i + 5 -> 24i + 16 -> 12i + 8 -> 6i + 4 -> 3i + 2
8i + 6 -> 4i + 3 -> 12i + 10 -> 6i + 5 -> 18i + 16 -> 9i + 8
8i + 7 -> 24i + 22 -> 12i + 11 -> 36i + 34 -> 18i + 17 -> 54i + 52 -> 27i + 26
iの係数が偶数である限り,定数項だけで全体の偶奇が決まります。
i>0ならばこのテーブルを参照することで一気に進めます。
これを実装したのが以下です。itaさんとはキャッシュミス時の対応が違いますね。こちらは自己再帰でキャッシュを計算させますが,itaさんはミスになってる箇所をマークしておき,計算後に遡って値を埋めていく方法のようです。
#include <stdio.h>
#include <stdlib.h>
#define TAB_BITS 12
#define TAB_NUM (1<<TAB_BITS)
#ifndef CACHE_SIZE
#define CACHE_SIZE 512*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 lookup_or_call(n) ((n) < cachesize ? (cache[(n)]?cache[(n)]:g(n)) : 0)
#define store(n,g) if ((n) < cachesize && g <= MAX_U16) cache[(n)] = (g)
typedef struct {
int coef;
int shift;
int step;
int last;
} table_entry;
table_entry table[TAB_NUM];
void dump_table(){
int i;
for(i = 0; i < TAB_NUM; ++i){
table_entry ent = table[i];
printf("%stab[%d]: %d * x + %d (%d step) last:%d\n",(ent.coef > TAB_NUM ? "* " : " "),
i, ent.coef, ent.shift, ent.step, ent.last);
}
}
void create_skip_table(int i){
table_entry ent, ent_min;
ent.coef = TAB_NUM;
ent.shift = i;
ent.step = 0;
ent.last = 0;
ent_min.coef = TAB_NUM*3+1;
while(!(ent.coef&1)){
++ent.step;
if(ent.shift&1){
ent.coef *= 3;
ent.shift *= 3;
++ent.shift;
}else{
ent.coef /= 2;
ent.shift /= 2;
}
if(ent_min.coef > ent.coef) ent_min = ent;
}
/*table[i] = ent_min;*/
table[i] = ent;
}
void create_last_table(){
int i;
for(i=1; i<TAB_NUM;++i){
if(table[i].last > 0) continue;
int n = i;
int step = 1;
while(n > 1){
n = (n&1) ? 3*n+1 : n>>1;
if(n<TAB_NUM && table[n].last > 0){
step += table[n].last;
break;
}
++step;
}
table[i].last = step;
}
}
void create_table(){
int i;
for(i = 0; i < TAB_NUM; ++i){
create_skip_table(i);
}
create_last_table();
}
I64 g(I64 n){
I64 result = lookup(n);
if (result) return result;
I64 i = n;
I64 j;
I64 l;
result = 0;
int gn = 0;
while (j = i>>TAB_BITS){
int k = i & (TAB_NUM-1);
table_entry ent = table[k];
gn += ent.step;
i = ent.coef * j + ent.shift;
l = lookup_or_call(i);
if (l){
result = gn + l;
break;
}
}
if (i < 1) fprintf(stderr, "overflow! g(%lld) = %lld\n", n, i);
if(!result) result = gn + table[i&(TAB_NUM-1)].last;
store(n, result);
return result;
}
I64 *h(I64 n){
I64 i, nmax = 0, gmax = 0, gnext = 0;
static I64 result[2];
for (i = n>>1; i <= n; i++){
gnext = g(i);
if (gnext > gmax){
nmax = i; gmax = gnext;
}
if (i & 0xffff) continue;
printf("h(%lld) = %qd where g(%lld:%llX) = %lld\r", i, nmax, nmax, nmax, gmax);
fflush(stdout);
}
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 = atoll(argv[1]);
}
init_cache(n);
create_table();
I64 *hg = h(n);
printf("h(%lld) = %lld where g(%lld:%llx) = %lld\n", n, hg[0], hg[0], hg[0],hg[1]);
return 0;
}
とりあえずここまで。続きはまた後で。(バグが見つかったため)
Trackback URL for this post:
http://old.typemiss.net/trackback/86
