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
ここで，8i+4と8i+5に注目です。この2つの値はどちらも4回 (経過を見ると3回) の操作で同じ値に合流しています。つまり，任意のi>0について，g(8i+4)=g(8i+5)が成り立ちますので，h(n)の計算の過程で8i+5は計算する必要がありません。こんなパターンはk（2^k*i+jの形でテーブルを作るときのk）が大きくなるとより多く出現します。例えばk=16なら65536個のjのうち，19275個を調べれば十分です。この処理を入れたのが以下のコードと前のエントリに載せた結果です。Danさんのコードは分かりやすくていじりやすくて良いですね。
#include <stdio.h>
#include <stdlib.h>
#define TAB_BITS 14
#define TAB_NUM (1<<TAB_BITS)
#define SKIP_TABLE_SIZE 16000 /* 注意:TAB_NUM/3程度が必要 */
#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 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;
} skip_table_entry;
int last_table[TAB_NUM];
skip_table_entry* skip_table[TAB_NUM];
char need_check_table[TAB_NUM];
skip_table_entry skip_table_pool[SKIP_TABLE_SIZE];
int skip_table_used = 0;
void dump_table(){
int i;
for(i = 0; i < TAB_NUM; ++i){
skip_table_entry ent = *skip_table[i];
printf("%stab[%d]: %d * x + %d (%d step) last:%d\n", (ent.coef > TAB_NUM ? "* " : " "),
i, ent.coef, ent.shift, ent.step, last_table[i]);
}
}
void create_last_table(){
int n;
for(n=0; n<TAB_NUM; ++n) last_table[n] = 0;
for(n=1; n<TAB_NUM;++n){
int i = n;
int step = 1;
while(i > 1){
i = (i&1) ? 3*i+1 : i>>1;
if(i<TAB_NUM && last_table[i] > 0){
step += last_table[i];
break;
}
++step;
}
last_table[n] = step;
}
}
typedef struct list_entry_t {
int value;
skip_table_entry *skip;
struct list_entry_t *next;
} list_entry;
void create_table(){
int i;
int pow3tab[TAB_BITS+1];
list_entry* heads[TAB_BITS+1];
int count = 0;
pow3tab[0] = 1;
heads[0] = NULL;
for(i=1; i<=TAB_BITS;++i){
pow3tab[i] = pow3tab[i-1]*3;
heads[i] = NULL;
}
for(i = 0; i < TAB_NUM; ++i){
int pow2 = TAB_BITS;
int pow3 = 0;
int shift = i;
int step = 0;
while(pow2>0){
++step;
if(shift&1){
shift = shift * 3 + 1;
++pow3;
}else{
shift>>=1;
--pow2;
}
}
list_entry *p = heads[pow3];
if(p == NULL){
p = (list_entry*) malloc(sizeof(list_entry));
p->next = NULL;
p->value = shift;
p->skip = NULL;
heads[pow3] = p;
}else{
while(p->next != NULL && p->next->value != shift) p = p->next;
if(p->next == NULL){
p->next = (list_entry*)malloc(sizeof(list_entry));
p = p->next;
p->next = NULL;
p->skip = NULL;
p->value = shift;
}else{
p = p->next;
}
}
if(p->skip == NULL){
need_check_table[i] = 1;
skip_table[i] = &skip_table_pool[skip_table_used++];
count++;
skip_table[i]->coef = pow3tab[pow3];
skip_table[i]->shift = shift;
skip_table[i]->step = step;
p->skip = skip_table[i];
}else{
need_check_table[i] = 0;
if(p->skip->step != step){
fprintf(stderr,"step mismatch. for i:%d step:%d expected step:%d\n",
i, step, p->skip->step);
exit(-1);
}
skip_table[i] = p->skip;
}
}
printf("need for check : %d / %d (ratio:%f)\n",count, TAB_NUM, 1.0*count/TAB_NUM);
create_last_table();
for(i = 0; i < TAB_BITS; ++i){
list_entry* p = heads[i];
while(p != NULL){
list_entry* p2 = p->next;
free(p);
p = p2;
}
}
}
int 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);
skip_table_entry* ent = skip_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);
return -1;
}
if(!result) result = gn + last_table[i&(TAB_NUM-1)];
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) & ~((long long)TAB_NUM-1) ; i <= n; i++){
if(need_check_table[i&(TAB_NUM-1)] == 0 ) continue;
gnext = g(i);
if (gnext > gmax){
nmax = i; gmax = gnext;
}
if (i & 0xffffff) continue;
printf("h(%llX:%lld) = %lld where g(%lld:%llX) = %lld\r", i,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();
dump_table(); return 0;
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;
}
これで10^10が6分ほど。もうちょっと整理すると速くなりそうですが（スキップテーブルを節約させてみたせいで遅くなっているとか，スキップテーブルと，チェックに使うテーブルを分けちゃえとか），とりあえずこのくらいで。
それと，増えるときが3倍ずつじゃなくなったので，64bitの範囲で収まるかどうかのチェックが0との比較じゃできなくなっているはずです。この辺の改良と多倍長化の工夫が今後の課題といったところでしょうか。とりあえず，5*10^10あたりで運良くチェックに引っかかることがありました。すり抜けてることもあるかもしれません。
Trackback URL for this post:
http://old.typemiss.net/trackback/87
