0

Computing floor(log_2(x)) can be done by counting the number of zeros for which there are many fast algorithms.

Are there any similar tricks for computing floor(log_{2^(1/4)}(x)) when x is an 64 bit unsigned int ranging over all possible values?

Since log_{2^(1/4)}(x)=4*log_2(x)=log_2(x^4), this would be equivalent to finding an efficient algorithm for floor(4*log_2(x)) or floor(log(x^4)).

Community
  • 1
  • 1
Frumple
  • 3
  • 1

1 Answers1

0

After we compute floor(log_2(x)), we can divide z = x / 2^floor(log_2(x)) and consider the problem of computing floor(log_{2^(1/4)}(z)) where z belongs to the range [1, 2), since floor(log_{2^(1/4)}(x)) = 4 floor(log_2(x)) + floor(log_{2^(1/4)}(z)). There are only four possibilities for floor(log_{2^(1/4)}(z)), so two comparisons (three for branchless) of z with the constants

(2^(1/4))^1, (2^(1/4))^2, (2^(1/4))^3

will suffice. To avoid floating-point arithmetic entirely, the "divide" can be implemented as a left shift that leaves (2^63) z to be compared with similarly represented constants.

Now with C code:

#include <assert.h>
#include <stdio.h>

static int mylog(unsigned long long x) {
  assert(x > 0ULL);
  /* compute n = floor(log2(x)) */
  unsigned long long y = x | (x >> 1);
  y |= y >>  2;
  y |= y >>  4;
  y |= y >>  8;
  y |= y >> 16;
  y |= y >> 32;
  y -= (y >> 1) & 0x5555555555555555ULL;
  y = (y & 0x3333333333333333ULL) + ((y >> 2) & 0x3333333333333333ULL);
  y = (y + (y >>  4)) & 0x0f0f0f0f0f0f0f0fULL;
  y = (y + (y >>  8)) & 0x00ff00ff00ff00ffULL;
  y = (y + (y >> 16)) & 0x0000ffff0000ffffULL;
  y = (y + (y >> 32)) & 0x00000000ffffffffULL;
  int n = (int)y - 1;
  /* normalize x and use binary search to find the last two bits of the log */
  x <<= 63 - n;
  n <<= 2;
  if (x < 0xb504f333f9de6485ULL) {
    return x < 0x9837f0518db8a970ULL ? n     : n + 1;
  } else {
    return x < 0xd744fccad69d6af5ULL ? n + 2 : n + 3;
  }
}

int main(void) {
  unsigned long long x;
  while (scanf("%llu", &x) == 1) {
    printf("%d\n", mylog(x));
  }
}

The Mathematica/Wolfram Alpha query to compute the constants is

BaseForm[Table[Ceiling[2^(63+i/4)],{i,1,3}],16] .
David Eisenstat
  • 64,237
  • 7
  • 60
  • 120
  • That is magical, thanks! As an added bonus, this is actually _more_ accurate than the naive `floor(4*log_2(x))` method due to rounding errors for large `x`. Performance wise, a straightforward translation into Go doesn't produce as big a performance difference as I expected. For random, uniformly distributed uint64 values, I am benchmarking `54.7ns/op` for the naive approach and `19.8ns/op` for the above. The `log_2` calculation portion alone benchmarks at `7.9ns/op`. – Frumple Aug 29 '14 at 13:32