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

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

前回のコードを適当にベクトライズするとこうなる。

// gcc rcollatz.c -O3 -march=haswell
// clang rcollatz.c -O3 -mavx
#include <stdint.h>
#include <stdio.h>
//#include <assert.h>
void loop(uint64_t x, int l){
//  assert(x&1);
  uint64_t xa[64];
  for (int i = 0; i < 63; i++) xa[i] = x;
  for (int i = 0; i < 63; i++) xa[i] = xa[i] << (i+1); // shift + add vs. fma
  for (int i = 0; i < 63; i++) xa[i] = xa[i] + 0xaaaaaaaaaaaaaaab * ((1 << (i+1)) - 1);
  for (int i = 0; i < 63; i++) {
    if (0x5555555555555555 >= xa[i] /*暫定→*/ && xa[i] < 0xffff){
      printf("l=%d::%0lx\n", l, xa[i]);
      if (l) loop(0xaaaaaaaaaaaaaaab * (xa[i] ^ 1), l-1);
    }
  }
}

int main(void){
  uint64_t s = 1, x = 0;
  printf("l=ss::1\n");
  for(int i = 1; i < 16; i++){
    s <<= 2; s++;
    printf("l=s::%0lx\n", s);

    x <<= 2; x += (0xaaaaaaaaaaaaaaab << 2);

    loop(x, 99);
  }
  return 0;
}

以下のコードはシフトより上と下に分けることが出来る。

bool ge_shift_add(uint32_t x){
  return (CMPVAL >= (x<<SHIFT) + ADDVAL);
}

そのためxのシフトを無くしてこう変形できるはず。

bool ge_shift_add2(uint32_t x){
  return ((((CMPVAL>>SHIFT) - !!((CMPVAL & ((1<<SHIFT)-1)) < (ADDVAL & ((1<<SHIFT)-1))))&(0xffffffff>>SHIFT)) >= ((x + (ADDVAL>>SHIFT))&(0xffffffff>>SHIFT)));
}

ADDVAL>>SHIFT側をpythonで計算するとこんな感じ。

>>> print [(0xaaaaaaaaaaaaaaab * ((1 << (i+1)) - 1) &0xffffffffffffffff) >> (i+1) for i in range(63)]
[6148914691236517205L, 0L, 1537228672809129301L, 0L, 384307168202282325L, 0L, 96076792050570581L, 0L, 24019198012642645L, 0L, 6004799503160661L, 0L, 1501199875790165L, 0L, 375299968947541L, 0L, 93824992236885L, 0L, 23456248059221L, 0L, 5864062014805L, 0L, 1466015503701L, 0L, 366503875925L, 0L, 91625968981L, 0L, 22906492245L, 0L, 5726623061L, 0L, 1431655765L, 0L, 357913941L, 0L, 89478485L, 0L, 22369621L, 0L, 5592405L, 0L, 1398101L, 0L, 349525L, 0L, 87381L, 0L, 21845L, 0L, 5461L, 0L, 1365L, 0L, 341L, 0L, 85L, 0L, 21L, 0L, 5L, 0L, 1L]

偶数が0、奇数が[(4**i)/3 for i in range(1, 33)]の逆順になってる。

((((CMPVAL>>SHIFT) - !!((CMPVAL & ((1<<SHIFT)-1)) < (ADDVAL & ((1<<SHIFT)-1))))&(0xffffffff>>SHIFT)) 側をpythonで計算するとこんな感じ。

>>> print [(lambda shift=i+1, mask=(1<<(i+1))-1:((0x5555555555555555 >> shift) - ((0x5555555555555555&mask) < ((0xaaaaaaaaaaaaaaab*mask) & mask)))&(0xffffffffffffffff>>shift))() for i in range(63)]
[3074457345618258602L, 1537228672809129301L, 768614336404564650L, 384307168202282325L, 192153584101141162L, 96076792050570581L, 48038396025285290L, 24019198012642645L, 12009599006321322L, 6004799503160661L, 3002399751580330L, 1501199875790165L, 750599937895082L, 375299968947541L, 187649984473770L, 93824992236885L, 46912496118442L, 23456248059221L, 11728124029610L, 5864062014805L, 2932031007402L, 1466015503701L, 733007751850L, 366503875925L, 183251937962L, 91625968981L, 45812984490L, 22906492245L, 11453246122L, 5726623061L, 2863311530L, 1431655765L, 715827882L, 357913941L, 178956970L, 89478485L, 44739242L, 22369621L, 11184810L, 5592405L, 2796202L, 1398101L, 699050L, 349525L, 174762L, 87381L, 43690L, 21845L, 10922L, 5461L, 2730L, 1365L, 682L, 341L, 170L, 85L, 42L, 21L, 10L, 5L, 2L, 1L, 0L]

つまり、[(2**(i+2)-3-(-1)**i)/6 for i in range(63)] の逆順になってる。

うーん、このあともシフト計算無くして行ける気がするけど、式が長くなって頭こんがらがっていくので今日はここまで。

別式の追加