快速计算 64 位整数的 log2答案

【问题标题】：Fast computing of log2 for 64-bit integers快速计算 64 位整数的 log2
【发布时间】：2012-07-07 17:51:17
【问题描述】：

一个很棒的编程资源 Bit Twiddling Hacks 提出 (here) 以下方法来计算 32 位整数的 log2：

#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];
}

并提到

查表法只需要大约7次操作即可找到日志一个 32 位的值。如果扩展到 64 位数量，则需要大约 9 次操作。

但是，遗憾的是，没有提供任何额外信息，说明应该采用哪种方式将算法扩展到 64 位整数。

关于这种 64 位算法的外观有什么提示吗？

【问题讨论】：

@dbaupp 我有一袋 ifs 的所有可能种类、种类和口味，哪一个最好？
这只是一个学术问题，对吧？否则只需使用 _BitScanReverse64 (msvc) 或 __builtin_clzll (gcc)
你已经拥有的那些。（使用最简单的扩展，它看起来像if (tt = v >> 48) { ... } else if (tt = v >> 32) { ... } ...，尽管这会比正确的二分搜索 Kendall 正确建议的效果稍差。）
@harold：不，这根本不是学术问题。即使有人决定使用特定于编译器的 intrisinc，他们也会进入特定于编译器的 #if 分支。当然，这并不能消除或多或少普遍实现的“默认”分支的需要。
@AndreyT 如果人们开始这样做会很有趣。代码实际上可能成为现实生活中可移植的，而不是 Ivory-Tower 可移植（不能依赖 int 的合理实现，但 gcc 特定的语言扩展可以）

标签： c 64-bit bit-manipulation 32bit-64bit lookup

【解决方案1】：

该算法基本上是找出哪个字节包含最高有效1位，然后在查找中查找该字节以找到该字节的日志，然后将其添加到该字节的位置。

这里是 32 位算法的简化版本：

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];
    }
}

这是等效的 64 位算法：

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];
        }
    }
}

我想出了一个算法，适用于我认为比原来更好的任何尺寸类型。

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;
    }
}

注意：b = sizeof(v) << 2 将 b 设置为 v 中位数的一半。我在这里使用移位而不是乘法（只是因为我喜欢它）。

您可以向该算法添加一个查找表以加快其速度，但这更像是一个概念验证。

【讨论】：

就个人而言，我认为更紧凑的三元版本是“更简单”：占用更少的空间。 :)
@dbaupp：取决于你的观点。我扩展了三元组，以便更容易看到发生了什么。
@KendallFrey 谢谢，但是第四次查表，如果从 64 位算法开始算起，会不会超出表的边界？
@DesmondHume：是的，我相信是的。此处复制粘贴错误。固定。

【解决方案2】：

这是一个非常紧凑的和快速扩展，不使用额外的临时：

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];
  }

这里是用ifs 链构建的，同样没有使用额外的临时对象。不过可能不是最快的。

  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];
    }

【讨论】：

（unsigned long long 或 uint64_t）
@dbaupp - 我懒得写main 并包含stdint.h。谢谢你的推动。同时我也试过了，效果很好。
@ArjunShankar 第一个算法非常棒。谢谢。

【解决方案3】：

如果您使用的是 GCC，则在这种情况下不需要查找表。

GCC 提供了一个内置函数来确定前导零的数量：

内置函数：int __builtin_clz (unsigned int x)
返回 x 中前导 0 位的数量，从最高有效位位置开始。如果 x 为 0，则结果未定义。

所以你可以定义：

#define LOG2(X) ((unsigned) (8*sizeof (unsigned long long) - __builtin_clzll((X)) - 1))

它适用于任何 unsigned long long int。结果四舍五入。

对于 x86 和 AMD64，GCC 会将其编译为 bsr 指令，因此解决方案非常快（比查找表快得多）。

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;
}

编译输出：https://godbolt.org/z/16GnjszMs

【讨论】：

如何处理“如果 x 为 0，则结果未定义”。在你的例子中？ :)
@ArjunShankar 实际上我想过，但想不出适合这种情况的整数。 ;) 我留给感兴趣的读者在宏中添加 if-else 案例。（也只有结果是未定义的 [最可能是 0]，但如果提供了 0，则不会发生崩溃。）
@kay 不知道bsr 指令。希望它更加独立于 CPU。谢谢。
@DesmondHume __builtin_clz 不是特定于处理器的。 GCC 会找到一组指令，这些指令适用于任何受支持的架构。
如果为 Haswell 或更高版本编译，GCC 和 Clang 选择 LZCNT 而不是 BSR。 Clang 然后甚至 generates 只是 lzcnt rax, rdi; xor eax, 63 而 GCC 或 ICC 为这个整数 log2 函数选择 1 或 2 条额外指令。

【解决方案4】：

内部函数非常快，但对于真正跨平台、独立于编译器的 log2 实现来说仍然不够。因此，如果有人感兴趣，这里是我在自己研究该主题时遇到的最快、无分支、类似 CPU 抽象的 DeBruijn 算法。

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];
}

向下舍入到 2 的下一个较低幂的部分取自 Power-of-2 Boundaries，获取尾随零数量的部分取自 BitScan（(bb & -bb) 代码用于挑出最右边的设置为 1 的位，在我们将该值向下舍入到 2 的下一个幂之后就不需要了。

顺便说一下，32 位实现是

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];
}

与任何其他计算方法一样，log2 要求输入值大于零。

【讨论】：

@ArjunShankar 不客气。但是，我仍然认为有一种方法可以通过生成另一个完美的散列函数来减少查表行中的这两个操作，即减法和右移。只要我的主要编译器是MSVC和GCC，不知道我是否有足够的时间来追求禅；）
@ArjunShankar 而且，为了清楚操作、表条目和常量的来源，我更新了答案以引用来源。
为了帮助在便携性和速度之间取得平衡，在 Intel(R) Xeon(R) CPU X5650 @ 2.67GHz 上，查找表平均约为 4 倍比内在慢（大约 13 个周期对 4 个周期）
@DesmondHume — 因为，正如你所说，输入值必须大于零，也许你想在函数的开头添加一行 assert(value > 0);。这不仅有助于发现任何错误，还有助于记录进入条件。
@DesmondHume — 另一个观察结果：与其让tab64[] 表包含int 值，不如让它包含8 位值并将表指定为const int8_t tab64[64]... 这种方式，该表只需要 64 个字节（而不是 256 个）。此外，如果您包含__attribute__((aligned(64)))，那么它也可以保证适合现代处理器上的单个高速缓存行。正如当前所写的那样（256 字节并且没有对齐），该表最多可以使用 5 个缓存行。

【解决方案5】：

我试图通过强制使用幻数将Find the log base 2 of an N-bit integer in O(lg(N)) operations with multiply and lookup 转换为 64 位。不用说这需要一段时间。

然后我找到了戴斯蒙德的答案，并决定尝试他的幻数作为起点。由于我有一个 6 核处理器，我从 0x07EDD5E59A4E28C2 / 6 倍数开始并行运行它。我很惊讶它立即发现了一些东西。结果 0x07EDD5E59A4E28C2 / 2 有效。

所以这里是 0x07EDD5E59A4E28C2 的代码，它为您节省了移位和减法：

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];
}

【讨论】：

您是否 100% 确定这在所有情况下都是正确的？我还没有真正了解它在内部是如何工作的，所以我不知道是否有办法证明正确性......
除 0 以外的所有输入都是正确的。它返回 0 表示 0，这可能对您使用它的目的有效。带有移位的行将 n 向上舍入到小于 2 的下一个幂的 1。它基本上将前 1 位之后的所有位设置为 1。这将所有可能的输入减少到 64 个可能的值：0x0、0x1、0x3、0x7、0xf , 0x1f, 0x3f 等。将这 64 个值与数字 0x03f6eaf2cd271461 相乘会在前 6 位中得到另外 64 个唯一值。移位 58 只是将这 6 位定位为用作表中的索引。
这很有意义。谢谢你。不知何故，我有一个脑块，我正在阅读代码为 n = n | (n >> 1) | (n >> 2) 等等，并试图弄清楚我需要检查多少个案例才能验证正确性...... 呃......总是很高兴知道为什么某些东西有效！ :) +1
@Avernar — 正如我刚刚在 Desmond Hume 在他几乎相同的解决方案中指出的那样，如果您让表包含类型为 int8_t 而不是 int 的元素，并将其对齐到 64 -byte 边界，则该表将只有 64 字节而不是 256 字节，并且将适合现代处理器上的单个缓存行。
我想这是可行的，因为您的“幻数”m 是一个DeBruin sequence B(2, 6)，它以六个0s 开头，后跟六个1s。使用这样的序列，与 00011111 形式的 ns 相乘而不是 00010000 仍然有效。现在操作是m * (2^(n+1) - 1)，或m * 2^(n+1) - m，而不是m * 2^n。由于该序列以六个0s 开头，因此乘以2^(n+1) 仍会产生完美的散列。由于后面有六个1s，因此减法总是会将进位位传播到结果的高六位。

【解决方案6】：

Base-2 整数对数

这是我为 64 位无符号整数所做的。这会计算以 2 为底的对数的下限，它相当于最高有效位的索引。对于大量数字，此方法快得冒烟，因为它使用了一个展开的循环，该循环始终以 log₂64 = 6 步执行。

本质上，它所做的是在序列 { 0 ≤ k ≤ 5: 2^(2^k) } = { 2³², 2¹⁶, 2⁸, 2⁴, 2², 2¹ } = { 4294967296, 65536, 256, 16, 4, 2, 1 } 并将减法值的指数 k 相加。

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
}

请注意，如果给定无效输入 0（这是初始 -(n == 0) 正在检查的内容），则返回 –1。如果您从没想过会使用n == 0 调用它，您可以用int i = 0; 替换初始化程序，并在函数入口处添加assert(n != 0);。

以 10 为底的整数对数

可以使用类似的方法计算以 10 为底的整数对数 - 要测试的最大平方为 10¹⁶，因为 log₁₀2⁶⁴ ≅ 19.2659...（注意：这不是完成以 10 为底的整数对数的最快方法，因为它使用整数除法，这本来就很慢。更快的实现是使用具有指数增长的累加器，并与累加器进行比较，实际上是进行一种二进制搜索。）

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
}

【讨论】：

【解决方案7】：

拿着这个：

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);
}

它是如何工作的：输入整数x 被转换为float，然后重新解释为位。 IEEE float 格式将 30-23 位中的指数存储为具有偏差 127 的整数，因此通过将其向右移动 23 位并减去偏差，我们得到 log2(x)。对于 64 位整数输入，x 被强制转换为 double，其指数在 62-52 位（右移 52 位），指数偏差为 1023。

【讨论】：

当unsigned 是 16 位时，不需要 2 个 typedef 并且 typedef unsigned int uint; 不正确。考虑使用uint64_t, uint32_t - 反正更清楚。

【解决方案8】：

这里是SPWorley on 3/22/2009稍作修改的一个（详情见帖子）

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

旁注：一些用户指出，在某些情况下编译时，这可能会导致错误的答案。

【讨论】：

【解决方案9】：

如果您正在寻找 c++ 的答案并且您来到这里，因为它归结为数零，那么您会得到 std::countl_zero，根据 godbolt.org 的说法，它称为 bsr。 std::countl_zero 在 C++20 中可用，您可能需要将 -std=gnu++2a 添加到编译器命令行

【讨论】：

这个相同的 C++ 版本还添加了函数std::bit_width(x)，如果x > 0 和0 如果x == 0，则导致1 + ⎣log2(x)⎦。