64-bit unsigned integer cube root

gistfile1.cpp

#include <stdio.h>

typedef unsigned int U32;
typedef unsigned __int64 U64;

// ---- actual cube root code

U32 icbrt64(U64 x) {
  int s;
  U32 y;
  U64 b;

  y = 0;
  for (s = 63; s >= 0; s -= 3) {
    y += y;
    b = 3*y*((U64) y + 1) + 1;
    if ((x >> s) >= b) {
      x -= b << s;
      y++;
    }
  }
  return y;
}

// ---- test driver

int main()
{
  U32 max = 2642245; // floor(cbrt(2^64))

  for (U32 i=1; i<=max; ++i) { // 
    U64 v = (U64) i*i*i;
    U64 cb = icbrt64(v);
    
    if (cb != i) {
      printf("err1! i=%d cb=%d\n", i, (int) cb);
      break;
    }
    
    cb = icbrt64(v-1);
    if (cb != i-1) {
      printf("err2! i-1=%d cb=%d\n", i-1, cb);
      break;
    }
  }
  
  // check max value
  U64 v = ~0ull;
  if (icbrt64(v) != max) {
    printf("err3!\n");
  }
  
  return 0;
}