どこまで行けるか? Collatz予想

Danさんのコードを元に改良をしていたのですが,itaさんがほとんど同じことをされていたので,もう一捻りしてみました。 その結果10^10がcygwin on Pentium4 2.6GHz上で6分程度で計算できるくらいには高速化しました。 $ 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