Compute fast log base 2 ceiling

Viewed 38787

What is a fast way to compute the (long int) ceiling(log_2(i)), where the input and output are 64-bit integers? Solutions for signed or unsigned integers are acceptable. I suspect the best way will be a bit-twiddling method similar to those found here, but rather than attempt my own I would like to use something that is already well tested. A general solution will work for all positive values.

For instance, the values for 2,3,4,5,6,7,8 are 1,2,2,3,3,3,3

Edit: So far the best route seems to be to compute the integer/floor log base 2 (the position of the MSB) using any number of fast existing bithacks or register methods, and then to add one if the input is not a power of two. The fast bitwise check for powers of two is (n&(n-1)).

Edit 2: A good source on integer logarithms and leading zeroes methods is Sections 5-3 and 11-4 in Hacker's Delight by Henry S. Warren. This is the most complete treatment I've found.

Edit 3: This technique looks promising: https://stackoverflow.com/a/51351885/365478

Edit 4: C23 is apparently adding stdc_first_(leading/trailing)_(one/zero)

15 Answers

The fastest approach I'm aware of uses a fast log2 that rounds down, combined unconditional adjustment of input value before and after to handle the rounding up case as in lg_down() shown below.

/* base-2 logarithm, rounding down */
static inline uint64_t lg_down(uint64_t x) {
  return 63U - __builtin_clzl(x);
}

/* base-2 logarithm, rounding up */
static inline uint64_t lg_up(uint64_t x) {
  return lg_down(x - 1) + 1;
}

Basically adding 1 to the rounded-down result is already correct for all values except exact powers of two (since in that case the floor and ceil approaches should return the same answer), so it is sufficient to subtract 1 from the input value to handle that case (it doesn't change the answer for the other cases) and add one to the result.

This is usually slightly faster than the approaches that adjust the value by explicitly checking for exact powers of two (e.g., adding a !!(x & (x - 1)) term). It avoids any comparisons and conditional operations or branches, is more likely to simply when inlining, is more amenable to vectorization, etc.

This relies on the "count leading bits" functionality offered by most CPUs using the clang/icc/gcc builtin __builtin_clzl, but other platforms offer something similar (e.g., the BitScanReverse intrinsic in Visual Studio).

Unfortunately, this many return the wrong answer for log(1), since that leads to __builtin_clzl(0) which is undefined behavior based on the gcc documentation. Of course, the general "count leading zeros" function has perfectly defined behavior at zero, but the gcc builtin is defined in this way because prior to the BMI ISA extension on x86, it would have been using the bsr instruction which itself has undefined behavior.

You could work around this if you know you have the lzcnt instruction by using the lzcnt intrinsic directly. Most platforms other than x86 never went through the bsr mistake in the first place, and probably also offer methods to access their "count leading zeros" instruction if they have one.

I will give you the fastest way for x86-64 at the time of writing, and a general technique if you have a fast floor that works for arguments < 2^63, if you care for the full range, then see below.

I am surprised about the low quality of the other answers because they tell you how to get the floor but transform the floor in a very expensive way (using conditionals and everything!) to the ceiling.

If you can get the floor of the logarithm quickly, for example, using __builtin_clzll, then the floor is very easily obtained like this:

unsigned long long log2Floor(unsigned long long x) {
    return 63 - __builtin_clzll(x);
}

unsigned long long log2Ceiling(unsigned long long x) {
    return log2Floor(2*x - 1);
}

It works because it adds 1 to the result unless x is exactly a power of 2.

See the x86-64 assembler difference at the compiler explorer for another implementation of ceiling like this:

auto log2CeilingDumb(unsigned long x) {
    return log2Floor(x) + (!!(x & (x - 1)));
}

Gives:

log2Floor(unsigned long): # @log2Floor(unsigned long)
  bsr rax, rdi
  ret
log2CeilingDumb(unsigned long): # @log2CeilingDumb(unsigned long)
  bsr rax, rdi
  lea rcx, [rdi - 1]
  and rcx, rdi
  cmp rcx, 1
  sbb eax, -1
  ret
log2Ceiling(unsigned long): # @log2Ceiling(unsigned long)
  lea rax, [rdi + rdi]
  add rax, -1
  bsr rax, rax
  ret

For the full range, it is in a previous answer: return log2Floor(x - 1) + 1, this is significantly slower since it uses in x86-64 four instructions compared to three above.

I have benchmarked several implementations of a 64 bit "highest bit". The most "branch free" code is not in fact the fastest.

This is my highest-bit.c source file:

int highest_bit_unrolled(unsigned long long n)
{
  if (n & 0xFFFFFFFF00000000) {
    if (n & 0xFFFF000000000000) {
      if (n & 0xFF00000000000000) {
        if (n & 0xF000000000000000) {
          if (n & 0xC000000000000000)
            return (n & 0x8000000000000000) ? 64 : 63;
          else
            return (n & 0x2000000000000000) ? 62 : 61;
        } else {
          if (n & 0x0C00000000000000)
            return (n & 0x0800000000000000) ? 60 : 59;
          else
            return (n & 0x0200000000000000) ? 58 : 57;
        }
      } else {
        if (n & 0x00F0000000000000) {
          if (n & 0x00C0000000000000)
            return (n & 0x0080000000000000) ? 56 : 55;
          else
            return (n & 0x0020000000000000) ? 54 : 53;
        } else {
          if (n & 0x000C000000000000)
            return (n & 0x0008000000000000) ? 52 : 51;
          else
            return (n & 0x0002000000000000) ? 50 : 49;
        }
      }
    } else {
      if (n & 0x0000FF0000000000) {
        if (n & 0x0000F00000000000) {
          if (n & 0x0000C00000000000)
            return (n & 0x0000800000000000) ? 48 : 47;
          else
            return (n & 0x0000200000000000) ? 46 : 45;
        } else {
          if (n & 0x00000C0000000000)
            return (n & 0x0000080000000000) ? 44 : 43;
          else
            return (n & 0x0000020000000000) ? 42 : 41;
        }
      } else {
        if (n & 0x000000F000000000) {
          if (n & 0x000000C000000000)
            return (n & 0x0000008000000000) ? 40 : 39;
          else
            return (n & 0x0000002000000000) ? 38 : 37;
        } else {
          if (n & 0x0000000C00000000)
            return (n & 0x0000000800000000) ? 36 : 35;
          else
            return (n & 0x0000000200000000) ? 34 : 33;
        }
      }
    }
  } else {
    if (n & 0x00000000FFFF0000) {
      if (n & 0x00000000FF000000) {
        if (n & 0x00000000F0000000) {
          if (n & 0x00000000C0000000)
            return (n & 0x0000000080000000) ? 32 : 31;
          else
            return (n & 0x0000000020000000) ? 30 : 29;
        } else {
          if (n & 0x000000000C000000)
            return (n & 0x0000000008000000) ? 28 : 27;
          else
            return (n & 0x0000000002000000) ? 26 : 25;
        }
      } else {
        if (n & 0x0000000000F00000) {
          if (n & 0x0000000000C00000)
            return (n & 0x0000000000800000) ? 24 : 23;
          else
            return (n & 0x0000000000200000) ? 22 : 21;
        } else {
          if (n & 0x00000000000C0000)
            return (n & 0x0000000000080000) ? 20 : 19;
          else
            return (n & 0x0000000000020000) ? 18 : 17;
        }
      }
    } else {
      if (n & 0x000000000000FF00) {
        if (n & 0x000000000000F000) {
          if (n & 0x000000000000C000)
            return (n & 0x0000000000008000) ? 16 : 15;
          else
            return (n & 0x0000000000002000) ? 14 : 13;
        } else {
          if (n & 0x0000000000000C00)
            return (n & 0x0000000000000800) ? 12 : 11;
          else
            return (n & 0x0000000000000200) ? 10 : 9;
        }
      } else {
        if (n & 0x00000000000000F0) {
          if (n & 0x00000000000000C0)
            return (n & 0x0000000000000080) ? 8 : 7;
          else
            return (n & 0x0000000000000020) ? 6 : 5;
        } else {
          if (n & 0x000000000000000C)
            return (n & 0x0000000000000008) ? 4 : 3;
          else
            return (n & 0x0000000000000002) ? 2 : (n ? 1 : 0);
        }
      }
    }
  }
}

int highest_bit_bs(unsigned long long n)
{
  const unsigned long long mask[] = {
    0x000000007FFFFFFF,
    0x000000000000FFFF,
    0x00000000000000FF,
    0x000000000000000F,
    0x0000000000000003,
    0x0000000000000001
  };
  int hi = 64;
  int lo = 0;
  int i = 0;

  if (n == 0)
    return 0;

  for (i = 0; i < sizeof mask / sizeof mask[0]; i++) {
    int mi = lo + (hi - lo) / 2;

    if ((n >> mi) != 0)
      lo = mi;
    else if ((n & (mask[i] << lo)) != 0)
      hi = mi;
  }

  return lo + 1;
}

int highest_bit_shift(unsigned long long n)
{
  int i = 0;
  for (; n; n >>= 1, i++)
    ; /* empty */
  return i;
}

static int count_ones(unsigned long long d)
{
  d = ((d & 0xAAAAAAAAAAAAAAAA) >>  1) + (d & 0x5555555555555555);
  d = ((d & 0xCCCCCCCCCCCCCCCC) >>  2) + (d & 0x3333333333333333);
  d = ((d & 0xF0F0F0F0F0F0F0F0) >>  4) + (d & 0x0F0F0F0F0F0F0F0F);
  d = ((d & 0xFF00FF00FF00FF00) >>  8) + (d & 0x00FF00FF00FF00FF);
  d = ((d & 0xFFFF0000FFFF0000) >> 16) + (d & 0x0000FFFF0000FFFF);
  d = ((d & 0xFFFFFFFF00000000) >> 32) + (d & 0x00000000FFFFFFFF);
  return d;
}

int highest_bit_parallel(unsigned long long n)
{
  n |= n >> 1;
  n |= n >> 2;
  n |= n >> 4;
  n |= n >> 8;
  n |= n >> 16;
  n |= n >> 32;
  return count_ones(n);
}

int highest_bit_so(unsigned long long x)
{
  static const unsigned long long t[6] = {
    0xFFFFFFFF00000000ull,
    0x00000000FFFF0000ull,
    0x000000000000FF00ull,
    0x00000000000000F0ull,
    0x000000000000000Cull,
    0x0000000000000002ull
  };

  int y = (((x & (x - 1)) == 0) ? 0 : 1);
  int j = 32;
  int i;

  for (i = 0; i < 6; i++) {
    int k = (((x & t[i]) == 0) ? 0 : j);
    y += k;
    x >>= k;
    j >>= 1;
  }

  return y;
}

int highest_bit_so2(unsigned long long value)
{
  int pos = 0;
  if (value & (value - 1ULL))
  {
    pos = 1;
  }
  if (value & 0xFFFFFFFF00000000ULL)
  {
    pos += 32;
    value = value >> 32;
  }
  if (value & 0x00000000FFFF0000ULL)
  {
    pos += 16;
    value = value >> 16;
  }
  if (value & 0x000000000000FF00ULL)
  {
    pos += 8;
    value = value >> 8;
  }
  if (value & 0x00000000000000F0ULL)
  {
    pos += 4;
    value = value >> 4;
  }
  if (value & 0x000000000000000CULL)
  {
    pos += 2;
    value = value >> 2;
  }
  if (value & 0x0000000000000002ULL)
  {
    pos += 1;
    value = value >> 1;
  }
  return pos;
}

This is highest-bit.h:

int highest_bit_unrolled(unsigned long long n);
int highest_bit_bs(unsigned long long n);
int highest_bit_shift(unsigned long long n);
int highest_bit_parallel(unsigned long long n);
int highest_bit_so(unsigned long long n);
int highest_bit_so2(unsigned long long n);

And the main program (sorry about all the copy and paste):

#include <stdlib.h>
#include <stdio.h>
#include <time.h>
#include "highest-bit.h"

double timedelta(clock_t start, clock_t end)
{
  return (end - start)*1.0/CLOCKS_PER_SEC;
}

int main(int argc, char **argv)
{
  int i;
  volatile unsigned long long v;
  clock_t start, end;

  start = clock();

  for (i = 0; i < 10000000; i++) {
    for (v = 0x8000000000000000; v; v >>= 1)
      highest_bit_unrolled(v);
  }

  end = clock();

  printf("highest_bit_unrolled = %6.3fs\n", timedelta(start, end));

  start = clock();

  for (i = 0; i < 10000000; i++) {
    for (v = 0x8000000000000000; v; v >>= 1)
      highest_bit_parallel(v);
  }

  end = clock();

  printf("highest_bit_parallel = %6.3fs\n", timedelta(start, end));

  start = clock();

  for (i = 0; i < 10000000; i++) {
    for (v = 0x8000000000000000; v; v >>= 1)
      highest_bit_bs(v);
  }

  end = clock();

  printf("highest_bit_bs = %6.3fs\n", timedelta(start, end));

  start = clock();

  for (i = 0; i < 10000000; i++) {
    for (v = 0x8000000000000000; v; v >>= 1)
      highest_bit_shift(v);
  }

  end = clock();

  printf("highest_bit_shift = %6.3fs\n", timedelta(start, end));

  start = clock();

  for (i = 0; i < 10000000; i++) {
    for (v = 0x8000000000000000; v; v >>= 1)
      highest_bit_so(v);
  }

  end = clock();

  printf("highest_bit_so = %6.3fs\n", timedelta(start, end));

  start = clock();

  for (i = 0; i < 10000000; i++) {
    for (v = 0x8000000000000000; v; v >>= 1)
      highest_bit_so2(v);
  }

  end = clock();

  printf("highest_bit_so2 = %6.3fs\n", timedelta(start, end));

  return 0;
}

I've tried this various Intel x86 boxes, old and new.

The highest_bit_unrolled (unrolled binary search) is consistently significantly faster than highest_bit_parallel (branch-free bit ops). This is faster than than highest_bit_bs (binary search loop), and that in turn is faster than highest_bit_shift (naive shift and count loop).

highest_bit_unrolled is also faster than the one in the accepted SO answer (highest_bit_so) and also one given in another answer (highest_bit_so2).

The benchmark cycles through a one-bit mask that covers successive bits. This is to try to defeat branch prediction in the unrolled binary search, which is realistic: in a real world program, the input cases are unlikely to exhibit locality of bit position.

Here are results on an old Intel(R) Core(TM)2 Duo CPU E4500 @ 2.20GHz:

$ ./highest-bit
highest_bit_unrolled =  6.090s
highest_bit_parallel =  9.260s
highest_bit_bs = 19.910s
highest_bit_shift = 21.130s
highest_bit_so =  8.230s
highest_bit_so2 =  6.960s

On a newer model Intel(R) Core(TM) i7-6700K CPU @ 4.00GHz:

highest_bit_unrolled =  1.555s
highest_bit_parallel =  3.420s
highest_bit_bs =  6.486s
highest_bit_shift =  9.505s
highest_bit_so =  4.127s
highest_bit_so2 =  1.645s

On the newer hardware, highest_bit_so2 comes closer to highest_bit_unrolled on the newer hardware. The order is not quite the same; now highest_bit_so really falls behind, and is slower than highest_bit_parallel.

The fastest one, highest_bit_unrolled contains the most code with the most branches. Every single return value reached by a different set of conditions with its own dedicated piece of code.

The intuition of "avoid all branches" (due to worries about branch mis-predictions) is not always right. Modern (and even not-so-modern any more) processors contain considerable cunning in order not to be hampered by branching.


P.S. the highest_bit_unrolled was introduced in the TXR language in December 2011 (with mistakes, since debugged).

Recently, I started wondering whether some nicer, more compact code without branches might not be faster.

I'm somewhat surprised by the results.

Arguably, the code should really be #ifdef-ing for GNU C and using some compiler primitives, but as far as portability goes, that version stays.

A sort of binary search for the highest order bit.

A good reference on bitwise operation hacks posted on stanford.edu: Bit Twiddling Hacks. There are a lot of examples in that reference for many things, including several algorithms for computing log2 values. This includes similar - and maybe better - approaches to this one.

I had the idea that a sort of binary search for the most significant bit would be faster than the max 60+ loops of the simple solution. The algorithm shifts n by half of it's last shift value each iteration to zero in on the msb.

For 64-bit values, it requires 6 iterations to determine the position of the most significant bit - which is the floor log2 of n.

The performance of this approach, as I measured it, wasn't much better than the simple loop in the 2nd example below.

// floor log2 in 6 max iterations for unsigned 64-bit values.
//
int floor_log2(unsigned long long n)
{
    int shval = 32;
    int msb   =  0;
    while (shval) {
        if (n >> shval) {
            msb += shval;
            n  >>= shval;
        }
        shval >>= 1;
    }
    return msb;
}

int ceil_log2(unsigned long long n) 
{
    return floor_log2(2 * n - 1);
}

Simple algorithm:

int floor_log2_simple(unsigned long long n)
{
    int c =  0;
    while (n) {
        n >>= 1;
        c++;
    }
    return c - 1;
}

These same algorithms implemented in Rust performed nearly identically. The top approach performed a tiny amount better - but by only by a couple nano seconds on average. No amount of tweaking the algorithm gave any better performance.

meh...

Related