找出正数范围内二进制表示中 1 的个数的算法
Algorithm to find out the number of 1's in the binary representation over a range of positive numbers
我刚遇到一个问题，我们应该计算一个大范围内数字的二进制表示中 1 的个数。是否有任何算法或技术可以轻松找到它？
例如，
对于输入 N = 6，其前面数字的二进制表示形式中 1 的数量。喜欢，
1 - 0001 - No. of 1's = 1;
2 - 0010 - No. of 1's = 1;
3 - 0011 - No. of 1's = 2;
4 - 0100 - No. of 1's = 1;
5 - 0101 - No. of 1's = 2;
约束：1 <= N <= 10 ^ 20
所以总数是 7(1+1+2+1+2)。还有其他技巧可以找出这个吗？
提前致谢！
是的。我们先分析一下1和2的幂2k之间的个数（下界包括在内，上限*不包括）。后面我们会根据这个方法来解决一般问题。
那么这意味着最后的 k 位最终会选择所有位组合（000 除外，但这不包含任何设置位）。事实上，对于 k=3，我们看到 001、010、011、100、101、110 和 111。所以平均来说，一半的位被设置。因此我们知道设置的总位数是：
k
2
---
\ k k-1
/ --- = 2 * k
--- 2
i=0
因此对于 1（或 0 之间的范围，但这没有区别，因为 0没有设置位）和2k，我们有2k-1×k 设置位。例如 k=3，我们计算 22×3=12 位，这确实是什么当我们手动枚举它时我们会看到。
这对一般情况有何帮助？
假设我们要计算在0和l和2[=123之间设置的位数=]kk+1 ，那么我们可以先统计一下总共设置的位数2 k，然后将其与2k[=142之间设置的总位数相加=] 和 l.
现在后者当然还有一个问题：因为我们不知道如何计算。但是我们可以执行一个"shift"：我们可以计算出0和l-2k[=143=之间的总位数]（我们知道怎么做），并添加 l-2k 额外的结果。我们用同样的方法计算0和l-2k之间的总位数，我们但是知道 l-2k 中两个的最大幂将小于 2k, 因为 2k 是 l[=142= 中两个的最大幂]，所以保证"progress"。
如何将 l-2k 添加到结果中？举个例子：如果我们要计算000和110之间（不包括）的集合位数，那么我们就得把000、[=19的位数相加=]、010、011，也就是第一个"iteration"。第二次迭代然后是在 100 和 110 之间设置的位，因此我们通过执行移位并计算 00 和 10 之间的元素数来做到这一点，但是"original" 数字中的每个数字都设置了一个额外的位：最高设置位，因此我们计算迭代的元素数量，从而补偿位丢失。
算法：我们现在可以推导出一个算法：
def count_bit_range(n):
if n <= 1:
return 0
k = n.bit_length()-1
pk = 1 << k
pk1 = 1 << (k-1)
return k * pk1 + (n-pk) + count_bit_range(n-pk)
或非递归方法：
def count_bit_range(n):
c = 0
while n > 1:
k = n.bit_length()-1
pk = 1 << k
pk1 = 1 << (k-1)
c += k * pk1 + n - pk
n -= pk
return c
例如：
>>> count_bit_range(0)
0
>>> count_bit_range(1)
0
>>> count_bit_range(2)
1
>>> count_bit_range(3)
2
>>> count_bit_range(4)
4
>>> count_bit_range(5)
5
>>> count_bit_range(6)
7
>>> count_bit_range(12)
20
>>> count_bit_range(28)
64
例如对于 12，我们得到：
0001 0010 0011 0100 0101 0110 0111
1000 1001 1010 1011
所以 20 个设置位。
或 28:
00001 00010 00011 00100 00101 00110 00111
01000 01001 01010 01011 01100 01101 01110 01111
10000 10001 10010 10011 10100 10101 10110 10111
11000 11001 11010 11011
确实是 64。
Benchmarks：如果我们运行算法具有上限（1020)，我们在本地机器上获得 11.9 微秒：
>>> timeit(partial(count_bit_range, 10**20), number=1000000)
11.911393816000782
这（可能）不是范围内最昂贵的数字但是，递归调用的数量与上限的设置位数成比例，因此是最昂贵的范围内的昂贵数字可能是 (1<<66)-1:
>>> timeit(partial(count_bit_range, (1<<66)-1), number=1000000)
32.43066442897543
但是 32.4 微秒对于计算 1 到 73'786'976'294'838'206'463 之间设置的位数来说看起来仍然合理。
在本地机器上，它给出非递归方法的即时结果高达 1020'0000.
时间复杂度：递归调用的次数与上界中设置的位数成比例：实际上每次迭代都会删除最高设置位，当上限命中一或零。对于 w 位数，该算法因此需要 O(w) 次递归调用（这是 not 本身就是 "basic operations").
的数量
确切的时间复杂度有点难计算，因为在调用过程中，我们对变量执行大量计算，严格来说，这些变量可以变得任意大，加法、减法等取非常量任意大数的时间。
我们可以假设大多数操作将 运行 与该数字的位数成线性关系（如 .bit_length() 和二进制移位），但乘法需要的时间超过线性时间（在长度）。特别是因为 k 的长度与上界的长度成对数，所以它是 "tricky".
如果我们假设递归步骤在上界的 length 中花费二次时间（这可能是高估了），那么我们得到的时间复杂度为 O(w3)，或者对于最大范围n，时间复杂度为O (log3n).
为值 0..2^P-1 创建一个 table，其中 P = 8
byte[] table = new byte[] {0,1,1,2,1,2,1,3, ... 7,8};
和所有长度单位P的掩码：
long mask = (1 << P)-1;
然后，以字节为单位拆分输入数字，并对每个字节求和：
int numUnits(long number) {
int sum=0;
for (int k=0; k<64/P, k++) {
sum += table[number & mask];
num = num >> P;
}
return sum;
}
您可以取 P = 4 或 16 而不是 8，具体取决于您可以为 table 提供多少内存。
设 S(n) 为数字 0 到 n 的集合（没有重复，但顺序任意）。然后 S(2n+1) = {2*s for s in S(n)} + {2*s+1 for s in S(n)}，和 S(2n) = {2*s for s in S(n)} + {2*s+1 for s in S(n-1)}.
两个例子：
S(7) = {2*s for s in S(3)} + {2*s+1 for s in S(3)}
= {0, 2, 4, 6} + {1, 3, 5, 7}
S(10) = {2*s for s in S(5)} + {2*s+1 for s in S(4)}
= {0, 2, 4, 6, 8, 10} + {1, 3, 5, 7, 9}
让a(n)定义为S(n)中所有数字的总位数，并使用S的公式，我们有a(2n+1) = 2a(n) + n+1，和 a(2n) = a(n) + a(n-1) + n。这是因为{2*s for s in S(n)}中设置的位数与S(n)中设置的位数相同，而{2*s+1 for s in S(n)}中设置的位数与[=中设置的位数相同=24=]对S(n)的每个元素加一（即：n+1）。
这些相同的方程出现在 https://oeis.org/A000788 上，归功于 Ralf Stephan：
a(0) = 0
a(2n) = a(n)+a(n-1)+n
a(2n+1) = 2a(n)+n+1
使用这个，可以写一个函数 B with B(N) = a(N), a(N-1):
def B(N):
if N == 0:
return 0, 0
r, s = B(N//2)
if N % 2:
return 2*r+N//2+1, r+s+N//2
else:
return r+s+N//2, 2*s+N//2
双 return 值是动态规划的一种形式，避免多次重新计算相同的值。
第二个 return 值是您感兴趣的值。例如：
>> print(B(7)[1])
9
>> print(B(28)[1])
64
>> print(B(10**20)[1])
3301678091638143975424
这显然在 O(log N) 算术运算中运行，并使用 O(log N) 堆栈。
变得恒定 space 复杂性
只要小心一点，就可以将 space 的复杂度降低到 O(1)。
我们可以将 Ralf Stephan 方程写成矩阵乘以向量的形式：
[ a(2n+1) ] = [2 0 1 1] [ a(n) ]
[ a(2n) ] [1 1 1 0] * [ a(n-1)]
[ 2n+1 ] [0 0 2 1] [ n ]
[ 1 ] [0 0 0 1] [ 1 ]
和
[ a(2n) ] = [1 1 1 0] [ a(n) ]
[ a(2n-1) ] [0 2 1 0] * [ a(n-1)]
[ 2n ] [0 0 2 0] [ n ]
[ 1 ] [0 0 0 1] [ 1 ]
重复应用这些规则中的一个或另一个，得到：
[ a(n) ] = M[0] * M[1] * ... * M[k] * [ a(0) ]
[ a(n-1)] [ a(-1)]
[ n ] [ 0 ]
[ 1 ] [ 1 ]
其中 M[0]、M[1]、...、M[k] 是出现在矩阵乘向量版本中的两个 4x4 矩阵中的一个或另一个Ralf Stephan 方程，取决于 n.
的第 k 位
因此：
def mat_mul(A, B):
C = [[0] * 4 for _ in range(4)]
for i in range(4):
for j in range(4):
for k in range(4):
C[i][k] += A[i][j] * B[j][k]
return C
M1 = [[2, 0, 1, 1], [1, 1, 1, 0], [0, 0, 2, 1], [0, 0, 0, 1]]
M0 = [[1, 1, 1, 0], [0, 2, 1, 0], [0, 0, 2, 0], [0, 0, 0, 1]]
def B2(N):
M = [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]
while N:
M = mat_mul(M, M1 if N%2 else M0)
N >>= 1
return M[1][3]
函数B2执行O(log n)算术运算，但使用常量space。
我们可以做得更好一点，注意 M 矩阵总是以下形式：
[ a b c d ]
[ a-1 b+1 c e ]
[ 0 0 a+b a-1 ]
[ 0 0 0 1 ]
然后，B3 以优化的方式执行 B2 的矩阵乘法，具体取决于 M 的观察结构：
def B3(N):
a, b, c, d, e = 1, 0, 0, 0, 0
while N:
if N%2:
a, c, d, e = 2*a+b, a+b+2*c, a+c+d, a+c+e-1
else:
b, c = a+2*b, a+b+2*c
N >>= 1
return e
这和这种方法所能带给我们的一样好：唯一的算术运算是加法、乘以二、除以二和测试最低位。 space 复杂度是恒定的。即使对于巨大的 N（例如 10^200），所花费的时间也可以忽略不计。
C 中的快速版本
对于速度，C 版本（使用 gcc 的 __int128 扩展）在我的机器上计算 b3(10**20) 大约需要 140 纳秒。该代码是 B3 python 函数的直接转换（注意 d 不是必需的），由于 C 中缺少多重赋值而受到轻微阻碍。
typedef unsigned __int128 uint128;
uint128 b3(uint128 n) {
uint128 a=1, b=0, c=0, e=0;
while (n) {
if (n&1) {
e = a+c+e-1;
c = a+b+2*c;
a = 2*a+b;
} else {
c = a+b+2*c;
b = a+2*b;
}
n >>= 1;
}
return e;
}
我刚遇到一个问题，我们应该计算一个大范围内数字的二进制表示中 1 的个数。是否有任何算法或技术可以轻松找到它？
例如，
对于输入 N = 6，其前面数字的二进制表示形式中 1 的数量。喜欢，
1 - 0001 - No. of 1's = 1;
2 - 0010 - No. of 1's = 1;
3 - 0011 - No. of 1's = 2;
4 - 0100 - No. of 1's = 1;
5 - 0101 - No. of 1's = 2;
约束：1 <= N <= 10 ^ 20
所以总数是 7(1+1+2+1+2)。还有其他技巧可以找出这个吗？ 提前致谢！
是的。我们先分析一下1和2的幂2k之间的个数（下界包括在内，上限*不包括）。后面我们会根据这个方法来解决一般问题。
那么这意味着最后的 k 位最终会选择所有位组合（000 除外，但这不包含任何设置位）。事实上，对于 k=3，我们看到 001、010、011、100、101、110 和 111。所以平均来说，一半的位被设置。因此我们知道设置的总位数是：
k
2
---
\ k k-1
/ --- = 2 * k
--- 2
i=0
因此对于 1（或 0 之间的范围，但这没有区别，因为 0没有设置位）和2k，我们有2k-1×k 设置位。例如 k=3，我们计算 22×3=12 位，这确实是什么当我们手动枚举它时我们会看到。
这对一般情况有何帮助？
假设我们要计算在0和l和2[=123之间设置的位数=]k
现在后者当然还有一个问题：因为我们不知道如何计算。但是我们可以执行一个"shift"：我们可以计算出0和l-2k[=143=之间的总位数]（我们知道怎么做），并添加 l-2k 额外的结果。我们用同样的方法计算0和l-2k之间的总位数，我们但是知道 l-2k 中两个的最大幂将小于 2k, 因为 2k 是 l[=142= 中两个的最大幂]，所以保证"progress"。
如何将 l-2k 添加到结果中？举个例子：如果我们要计算000和110之间（不包括）的集合位数，那么我们就得把000、[=19的位数相加=]、010、011，也就是第一个"iteration"。第二次迭代然后是在 100 和 110 之间设置的位，因此我们通过执行移位并计算 00 和 10 之间的元素数来做到这一点，但是"original" 数字中的每个数字都设置了一个额外的位：最高设置位，因此我们计算迭代的元素数量，从而补偿位丢失。
算法：我们现在可以推导出一个算法：
def count_bit_range(n):
if n <= 1:
return 0
k = n.bit_length()-1
pk = 1 << k
pk1 = 1 << (k-1)
return k * pk1 + (n-pk) + count_bit_range(n-pk)
或非递归方法：
def count_bit_range(n):
c = 0
while n > 1:
k = n.bit_length()-1
pk = 1 << k
pk1 = 1 << (k-1)
c += k * pk1 + n - pk
n -= pk
return c
例如：
>>> count_bit_range(0)
0
>>> count_bit_range(1)
0
>>> count_bit_range(2)
1
>>> count_bit_range(3)
2
>>> count_bit_range(4)
4
>>> count_bit_range(5)
5
>>> count_bit_range(6)
7
>>> count_bit_range(12)
20
>>> count_bit_range(28)
64
例如对于 12，我们得到：
0001 0010 0011 0100 0101 0110 0111
1000 1001 1010 1011
所以 20 个设置位。
或 28:
00001 00010 00011 00100 00101 00110 00111
01000 01001 01010 01011 01100 01101 01110 01111
10000 10001 10010 10011 10100 10101 10110 10111
11000 11001 11010 11011
确实是 64。
Benchmarks：如果我们运行算法具有上限（1020)，我们在本地机器上获得 11.9 微秒：
>>> timeit(partial(count_bit_range, 10**20), number=1000000)
11.911393816000782
这（可能）不是范围内最昂贵的数字但是，递归调用的数量与上限的设置位数成比例，因此是最昂贵的范围内的昂贵数字可能是 (1<<66)-1:
>>> timeit(partial(count_bit_range, (1<<66)-1), number=1000000)
32.43066442897543
但是 32.4 微秒对于计算 1 到 73'786'976'294'838'206'463 之间设置的位数来说看起来仍然合理。
在本地机器上，它给出非递归方法的即时结果高达 1020'0000.
时间复杂度：递归调用的次数与上界中设置的位数成比例：实际上每次迭代都会删除最高设置位，当上限命中一或零。对于 w 位数，该算法因此需要 O(w) 次递归调用（这是 not 本身就是 "basic operations").
的数量确切的时间复杂度有点难计算，因为在调用过程中，我们对变量执行大量计算，严格来说，这些变量可以变得任意大，加法、减法等取非常量任意大数的时间。
我们可以假设大多数操作将 运行 与该数字的位数成线性关系（如 .bit_length() 和二进制移位），但乘法需要的时间超过线性时间（在长度）。特别是因为 k 的长度与上界的长度成对数，所以它是 "tricky".
如果我们假设递归步骤在上界的 length 中花费二次时间（这可能是高估了），那么我们得到的时间复杂度为 O(w3)，或者对于最大范围n，时间复杂度为O (log3n).
为值 0..2^P-1 创建一个 table，其中 P = 8
byte[] table = new byte[] {0,1,1,2,1,2,1,3, ... 7,8};
和所有长度单位P的掩码：
long mask = (1 << P)-1;
然后，以字节为单位拆分输入数字，并对每个字节求和：
int numUnits(long number) {
int sum=0;
for (int k=0; k<64/P, k++) {
sum += table[number & mask];
num = num >> P;
}
return sum;
}
您可以取 P = 4 或 16 而不是 8，具体取决于您可以为 table 提供多少内存。
设 S(n) 为数字 0 到 n 的集合（没有重复，但顺序任意）。然后 S(2n+1) = {2*s for s in S(n)} + {2*s+1 for s in S(n)}，和 S(2n) = {2*s for s in S(n)} + {2*s+1 for s in S(n-1)}.
两个例子：
S(7) = {2*s for s in S(3)} + {2*s+1 for s in S(3)}
= {0, 2, 4, 6} + {1, 3, 5, 7}
S(10) = {2*s for s in S(5)} + {2*s+1 for s in S(4)}
= {0, 2, 4, 6, 8, 10} + {1, 3, 5, 7, 9}
让a(n)定义为S(n)中所有数字的总位数，并使用S的公式，我们有a(2n+1) = 2a(n) + n+1，和 a(2n) = a(n) + a(n-1) + n。这是因为{2*s for s in S(n)}中设置的位数与S(n)中设置的位数相同，而{2*s+1 for s in S(n)}中设置的位数与[=中设置的位数相同=24=]对S(n)的每个元素加一（即：n+1）。
这些相同的方程出现在 https://oeis.org/A000788 上，归功于 Ralf Stephan：
a(0) = 0
a(2n) = a(n)+a(n-1)+n
a(2n+1) = 2a(n)+n+1
使用这个，可以写一个函数 B with B(N) = a(N), a(N-1):
def B(N):
if N == 0:
return 0, 0
r, s = B(N//2)
if N % 2:
return 2*r+N//2+1, r+s+N//2
else:
return r+s+N//2, 2*s+N//2
双 return 值是动态规划的一种形式，避免多次重新计算相同的值。
第二个 return 值是您感兴趣的值。例如：
>> print(B(7)[1])
9
>> print(B(28)[1])
64
>> print(B(10**20)[1])
3301678091638143975424
这显然在 O(log N) 算术运算中运行，并使用 O(log N) 堆栈。
变得恒定 space 复杂性
只要小心一点，就可以将 space 的复杂度降低到 O(1)。
我们可以将 Ralf Stephan 方程写成矩阵乘以向量的形式：
[ a(2n+1) ] = [2 0 1 1] [ a(n) ]
[ a(2n) ] [1 1 1 0] * [ a(n-1)]
[ 2n+1 ] [0 0 2 1] [ n ]
[ 1 ] [0 0 0 1] [ 1 ]
和
[ a(2n) ] = [1 1 1 0] [ a(n) ]
[ a(2n-1) ] [0 2 1 0] * [ a(n-1)]
[ 2n ] [0 0 2 0] [ n ]
[ 1 ] [0 0 0 1] [ 1 ]
重复应用这些规则中的一个或另一个，得到：
[ a(n) ] = M[0] * M[1] * ... * M[k] * [ a(0) ]
[ a(n-1)] [ a(-1)]
[ n ] [ 0 ]
[ 1 ] [ 1 ]
其中 M[0]、M[1]、...、M[k] 是出现在矩阵乘向量版本中的两个 4x4 矩阵中的一个或另一个Ralf Stephan 方程，取决于 n.
k 位
因此：
def mat_mul(A, B):
C = [[0] * 4 for _ in range(4)]
for i in range(4):
for j in range(4):
for k in range(4):
C[i][k] += A[i][j] * B[j][k]
return C
M1 = [[2, 0, 1, 1], [1, 1, 1, 0], [0, 0, 2, 1], [0, 0, 0, 1]]
M0 = [[1, 1, 1, 0], [0, 2, 1, 0], [0, 0, 2, 0], [0, 0, 0, 1]]
def B2(N):
M = [[1, 0, 0, 0], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1]]
while N:
M = mat_mul(M, M1 if N%2 else M0)
N >>= 1
return M[1][3]
函数B2执行O(log n)算术运算，但使用常量space。
我们可以做得更好一点，注意 M 矩阵总是以下形式：
[ a b c d ]
[ a-1 b+1 c e ]
[ 0 0 a+b a-1 ]
[ 0 0 0 1 ]
然后，B3 以优化的方式执行 B2 的矩阵乘法，具体取决于 M 的观察结构：
def B3(N):
a, b, c, d, e = 1, 0, 0, 0, 0
while N:
if N%2:
a, c, d, e = 2*a+b, a+b+2*c, a+c+d, a+c+e-1
else:
b, c = a+2*b, a+b+2*c
N >>= 1
return e
这和这种方法所能带给我们的一样好：唯一的算术运算是加法、乘以二、除以二和测试最低位。 space 复杂度是恒定的。即使对于巨大的 N（例如 10^200），所花费的时间也可以忽略不计。
C 中的快速版本
对于速度，C 版本（使用 gcc 的 __int128 扩展）在我的机器上计算 b3(10**20) 大约需要 140 纳秒。该代码是 B3 python 函数的直接转换（注意 d 不是必需的），由于 C 中缺少多重赋值而受到轻微阻碍。
typedef unsigned __int128 uint128;
uint128 b3(uint128 n) {
uint128 a=1, b=0, c=0, e=0;
while (n) {
if (n&1) {
e = a+c+e-1;
c = a+b+2*c;
a = 2*a+b;
} else {
c = a+b+2*c;
b = a+2*b;
}
n >>= 1;
}
return e;
}