Prime bits

OEIS A000788: Total number of 1's in binary expansions of 0, ..., n

One of the given formulas there is very interesting:

a(0) = 0;
a(2n) = a(n)+a(n-1)+n;
a(2n+1) = 2a(n)+n+1

I've still yet to decypher how that formula is derived, but here's mine in Maple:

full := n -> n * 2^n / 2:

flog2 := n -> floor(log[2](n)):

f := n -> `if`(n = 0, 0,
  full(flog2(n)) + (1 + n - 2^flog2(n)) + f(n - 2^flog2(n))
):

My formula basically partitions the rectangle of bits into 3 parts: top (first 2^k lines, explicit formula trivial), bottom-left (a single column of 1's) and bottom-right (the rest, counted recursively).

---

Ah, got it. The first formula reduces the problem by separating out the rightmost column and deinterleaving what's left. In any case, this problem is inspired by a Facebook puzzle, whose exact definition I can no longer find. I can assume it to be the following:

Given a, b, count how many numbers in the interval [a, b] have prime number of 1's in their binary representation.

My solution in Maple:

row := k -> add(binomial(k, i)*x^i, i=0..k):

g := (n, k) -> `if`(n = 0, 1,
  `if`(2^k > n, 
    g(n, k - 1),
    row(k) + x * g(n - 2^k, k - 1)
  )
):

f := n -> expand(sort(
  g(n, floor(evalf(log[2](n))))
)):

> f(100);
                   6       5       4       3       2
                2 x  + 11 x  + 26 x  + 33 x  + 21 x  + 7 x + 1

Basically, a term c*x^i in the polynomial f(n) means that among all the numbers in the interval [0..n], exactly c have i 1's in their binary representation.

The rest, of course, is straightforward:

frange := (lo, hi) -> f(hi) - f(lo - 1):

addprimecoeff := poly -> add(
  `if`(isprime(i), coeff(poly, x, i), 0),
  i=0..degree(poly)
):

> addprimecoeff(frange(2^50, 2^64));
                              4358342664758484397

Alternatively, we can follow the admittedly more elegant partitioning scheme given in the first formula and generate the polynomial as follows:

h := n -> `if`(n = 0, 1,
  `if`(type(n, even),
    x * h(n / 2 - 1) + h(n / 2),
    (1 + x) * h((n - 1) / 2)
  )
):

f := n -> expand(sort(
  h(n)
)):

> addprimecoeff(frange(999, 999999999));
                                   328122619

---

This agrees with somebody else's solution (no source code).