Fast computing of log2 for 64-bit integers
C64 BitBit Manipulation32bit 64bitLookupC Problem Overview
A great programming resource, Bit Twiddling Hacks, proposes (here) the following method to compute log2 of a 32-bit integer:
#define LT(n) n, n, n, n, n, n, n, n, n, n, n, n, n, n, n, n
static const char LogTable256[256] =
{
-1, 0, 1, 1, 2, 2, 2, 2, 3, 3, 3, 3, 3, 3, 3, 3,
LT(4), LT(5), LT(5), LT(6), LT(6), LT(6), LT(6),
LT(7), LT(7), LT(7), LT(7), LT(7), LT(7), LT(7), LT(7)
};
unsigned int v; // 32-bit word to find the log of
unsigned r; // r will be lg(v)
register unsigned int t, tt; // temporaries
if (tt = v >> 16)
{
r = (t = tt >> 8) ? 24 + LogTable256[t] : 16 + LogTable256[tt];
}
else
{
r = (t = v >> 8) ? 8 + LogTable256[t] : LogTable256[v];
}
and mentions that
> The lookup table method takes only about 7 operations to find the log > of a 32-bit value. If extended for 64-bit quantities, it would take > roughly 9 operations.
but, alas, doesn't give any additional info about which way one should actually go to extend the algorithm to 64-bit integers.
Any hints about how a 64-bit algorithm of this kind would look like?
C Solutions
Solution 1 - C
Intrinsic functions are really fast, but still are insufficient for a truly cross-platform, compiler-independent implementation of log2. So in case anyone is interested, here is the fastest, branch-free, CPU-abstract DeBruijn-like algorithm I've come to while researching the topic on my own.
const int tab64[64] = {
63, 0, 58, 1, 59, 47, 53, 2,
60, 39, 48, 27, 54, 33, 42, 3,
61, 51, 37, 40, 49, 18, 28, 20,
55, 30, 34, 11, 43, 14, 22, 4,
62, 57, 46, 52, 38, 26, 32, 41,
50, 36, 17, 19, 29, 10, 13, 21,
56, 45, 25, 31, 35, 16, 9, 12,
44, 24, 15, 8, 23, 7, 6, 5};
int log2_64 (uint64_t value)
{
value |= value >> 1;
value |= value >> 2;
value |= value >> 4;
value |= value >> 8;
value |= value >> 16;
value |= value >> 32;
return tab64[((uint64_t)((value - (value >> 1))*0x07EDD5E59A4E28C2)) >> 58];
}
The part of rounding down to the next lower power of 2 was taken from Power-of-2 Boundaries and the part of getting the number of trailing zeros was taken from BitScan (the (bb & -bb) code there is to single out the rightmost bit that is set to 1, which is not needed after we've rounded the value down to the next power of 2).
And the 32-bit implementation, by the way, is
const int tab32[32] = {
0, 9, 1, 10, 13, 21, 2, 29,
11, 14, 16, 18, 22, 25, 3, 30,
8, 12, 20, 28, 15, 17, 24, 7,
19, 27, 23, 6, 26, 5, 4, 31};
int log2_32 (uint32_t value)
{
value |= value >> 1;
value |= value >> 2;
value |= value >> 4;
value |= value >> 8;
value |= value >> 16;
return tab32[(uint32_t)(value*0x07C4ACDD) >> 27];
}
As with any other computational method, log2 requires the input value to be greater than zero.
Solution 2 - C
If you are using GCC, a lookup table is unnecessary in this case.
GCC provides a builtin function to determine the amount of leading zeros:
> Built-in Function: int __builtin_clz (unsigned int x)
> Returns the number of leading 0-bits in x, starting at the most significant bit position. If x is 0, the result is undefined.
So you can define:
#define LOG2(X) ((unsigned) (8*sizeof (unsigned long long) - __builtin_clzll((X)) - 1))
and it will work for any unsigned long long int. The result is rounded down.
For x86 and AMD64 GCC will compile it to a bsr instruction, so the solution is very fast (much faster than lookup tables).
#include <stdio.h>
#define LOG2(X) ((unsigned) (8*sizeof (unsigned long long) - __builtin_clzll((X)) - 1))
int main(void) {
unsigned long long input;
while (scanf("%llu", &input) == 1) {
printf("log(%llu) = %u\n", input, LOG2(input));
}
return 0;
}
Compiled output: https://godbolt.org/z/16GnjszMs
Solution 3 - C
I was trying to convert Find the log base 2 of an N-bit integer in O(lg(N)) operations with multiply and lookup to 64-bit by brute forcing the magic number. Needless to say it was taking a while.
I then found Desmond's answer and decided to try his magic number as a start point. Since I have a 6 core processor I ran it in parallel starting at 0x07EDD5E59A4E28C2 / 6 multiples. I was surprised it found something immediately. Turns out 0x07EDD5E59A4E28C2 / 2 worked.
So here is the code for 0x07EDD5E59A4E28C2 which saves you a shift and subtract:
int LogBase2(uint64_t n)
{
static const int table[64] = {
0, 58, 1, 59, 47, 53, 2, 60, 39, 48, 27, 54, 33, 42, 3, 61,
51, 37, 40, 49, 18, 28, 20, 55, 30, 34, 11, 43, 14, 22, 4, 62,
57, 46, 52, 38, 26, 32, 41, 50, 36, 17, 19, 29, 10, 13, 21, 56,
45, 25, 31, 35, 16, 9, 12, 44, 24, 15, 8, 23, 7, 6, 5, 63 };
n |= n >> 1;
n |= n >> 2;
n |= n >> 4;
n |= n >> 8;
n |= n >> 16;
n |= n >> 32;
return table[(n * 0x03f6eaf2cd271461) >> 58];
}
Solution 4 - C
Base-2 Integer Logarithm
Here is what I do for 64-bit unsigned integers. This calculates the floor of the base-2 logarithm, which is equivalent to the index of the most significant bit. This method is smokingly fast for large numbers because it uses an unrolled loop that executes always in log₂64 = 6 steps.
Essentially, what it does is subtracts away progressively smaller squares in the sequence { 0 ≤ k ≤ 5: 2^(2^k) } = { 2³², 2¹⁶, 2⁸, 2⁴, 2², 2¹ } = { 4294967296, 65536, 256, 16, 4, 2, 1 } and sums the exponents k of the subtracted values.
int uint64_log2(uint64_t n)
{
#define S(k) if (n >= (UINT64_C(1) << k)) { i += k; n >>= k; }
int i = -(n == 0); S(32); S(16); S(8); S(4); S(2); S(1); return i;
#undef S
}
Note that this returns –1 if given the invalid input of 0 (which is what the initial -(n == 0) is checking for). If you never expect to invoke it with n == 0, you could substitute int i = 0; for the initializer and add assert(n != 0); at entry to the function.
Base-10 Integer Logarithm
Base-10 integer logarithms can be calculated using similarly — with the largest square to test being 10¹⁶ because log₁₀2⁶⁴ ≅ 19.2659... (Note: This is not the fastest way to accomplish a base-10 integer logarithm, because it uses integer division, which is inherently slow. A faster implementation would be to use an accumulator with values that grow exponentially, and compare against the accumulator, in effect doing a sort of binary search.)
int uint64_log10(uint64_t n)
{
#define S(k, m) if (n >= UINT64_C(m)) { i += k; n /= UINT64_C(m); }
int i = -(n == 0);
S(16,10000000000000000); S(8,100000000); S(4,10000); S(2,100); S(1,10);
return i;
#undef S
}
Solution 5 - C
Here's a pretty compact and fast extension, using no additional temporaries:
r = 0;
/* If its wider than 32 bits, then we already know that log >= 32.
So store it in R. */
if (v >> 32)
{
r = 32;
v >>= 32;
}
/* Now do the exact same thing as the 32 bit algorithm,
except we ADD to R this time. */
if (tt = v >> 16)
{
r += (t = tt >> 8) ? 24 + LogTable256[t] : 16 + LogTable256[tt];
}
else
{
r += (t = v >> 8) ? 8 + LogTable256[t] : LogTable256[v];
}
Here is one built with a chain of ifs, again using no additional temporaries. Might not be the fastest though.
if (tt = v >> 48)
{
r = (t = tt >> 8) ? 56 + LogTable256[t] : 48 + LogTable256[tt];
}
else if (tt = v >> 32)
{
r = (t = tt >> 8) ? 40 + LogTable256[t] : 32 + LogTable256[tt];
}
else if (tt = v >> 16)
{
r = (t = tt >> 8) ? 24 + LogTable256[t] : 16 + LogTable256[tt];
}
else
{
r = (t = v >> 8) ? 8 + LogTable256[t] : LogTable256[v];
}
Solution 6 - C
If you were looking for the c++ answer and you got here and since it boils down to counting zeroes then you got the std::countl_zero which according to godbolt.org calls bsr.
std::countl_zero is available from C++20, you may need to add -std=gnu++2a to your compiler command line
Solution 7 - C
The algorithm basically finds out which byte contains the most significant 1 bit, and then looks up that byte in the lookup to find the log of the byte, then adds it to the position of the byte.
Here is a somewhat simplified version of the 32-bit algorithm:
if (tt = v >> 16)
{
if (t = tt >> 8)
{
r = 24 + LogTable256[t];
}
else
{
r = 16 + LogTable256[tt];
}
}
else
{
if (t = v >> 8)
{
r = 8 + LogTable256[t];
}
else
{
r = LogTable256[v];
}
}
This is the equivalent 64-bit algorithm:
if (ttt = v >> 32)
{
if (tt = ttt >> 16)
{
if (t = tt >> 8)
{
r = 56 + LogTable256[t];
}
else
{
r = 48 + LogTable256[tt];
}
}
else
{
if (t = ttt >> 8)
{
r = 40 + LogTable256[t];
}
else
{
r = 32 + LogTable256[ttt];
}
}
}
else
{
if (tt = v >> 16)
{
if (t = tt >> 8)
{
r = 24 + LogTable256[t];
}
else
{
r = 16 + LogTable256[tt];
}
}
else
{
if (t = v >> 8)
{
r = 8 + LogTable256[t];
}
else
{
r = LogTable256[v];
}
}
}
I came up with an algorithm for any size types that I think is nicer than the original.
unsigned int v = 42;
unsigned int r = 0;
unsigned int b;
for (b = sizeof(v) << 2; b; b = b >> 1)
{
if (v >> b)
{
v = v >> b;
r += b;
}
}
Note: b = sizeof(v) << 2 sets b to half the number of bits in v. I used shifting instead of multiplication here (just because I felt like it).
You could add a lookup table to that algorithm to speed it up possibly, but it's more a proof-of-concept.
Solution 8 - C
Take this:
typedef unsigned int uint;
typedef uint64_t ulong;
uint as_uint(const float x) {
return *(uint*)&x;
}
ulong as_ulong(const double x) {
return *(ulong*)&x;
}
uint log2_fast(const uint x) {
return (uint)((as_uint((float)x)>>23)-127);
}
uint log2_fast(const ulong x) {
return (uint)((as_ulong((double)x)>>52)-1023);
}
How it works:
The input integer x is casted to float and then reinterpret as bits. The IEEE float format stores the exponent in bits 30-23 as an integer with bias 127, so by shifting it 23 bits to the right and subtracting the bias, we get log2(x).
For a 64-bit integer input, x is casted to double, for which the exponent is in bits 62-52 (shift 52 bits to the right) and the exponent bias is 1023.
Solution 9 - C
Here is a slightly modified one from SPWorley on 3/22/2009 (see post for details)
double ff=(double)(v|1);
return ((*(1+(uint32_t *)&ff))>>52)-1023; // assumes x86 endianness
Side-note: some users pointed out that when compiled in some situations this may result in an incorrect answer.