计算连续整数上第一个设置位之和的最快方法?

Fastest way to compute sum of first set bit over consecutive integers?

编辑:我希望让我接受 2 个答案,因为缺一不可。我建议两者都读!

我正在尝试快速实现一个函数,该函数给定一个无符号 32 位整数 x returns 2^trailing_zeros(i)i=1..x-1 的总和,其中 trailing_zeroscount trailing zeros operation ,它被定义为在最低有效 1 位之后返回 0 位。这似乎是一种应该适用于巧妙的位操作实现的问题,无论输入如何,它都采用相同数量的指令,但我无法推导出它。

从数学上讲,2^trailing_zeros(i) 等于 2 的最大因数正好整除 i。因此,我们将 1..x-1.

的最大因素相加
i                   | 1    2     3    4    5    6    7    8    9    10
-----------------------------------------------------------------------
2^trailing_zeroes(i) | 1    2     1    4    1    2    1    8    1    2
-----------------------------------------------------------------------
Sum (desired value) | 0    1     3    4    8    9    11   12   20   21

如果我们 'plot' 值 -- 水平位置从左到右递增对应于 i 并且垂直位置从从上到下对应 trailing_zeroes(i).

 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
  2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2   2 
    4       4       4       4       4       4       4       4       4       4       4       4       4       4       4       4       4       4       4       4       4       4       4       4       4       4       4       4       4       4       4       4   
        8               8               8               8               8               8               8               8               8               8               8               8               8               8               8               8       
                16                               16                               16                               16                               16                               16                               16                               16               
                                32                                                               32                                                               32                                                               32                               
                                                                64                                                                                                                               64                                                               

在这里更容易看出 2 总是相隔 4 个,8 总是相隔 16 个等等的模式。但是,每个模式开始的时间不同 -- 8 直到 i=8 才开始, 16 直到 i=16 才开始,等等。如果您没有考虑到模式不会立即开始,您可能会想出不起作用的公式——例如,您可能会认为要确定进入总数的 8 的数量,您应该只计算 floor(x/16),但是 i=25 已经足够向右包括前两个 8

到目前为止我想到的最好的解决方案是:

对于每个幂的工作方式,它计算图上 x 和该幂的第一次出现之间的水平距离,例如x 和第一个 8 之间的距离是 (x-8),然后除以该幂的重复实例之间的距离,例如floor((x-8)/16),这给了我们该力量出现的次数,我们对该力量的总和,例如floor((x-8)/16)*8。然后我们添加给定功率的一个实例,因为该计算不包括该功率第一次出现的时间。

在实践中,这个实现应该非常快,因为 division/floor 可以通过向右移位来完成,而 2 的幂可以通过向左移位 1 来完成。然而,似乎仍然可以做得更好。对于更大的输入,此实现将循环更多,最多 32 次(它是 O(log2(n)),理想情况下我们希望 O(1) 没有巨大的查找 table 耗尽所有 CPU 缓存)。我一直在关注 BMI/BMI2 intrinsics,但我没有看到应用它们的明显方法。

虽然我的目标是用 C++ 或 Rust 等编译语言实现它,并具有真正的位移和内在函数,但我一直在 Python 中制作原型。下面包括我的脚本,其中包括我描述的实现,z(x),以及生成绘图的代码,tower(x)

#!/usr/bin/env python
# -*- coding: utf-8 -*-

from math import pow, floor, log, ceil

def leading_zeros(x):
    return len(bin(x).split('b')[-1].split('1')[-1])

def f(x):
    s = 0
    for c, i in enumerate(range(1,x)):
        a = pow(2, len(bin(i).split('b')[-1].split('1')[-1]))
        s += a
    return s

def g(x): return sum([pow(2,i)*floor((x+pow(2,i)-1)/pow(2,i+1)) for i in range(0,32)])

def h(x):
    s = 0
    extra = 0
    extra_s = 0
    for i in range(0,32):
        num = (x+pow(2,i)-1)
        den = pow(2,i+1)
        fraction = num/den
        floored = floor(num/den)
        power = pow(2,i)
        product = power*floored
        if product == 0:
            break
        s += product
        extra += (fraction - floored)
        extra_s += power*fraction
        #print(f"i={i} s={s} num={num} den={den} fraction={fraction} floored={floored} power={power} product={product} extra={extra} extra_s={extra_s}")
    return s

def z(x):
    upper_bound = floor(log(x,2)) if x > 0 else 0
    s = 0
    for i in range(upper_bound+1):
        num = (x - pow(2,i))
        den = pow(2,i+1)
        fraction = num/den
        floored = floor(fraction)
        added = pow(2,i)
        s += floored * added
        s += added
        print(f"i={i} s={s} upper_bound={upper_bound} num={num} den={den} floored={floored} added={added}")
    return s
#    return sum([floor((x - pow(2,i))/pow(2,i+1) + pow(2,i)) for i in range(floor(log(x, 2)))])

def tower(x):
    table = [[" " for i in range(x)] for j in range(ceil(log(x,2)))]
    for i in range(1,x):
        p = leading_zeros(i)
        table[p][i] = 2**p
    for row in table:
        for col in row:
            print(col,end='')
        print()


# h(9000)
for i in range(1,16):
    tower(i)
    print((i, f(i), g(i), h(i), z(i-1)))

观察如果我们从 1 数到 x 而不是 x−1,我们有一个模式:

x sum sum/x
1 1 1
2 3 1.5
4 8 2
8 20 2.5
16 48 3

因此我们可以很容易地计算出任何两个 p 的幂的总和为 p • (1 + ½ b),其中 b 是幂(等效于设置的位数或 log2力量)。通过归纳我们可以看出:如果1到2b的和是2b•(1+½b)(对于b=0),然后从1求和到 2b+1 重复两次单独的术语贡献,除了最后一个术语添加 2b+1而不是2b,所以和是2•2b•(1+½b) − 2b + 2b+1 = 2b+1•(1+½b) + ½•2b+1 = 2 b+1•(1+½(b+1)).

此外,在任意两个 2 的幂之间,较低的位重复前面的部分和。因此,对于任何 x,我们可以通过对其中的设置位求和来计算尾随零的累积数量。回想一下,这提供了从 1 到 x 的数字的总和,我们调整以得到从 1 到 x−1 的所需总和,从 x 计算前:

unsigned CountCumulative(unsigned x)
{
    --x;
    unsigned sum = 0;
    for (unsigned bit = 0; bit < sizeof x * CHAR_BIT; ++bit)
        sum += (x & 1u << bit) * (1 + bit * .5);
    return sum;
}

我们可以在x耗尽时终止循环:

unsigned CountCumulative(unsigned x)
{
    --x;
    unsigned sum = 0;
    for (unsigned bit = 0; x; ++bit, x >>= 1)
        sum += ((x & 1) << bit) * (1 + bit * .5);
    return sum;
}

正如 harold 指出的那样,我们可以分解出 1,因为 x 的每一位的值相加等于 x:

unsigned CountCumulative(unsigned x)
{
    --x;
    unsigned sum = x;
    for (unsigned bit = 0; x; ++bit, x >>= 1)
        sum += ((x & 1) << bit) * bit * .5;
    return sum;
}

然后消除浮点数:

unsigned CountCumulative(unsigned x)
{
    unsigned sum = --x;
    for (unsigned bit = 0; x; ++bit, x >>= 1)
        sum += ((x & 1) << bit) / 2 * bit;
    return sum;
}

请注意,当 bit 为零时,((x & 1) << bit) / 2 将丢失分数,但这无关紧要,因为 * bit 无论如何都会使贡献为零。对于 bit 的所有其他值,(x & 1) << bit 是偶数,因此除法不会丢失任何内容。

这会在某些时候溢出 unsigned,因此可能需要使用更宽的类型进行计算。

更多代码高尔夫

另一种根据位位置重复添加 x 位值一半的方法是移位 x(将其位值减半),然后重复添加,同时删除连续的位数从低到高:

unsigned CountCumulative(unsigned x)
{
    unsigned sum = --x;
    for (unsigned bit = 0; x >>= 1; ++bit)
        sum += x << bit;
    return sum;
}

基于 Eric Postpischil 的方法,这里有一种不用循环的方法。

请注意,每一位都乘以其位置,并将结果相加(有点,除了其中还有一个0.5的因子,让我们暂时搁置)。让我们称这些加起来的值为“部分产品”只是为了称呼它们,这样称呼它们并不准确,我想不出更好的东西。如果我们转置那么一点点,那么它是这样构建的:每个部分积的最低位是每个位乘以该位的位置的最低位。一位积是按位与,位置最低位的值为0,1,0,1等,所以得到x & 0xAAAAAAAA,每个部分积的第二位是x & 0xCCCCCCCC(并且“权重”为 2,因此必须乘以 2)等等

然后整个事情需要右移1,以考虑0.5的因素

总计:

unsigned CountCumulativeTrailingZeros(unsigned x)
{
    --x;
    unsigned sum = x;
    sum += (x >> 1) & 0x55555555;
    sum += x & 0xCCCCCCCC;
    sum += (x & 0xF0F0F0F0) << 1;
    sum += (x & 0xFF00FF00) << 2;
    sum += (x & 0xFFFF0000) << 3;
    return sum;
}

为了进一步说明,这里有一个更形象的例子。让我们再次暂时降低 0.5 的系数,它不会从根本上改变算法,但会增加一些复杂性。

首先,我在 v(一些示例值)的每一位上方写下该位在二进制中的位置(p0 是该位置的最低有效位,p1第二位等)。垂直阅读ps ,每一列都是一个数字:

p0: 10101010101010101010101010101010
p1: 11001100110011001100110011001100
p2: 11110000111100001111000011110000
p3: 11111111000000001111111100000000
p4: 11111111111111110000000000000000
v : 00000000100001000000001000000000

因此,例如第 9 位已设置,并且其上方(从下到上读取)为 01001(二进制为 9)。

我们想要做的(Eric 的回答已经解释了为什么这有效),是获取设置的位的索引,将它们移动到相应的位置,然后添加它们。在这种情况下,他们已经他们自己的位置(通过构造,数字写在他们自己的位置),所以没有移位,但他们仍然需要过滤所以只对应于设置位的数字继续存在。这就是我所说的“单位乘积”的意思:取 v 的一位并将其乘以 p0p1 等的相应位

您也可以将其视为将位值​​乘以其索引,因此 2^bit * bit 如评论中所述。这不是如何在这里完成的,但实际上是什么完成了。

回到示例,在这些部分产品中应用按位与结果:

pp0: 00000000100000000000001000000000
pp1: 00000000100001000000000000000000
pp2: 00000000100000000000000000000000
pp3: 00000000000000000000001000000000
pp4: 00000000100001000000000000000000
v  : 00000000100001000000001000000000

唯一剩下的值是01001、10010、10111,它们在它们对应的位置(所以,已经转移到它们需要去的地方)。

必须添加这些值,同时将它们保持在原位。它们不需要从它们所处的奇怪形式中提取出来,加法可以自由重新排序(结合和交换)所以可以先将部分产品的所有最低有效位添加到总和,然后是所有秒位, 等等。但是他们必须添加正确的“权重”,毕竟 pp0 中的设置位对应于该位置的 1 但 pp1 中的设置位实际上对应于该位置的 2 (因为它是它所属的数字的第二位)。所以直接用pp0,但是pp1左移1,pp2左移2等

仍然必须考虑 0.5 的因素,我主要是通过将部分乘积的位移动比它们的权重所暗示的少一位。 pp0 被左移了 0,所以现在必须 移 1。这可以通过将 return sum >> 1; 放在末尾来以较少的复杂性来完成,但这会减少函数在 运行 之前可以处理的值的范围进入整数包装模 232(这也会花费额外的操作,而以奇怪的方式进行操作则不会)。