## Monday, January 4, 2010

### 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(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(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).