コラッツ予想 (Collatz conjecture) その二

前回（↓）の続き。
nazodane.hatenadiary.org

コラッツ問題の計算で、最後の手前の数字は「3の倍数かつ1を足すと2の累乗数になる数」になるはず（偶数処理を遅延した場合はそれに2の累乗数を掛ける必要あるけど、今回は遅延しないで考えることとする）。

1を足すと2の累乗数になる数というのは二進法で1が連続してる数のこと（e.g. 0b1、0b11、0b111、...）であり、
その中で3の倍数になるのは1の連続が偶数個の場合となっている (以下の検証はPython)。

>>> print [(2**x-1) for x in range(1, 32)]
[1, 3, 7, 15, 31, 63, 127, 255, 511, 1023, 2047, 4095, 8191, 16383, 32767, 65535, 131071, 262143, 524287, 1048575, 2097151, 4194303, 8388607, 16777215, 33554431, 67108863, 134217727, 268435455, 536870911, 1073741823, 2147483647]
>>> print [(2**x-1)%3 for x in range(1, 32)] # 0の時に割り切れる
[1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1, 0, 1]
>>> print [(2**(x*2)-1)%3 for x in range(1, 32)]
[0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]

3を掛ける前はどうなるかというと以下のようになる。

>>> print [(2**(x*2)-1)/3 for x in range(1, 32)]
[1, 5, 21, 85, 341, 1365, 5461, 21845, 87381, 349525, 1398101, 5592405, 22369621, 89478485, 357913941, 1431655765, 5726623061, 22906492245, 91625968981, 366503875925, 1466015503701, 5864062014805, 23456248059221, 93824992236885, 375299968947541, 1501199875790165, 6004799503160661, 24019198012642645, 96076792050570581, 384307168202282325, 1537228672809129301]
>>> print [bin((2**(x*2)-1)/3) for x in range(1, 32)]
['0b1', '0b101', '0b10101', '0b1010101', '0b101010101', '0b10101010101', '0b1010101010101', '0b101010101010101', '0b10101010101010101', '0b1010101010101010101', '0b101010101010101010101', '0b10101010101010101010101', '0b1010101010101010101010101', '0b101010101010101010101010101', '0b10101010101010101010101010101', '0b1010101010101010101010101010101', '0b101010101010101010101010101010101', '0b10101010101010101010101010101010101', '0b1010101010101010101010101010101010101', '0b101010101010101010101010101010101010101', '0b10101010101010101010101010101010101010101', '0b1010101010101010101010101010101010101010101', '0b101010101010101010101010101010101010101010101', '0b10101010101010101010101010101010101010101010101', '0b1010101010101010101010101010101010101010101010101', '0b101010101010101010101010101010101010101010101010101', '0b10101010101010101010101010101010101010101010101010101', '0b1010101010101010101010101010101010101010101010101010101', '0b101010101010101010101010101010101010101010101010101010101', '0b10101010101010101010101010101010101010101010101010101010101', '0b1010101010101010101010101010101010101010101010101010101010101']

ここから2**nを掛けて（右側に0を増やす）、1を引いて（LSBをクリアしてその下のビットを全て立てる）、3の倍数を作っては3で割るというのを繰り返していけば良いはず。
前二つの操作は合わせると0b1をクリアして右側に1を増やしていくということになる。

3で割り切れるかどうかは以下みたいに変形できるらしい（gcc -O3の結果より）。

#include <stdint.h>
uint64_t is_divable3(uint64_t x){
  return (x % 3 == 0);
}

uint64_t is_divable3_opt(uint64_t x){
  return 0x5555555555555555 >= 0xaaaaaaaaaaaaaaab * x;
}

上記を全て合わせると逆から計算するこんなコードが出来上がったけど上手く動いているかは未確認。

#include <stdint.h>
#include <stdio.h>
//#include <assert.h>
void loop(uint64_t x, int l){
//  assert(x&1);
  x = 0xaaaaaaaaaaaaaaab * (x ^ 1);
  for (int i = 1; i < 64; i ++){
    x <<= 1; x += 0xaaaaaaaaaaaaaaab;
    if (0x5555555555555555 >= x /*暫定→*/ && x < 0xffff){
      printf("l=%d::%0lx\n", l, x);
      if (l) loop(x, l-1);
    }
  }
}
int main(void){
  uint64_t s = 1;
  printf("l=ss::1\n");
  for(int i = 1; i < 16; i++){
    s <<= 2; s++;
    printf("l=s::%0lx\n", s);
    loop(s, 99);
  }
  return 0;
}

偶数のケースと1のケースを修正。

同じ値でループするケースを削除。

重複無くしたはず。