Counting Integers with a Particular Binary Property

The Hamming weight or popcount of $n$ is the number of 1s in the binary representation of $n$.

Question: how many positive integers less than $N$ have a popcount divisible by 3?

The naive $O(N \log N)$ solution is to loop through each integer $n < N$ and count the number of bits that are set. This is slow.

Don't fret! There is an $O(\log^2 N)$ solution.

Suppose that $N$ is a power of 2. Then its binary representation is 1 followed by $\lg(N)$ 0s. We can count the number of positive $n\leq N$ with $\mathrm{popcount}(n)$ by simply doing $C(\lg(N), 3) + C(\lg(N), 6) + ... + C(\lg(N), 3k)$ while $3k \leq \lg(N)$. For instance, when $N=16$, $C(4,3) = 4$ which is correct.

Next, suppose that $\mathrm{popcount}(N) = 2$, (i.e. $N$ has exactly two bits set, i.e. N$ $is the sum of a pair of distinct powers of 2). E.g. say $N = 20$ (binary 10100). The solution that worked for $N=16$ will still work. Now in a sense we've "considered" that power of 2. Then consider the next power of 2 that makes up the number, 4. 4 has 2 0's following, so we use the same principle as before! However, we remember that the bit representing 16 has already been set, so we're no longer looking for popcounts that are multiples of 3; instead, we're looking for popcounts which give a remainder of 2 modulo 3. Mathematically, we do $C(\lg(4), 2) + ... + C(\lg(4), 3k-1)$ while $3k-1 \leq 4$. In this case, we have only $C(2,2) = 1$. This corresponds to $n = 19$ (binary 10011). So for $N=20$, the answer is $4+1=5$.

OK now we generalise the pattern.

from euler import popcount, binom
from numpy import log2

def solve(N):
    ans = 0
    num_considered = 0
    p = 2**60
    oldN = N
    while p >= 1:
        if N & p:
            if (num_considered + 1) % 3 == 0:
                ans += 1
            mult = 3 - (num_considered % 3)
            while mult <= log2(p):
                ans += binom(log2(p), mult)
                mult += 3
            num_considered += 1
        p //= 2
    return ans

I'm not even sure that this is $O(\log^2 N)$... there are two nested loops that each run $O(\log n)$ times, and inside the nested loop we call a binomial function which in theory is $O(nk)$, and we call it with $n,k \leq \log N$ producing an extra $O(\log^2 N)$ multiplier. But due to the ranges of the arguments we use for the binomial function, I'm pretty sure each call is amortized constant time. Maybe.

euler is my Project Euler helper functions module which is available on my Github.