60

一个很棒的编程资源 Bit Twiddling Hacks 提出(在这里)以下方法来计算 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 位算法的外观有什么提示吗?

4

9 回答 9

86

内在函数非常快,但对于真正的跨平台、独立于编译器的 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 要求输入值大于零。

于 2012-07-09T15:55:28.340 回答
59

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

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

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

所以你可以定义:

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

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

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

工作示例

#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

于 2012-07-07T16:30:41.923 回答
19

我试图通过强制使用幻数将 O(lg(N)) 操作中的 Find the log base 2 of an N-bit integer 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];
}
于 2014-04-10T22:59:05.430 回答
11

以 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);at entry 添加到函数中。

以 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
}
于 2014-07-15T01:59:46.387 回答
5

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

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];
    }
于 2012-07-07T15:56:24.497 回答
5

如果您正在寻找 c++ 的答案并且您到了这里,并且由于归结为计数零,那么您std::countl_zero根据 godbolt.org 的调用得到了 which bsrstd::countl_zero可从 C++20 获得,您可能需要添加-std=gnu++2a到编译器命令行

于 2020-10-10T20:07:11.380 回答
2

该算法基本上找出哪个字节包含最高有效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 中位数的一半。我在这里使用移位而不是乘法(只是因为我喜欢它)。

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

于 2012-07-07T15:42:55.103 回答
1

拿着这个:

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为位,然后重新解释为位。IEEEfloat格式将 30-23 位中的指数存储为具有偏差 127 的整数,因此通过将其向右移动 23 位并减去偏差,我们得到 log2(x)。对于 64 位整数输入,x转换为double,其指数在 62-52 位(右移 52 位),指数偏差为 1023。

于 2019-07-05T15:08:10.473 回答
0

这是SPWorley 于 2009 年 3 月 22 日稍作修改的版本(详见帖子)

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

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

于 2019-10-20T15:19:35.720 回答