ツリー状の自然数の列(2)

前回、64ビット整数(__int64)を使ったが、


n /= 2;

ここを、


n >>= 1;

としたら、全体で4倍速くなった。
アセンブリ見てみるもんだね。
しかし、符号なし32ビット整数を2つ使ったら、さらに1.8倍になった。
こんな感じ。


// n = n * 3 + 1
void triple(uint n[2]) {
n[0] = n[0] * 3 + 1;
n[1] = n[1] * 3;
uint over = n[0] >> 30;
n[1] += over;
n[0] &= (1 << 30) - 1;
}

// n /= 2
void half(uint n[2]) {
n[0] /= 2;
uint bit = n[1] & 1;
n[1] /= 2;
if(bit)
n[0] += 1 << 29;
}

上の2ビットを桁あふれ用に使っているので、
ちょっと表現できる範囲は狭くなっている。


速くなったので、前回より1桁データを増やしてみた。

やはりなんとも言えない。


あまりうまくいかないので、
今度は1に到達するまでの平均回数を調べた。
これは、1からその数まででの平均なので、
その数付近での平均は、その数を掛けて微分すればよい。
結果的に、近似曲線のLnにかかっている係数を足せばよい。

これだとかなりよく直線に乗る。
数値にすると、大きなところはかなり直線に乗る。


1048576 131.893 7.180
2097152 139.064 7.171
4194304 146.248 7.184
8388608 153.433 7.185
16777216 160.634 7.201
33554432 167.847 7.213
67108864 175.062 7.215

この7.215という幅(もうちょっと大きくなりそう)は、
あまり根拠はないが、次のように考えられる。
1から逆に辿ったとき、2倍していくと辿れるが、
3の剰余が1の場合だけ、枝分かれがある。
すなわち、1に到達するまでの回数(距離と表現しよう)が1増すにつれて、
その距離になるの数の個数は4/3倍になる、と思われる。
とすると、個数が2倍になるには、
log4/32 = 2.409距離が増す必要がある。
2.409 * 3 = 7.227 ≒ 7.215だから、
実際にはこの3倍かかっている、ということになる。
ちゃんとパックしたときより大きな数がこの距離内に入っている、
ということになるのではないだろうか。


(追記)
もう少しやってみた。


134217728 182.287 7.225
268435456 189.515 7.228
536870912 196.741 7.226



#include

using namespace std;

typedef unsigned int uint;
int calc_length(int m);
void triple(uint n[2]);
void half(uint n[2]);

void main() {
const int N = (int)7e7;

__int64 sum = 0;
int ptPrint = 1;
for(int i = 1; i <= N; i++) {
int length = calc_length(i);
sum += length;
if(i == ptPrint) {
cout << i << " " << sum / (double)i << endl;
cerr << i << " " << sum / (double)i << endl;
ptPrint *= 2;
}
}
}

int calc_length(int m) {
const uint limit = 0xFFFFFFFE / 3;
uint n[2];
n[0] = m; n[1] = 0;
int len;
for(len = 0; n[0] != 1 || n[1] != 0; len++) {
if(n[0] & 1) {
if(n[1] > limit) {
cerr << "overflow." << endl;
exit(1);
}
triple(n);
}
else {
half(n);
}
}
return len;
}

void triple(uint n[2]) {
n[0] = n[0] * 3 + 1;
n[1] = n[1] * 3;
uint over = n[0] >> 30;
n[1] += over;
n[0] &= (1 << 30) - 1;
}

void half(uint n[2]) {
n[0] /= 2;
uint bit = n[1] & 1;
n[1] /= 2;
if(bit)
n[0] += 1 << 29;
}