找出正数范围内二进制表示中 1 的个数的算法

Question

我刚刚遇到一个问题，我们应该计算一个大范围内数字的二进制表示形式中 1 的个数。是否有任何算法或技术可以轻松找到它？例如， 对于输入 N = 6，其前面数字的二进制表示形式中 1 的数量。喜欢，1 - 0001 - 1 的数量 = 1;2 - 0010 - 1 的数量 = 1;3 - 0011 - 1 的数量 = 2;4 - 0100 - 1 的数量 = 1;
5 - 0101 - 1 的数量 = 2;

约束:1 <= N <= 10 ^ 20

所以总数是 7(1+1+2+1+2)。还有其他技巧可以找出这个吗？提前致谢!

Answer 1

设 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)} 中设置的位数是 S(n) 中设置的位数加上每个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，其中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

双返回值是动态规划的一种形式，避免多次重新计算相同的值。

第二个返回值是你感兴趣的。例如:

>> print(B(7)[1])
9

>> print(B(28)[1])
64

>> print(B(10**20)[1])
3301678091638143975424

这显然在 O(log N) 算术运算中运行，并使用 O(log N) 堆栈。

达到常量空间复杂度

只要小心一点，就可以将空间复杂度降低到 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] 是其中之一出现在 Ralf Stephan 方程的矩阵乘向量版本中的两个 4x4 矩阵，取决于 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) 算术运算，但使用常量空间。

我们可以做得更好一点，注意 M 矩阵总是以下形式:

[ a   b   c   d   ]
[ a-1 b+1 c   e   ]
[ 0   0   a+b a-1 ]
[ 0   0   0   1   ]

然后，B3 根据观察到的 M 结构，以优化的方式执行 B2 的矩阵乘法:

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

这与这种方法所能带给我们的一样好:唯一的算术运算是加法、乘以二、除以二和测试最低位。空间复杂度是常数。即使对于巨大的 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;
}