C语言 64位整数log2的快速计算

ss2ws0br  于 2023-02-21  发布在  其他
关注(0)|答案(9)|浏览(687)

一个很棒的编程资源,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位算法的提示?

vyu0f0g1

vyu0f0g11#

内部函数确实很快,但仍然不足以实现真正的跨平台、独立于编译器的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要求输入值大于零。

5ssjco0h

5ssjco0h2#

如果您使用的是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

eiee3dmh

eiee3dmh3#

我试图通过对幻数进行暴力处理,将带有乘法和查找的O(lg(N))操作中的N位整数的以2为底的对数转换为64位,不用说,这需要一段时间。
然后我找到了Desmond的答案,并决定尝试他的神奇数字作为起点。因为我有一个6核处理器,我并行运行它开始0x 07 EDD 5E 59 A4 E28 C2/ 6倍。我很惊讶它立即找到了一些东西。原来0x 07 EDD 5E 59 A4 E28 C2/ 2工作。
下面是0x 07 EDD 5E 59 A4 E28 C2的代码,它节省了移位和减法操作:

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

o2gm4chl4#

以2为底的整数对数

下面是我对64位无符号整数所做的,计算以2为底的对数的下限,它相当于最高有效位的下标,这个方法对大数来说是“非常快”的,因为它使用了一个展开的循环,总是以log 64 = 6个步骤执行。
本质上,它所做的是减去序列{ 0 ≤ k ≤ 5:2^(2^k)} = { 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
}
c90pui9n

c90pui9n5#

下面是一个非常紧凑 * 和 * 快速的扩展,没有使用额外的临时变量:

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

下面是一个用if链构建的示例,同样没有使用额外的临时变量,但可能不是最快的。

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

jaql4c8m6#

如果你在寻找c的答案,你会得到这里,因为它归结为计数零,那么你会得到std::countl_zero,根据godbolt.org被称为bsrstd::countl_zero可以从C20中获得,你可能需要在编译器命令行中添加-std=gnu++2a

hzbexzde

hzbexzde7#

该算法基本上找出哪个字节包含最高有效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中位数的一半。我在这里使用了移位而不是乘法(只是因为我喜欢这样)。
您可以向该算法添加一个查找表,以尽可能地加快它的速度,但它更多的是一个概念验证。

5hcedyr0

5hcedyr08#

拿着这个

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对于64位整数输入,x被强制转换为double,对于double,指数在位62 - 52中(向右移位52位)并且指数偏移为1023。

chy5wohz

chy5wohz9#

这是SPWorley on 3/22/2009的修改版本
支持〉55位,0返回0。要使其返回负数而不是0,请删除"|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;  
}

一个测试函数和输出...

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

相关问题