67

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?

7
  • 4
    @dbaupp I've got bags of ifs of all possible kinds, sorts, and taste, which ones would do best? Jul 7, 2012 at 15:41
  • 10
    That's just an academical question, right? Otherwise just use _BitScanReverse64 (msvc) or __builtin_clzll (gcc)
    – harold
    Jul 7, 2012 at 15:42
  • Ones like the ones you already have. (Using the most naive extension, it'll look something like if (tt = v >> 48) { ... } else if (tt = v >> 32) { ... } ..., although this will perform marginally worse than the proper binary search Kendall correctly suggests.)
    – huon
    Jul 7, 2012 at 15:43
  • 3
    @harold: No, this is not an academical question at all. Even if someone decides to use compiler-specific intrisincs, they will go into compiler-specific #if branches. This, of course, does not eliminate the need for a "default" branch implemented more or less universally. Jul 7, 2012 at 16:39
  • 1
    @AndreyT it would be interesting if people started doing that. Code might actually become real-life-portable, instead of Ivory-Tower-portable (where a sensible implementation of int can not be relied upon, but gcc-specific language extenstion can)
    – harold
    Jul 7, 2012 at 16:59

9 Answers 9

95

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.

10
  • 2
    @ArjunShankar You're welcome. However, I still think there is a way of shaving off those two ops in the table lookup line, namely subtraction and right shift, by means of generating another perfect hashing function. Don't know if I'll have enough time for this pursuing of zen as long as my main compilers are MSVC and GCC ;) Jul 9, 2012 at 16:29
  • 1
    @ArjunShankar And, for the clarity of where the operations, table entries, and constants come from, I've updated the answer to cite the sources. Jul 9, 2012 at 16:50
  • 7
    To help with the trade off between portability and speed, on an Intel(R) Xeon(R) CPU X5650 @ 2.67GHz the lookup table is on average about 4 times slower than the intrinsic (about 13 cycles vs 4 cycles)
    – Come Raczy
    Aug 15, 2013 at 13:23
  • 1
    @DesmondHume — Since, as you say, the input value must be greater than zero, maybe you want to add a line at the beginning of the function that says assert(value > 0);. This not only helps catch any errors but also serves to help document the entry conditions. Apr 21, 2016 at 14:39
  • 6
    @DesmondHume — Another observation: Instead of having the tab64[] table contain int values, why not have it contain 8-bit values and specify the table as const int8_t tab64[64]... This way, the table requires only 64 bytes (instead of 256). Additionally, if you include __attribute__((aligned(64))), then it is also guaranteed to fit within a single cache line on a modern processor. As currently written (with 256 bytes and no alignment), the table could be using up to 5 cache lines. Apr 21, 2016 at 14:44
66

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).

Working example:

#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

8
  • 2
    How about also handling "If x is 0, the result is undefined." in your example? :) Jul 7, 2012 at 16:41
  • 2
    @ArjunShankar actually I thought about it, but couldn't think of an appropriate integer for that case. ;) I leave it to the interested reader to add an if-else case to the macro. (Also only the result will be undefined [most likely 0], but there won't be a crash if a 0 was supplied.)
    – Kijewski
    Jul 7, 2012 at 16:46
  • @kay Didn't know about bsr instruction. Wanted it to be more CPU-independent tho. Thanks. Jul 7, 2012 at 17:27
  • 15
    @DesmondHume __builtin_clz is not processor specific. GCC will find a set of instructions for that will perform well for any supported architecture.
    – Kijewski
    Jul 7, 2012 at 17:36
  • 1
    If you compile for Haswell or later, GCC and Clang select LZCNT instead of BSR. Clang then even generates just lzcnt rax, rdi; xor eax, 63 whereas GCC or ICC select 1 or 2 extra instructions, for this integer log2 function. Sep 24, 2020 at 10:49
22

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];
}
7
  • Are you 100% sure this is correct in all cases? I haven't really wrapped my head around how this works internally yet, so I don't know if there is a way to prove correctness...
    – Markus A.
    Oct 3, 2014 at 17:54
  • 3
    It's correct for all inputs except 0. It returns 0 for 0 which may be valid for what you're using it for. The lines with the shifts round n up to 1 less than the next power of 2. It basically sets all bits after the leading 1 bit to 1. This reduces all possible inputs to 64 possible values: 0x0, 0x1, 0x3, 0x7, 0xf, 0x1f, 0x3f, etc. Multiplying those 64 values with the number 0x03f6eaf2cd271461 gives you another 64 unique values in the top 6 bits. The shift by 58 just positions those 6 bits for use as an index into table.
    – Avernar
    Oct 7, 2014 at 1:27
  • That makes perfect sense. Thank you. Somehow I had a brain-block and I was reading the code as n = n | (n >> 1) | (n >> 2) and so forth and was trying to figure out how many cases I would have to check to verify correctness... D'Uh... Always good to know why something works! :) +1
    – Markus A.
    Oct 7, 2014 at 3:24
  • @Avernar — As I was just pointing out to Desmond Hume in his nearly identical solution, if you make your table contain elements of type int8_t instead of int, and align it to a 64-byte boundary, then the table will be only 64 bytes instead of 256 bytes and will fit within a single cache line on modern processors. Apr 21, 2016 at 14:49
  • 3
    I guess this works because your "magic number" m is a DeBruin sequence B(2, 6) that starts with six 0s followed by six 1s. With such a sequence, multiplying with ns of the form 00011111 instead of 00010000 still works. Now the operation is m * (2^(n+1) - 1), or m * 2^(n+1) - m, instead of m * 2^n. Since the sequence starts with six 0s, multiplying by 2^(n+1) still results in a perfect hash. Because of the following six 1s, the subtraction always propagates a carry bit into the upper six bits of the result.
    – nwellnhof
    May 30, 2017 at 0:29
11

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
}
2
  • What confused me about this answer is, you are saying you subtract squares, which wouldn't make sense. What this does is it divides by squares. For anyone who is also confused by this: Say your number is n = 2^49 => S(32) divides by 2^32, now i = 32 and n = 2^17 => S(16) divides by 2^16, now i = 48 and n = 2^1 => S(8), S(4) and S(2) do nothing. => S(1) divides by 2^1, now i = 49 and n = 2^0 => return i Jul 7, 2022 at 9:42
  • Subtracting is actually correct (except it's actually doing addition up to a ceiling rather than actual subtraction); there is no dividing by squares. The only "dividing" that happens is that n is divided by a power of 2 at each step. Where the "squares" come in is that, yes, it's dividing by powers of 2, but each successive power of 2 that it divides by is the square of the previous power of 2, e.g., 2^1, 2^2, 2^4, 2^8, 2^16, 2^32, etc. Hence, in effect, it's subtracting (but actually adding) squares that happen to also be powers of two. Jul 7, 2022 at 17:27
5

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];
    }
3
  • 2
    (unsigned long long or uint64_t)
    – huon
    Jul 7, 2012 at 15:58
  • 1
    @dbaupp - I was just too lazy to write a main and include stdint.h. Thanks for the nudge. I tried it in the meanwhile, and it's working fine. Jul 7, 2012 at 16:02
  • @ArjunShankar The first algorithm is hands down awesome. Thank you. Jul 7, 2012 at 16:53
5

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

1
  • 3
    This same C++ version also adds the function std::bit_width(x) which results in 1 + ⎣log2(x)⎦ if x > 0 and 0 if x == 0.
    – Lukas
    Nov 29, 2020 at 21:19
3

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.

4
  • 1
    Just personally, I think the more compact ternary version is "simpler": takes up less space. :)
    – huon
    Jul 7, 2012 at 15:45
  • @dbaupp: Depends on your point of view. I expanded the ternary so it was easier to see what was going on. Jul 7, 2012 at 15:48
  • @KendallFrey Thank you, but would the fourth table look-up, if count from the beginning of the 64-bit algorithm, be able to overrun the boundaries of the table? Jul 7, 2012 at 15:53
  • @DesmondHume: Yes, I believe so. Copy-n-paste error here. Fixed. Jul 7, 2012 at 15:56
1

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.

2
  • The 2 typedef are not needed and typedef unsigned int uint; is incorrect when unsigned is 16-bit. Consider using uint64_t, uint32_t - more clear anyways. Jan 2, 2021 at 22:48
  • 1
    The biggest problem with this is that you can only calculate the log of integer numbers up to 55 bits (in case of ulong) and 25 bits (in case of int). Feb 18 at 3:26
0

Here is a modified one from SPWorley on 3/22/2009

supports >55 bits and 0 returns 0. To have it return a negitive number instead of 0 remove the "| 1".

int Log2(uint64_t  v) // assumes x86 endianness
{
    int extra = 0;
    if (v > 1023)
    {
        v >>= 10;
        extra = 10;
    }

    double ff = (double)(v | 1);
    uint32_t result = ((*(1 + (uint32_t*)&ff)) >> 52) - 1023;  
    
    return result + extra;
}

int Log2(uint32_t  v) // assumes x86 endianness
{
    double ff = (double)(v | 1);
    return ((*(1 + (uint32_t*)&ff)) >> 52) - 1023;  
}

A testing function and the output...

int main()
{
    std::cout << ((uint64_t)0) << ": " << Log2((uint64_t)0) << std::endl;

    for (size_t i = 1; i < 64; i++)
        for (size_t j = 2; j-- > 0;)
            std::cout << ((uint64_t)1 << i) - j << ": " << Log2(((uint64_t)1 << i) - j) << std::endl;
}

/**** Output  ****/
0: 0
1: 0
2: 1
3: 1
4: 2
7: 2
8: 3
15: 3
16: 4
31: 4
32: 5
63: 5
64: 6
127: 6
128: 7
255: 7
256: 8
511: 8
512: 9
1023: 9
1024: 10
2047: 10
2048: 11
4095: 11
4096: 12
8191: 12
8192: 13
16383: 13
16384: 14
32767: 14
32768: 15
65535: 15
65536: 16
131071: 16
131072: 17
262143: 17
262144: 18
524287: 18
524288: 19
1048575: 19
1048576: 20
2097151: 20
2097152: 21
4194303: 21
4194304: 22
8388607: 22
8388608: 23
16777215: 23
16777216: 24
33554431: 24
33554432: 25
67108863: 25
67108864: 26
134217727: 26
134217728: 27
268435455: 27
268435456: 28
536870911: 28
536870912: 29
1073741823: 29
1073741824: 30
2147483647: 30
2147483648: 31
4294967295: 31
4294967296: 32
8589934591: 32
8589934592: 33
17179869183: 33
17179869184: 34
34359738367: 34
34359738368: 35
68719476735: 35
68719476736: 36
137438953471: 36
137438953472: 37
274877906943: 37
274877906944: 38
549755813887: 38
549755813888: 39
1099511627775: 39
1099511627776: 40
2199023255551: 40
2199023255552: 41
4398046511103: 41
4398046511104: 42
8796093022207: 42
8796093022208: 43
17592186044415: 43
17592186044416: 44
35184372088831: 44
35184372088832: 45
70368744177663: 45
70368744177664: 46
140737488355327: 46
140737488355328: 47
281474976710655: 47
281474976710656: 48
562949953421311: 48
562949953421312: 49
1125899906842623: 49
1125899906842624: 50
2251799813685247: 50
2251799813685248: 51
4503599627370495: 51
4503599627370496: 52
9007199254740991: 52
9007199254740992: 53
18014398509481983: 53
18014398509481984: 54
36028797018963967: 54
36028797018963968: 55
72057594037927935: 55
72057594037927936: 56
144115188075855871: 56
144115188075855872: 57
288230376151711743: 57
288230376151711744: 58
576460752303423487: 58
576460752303423488: 59
1152921504606846975: 59
1152921504606846976: 60
2305843009213693951: 60
2305843009213693952: 61
4611686018427387903: 61
4611686018427387904: 62
9223372036854775807: 62
9223372036854775808: 63
2
  • 2
    It cannot compute log of numbers above 55 bits. Feb 18 at 3:27
  • 1
    @HenriqueBucher - Thanks for the suggestion. I added support for over 55 bits in length by temporarily downshifting. It also used the same 1023 so the compiler just needs to store/load it once. It also handles a zero (by returning a zero). Feb 18 at 7:14

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.