Detailed mathematical approach

Problem Summary

We must evaluate

$S(N)=\sum_{n=1}^{N} f(n)^2 \pmod{10^9+7},$

where

$f(1)=1,\qquad f(2m)=2f(m),\qquad f(2m+1)=2m+1+2f(m)+\frac{f(m)}{m}\quad (m\ge 1).$

A direct recurrence evaluation up to \(N=10^{16}\) is impossible. The key is to identify a closed form for \(f(n)\), then sum that closed form with a binary prefix method instead of iterating over all integers.

Mathematical Approach

Write \(w(n)=\operatorname{popcount}(n)\), the number of ones in the binary expansion of \(n\). The whole solution is driven by the identity \(f(n)=n\,w(n)\).

Step 1: Prove the closed form for \(f(n)\)

The binary weight satisfies

$w(2m)=w(m),\qquad w(2m+1)=w(m)+1.$

Assume inductively that \(f(m)=m\,w(m)\). Then for even arguments,

$f(2m)=2f(m)=2m\,w(m)=2m\,w(2m).$

For odd arguments, the division term becomes harmless because \(f(m)/m=w(m)\):

$\begin{aligned} f(2m+1)&=2m+1+2f(m)+\frac{f(m)}{m}\\ &=2m+1+2m\,w(m)+w(m)\\ &=(2m+1)(w(m)+1)\\ &=(2m+1)w(2m+1). \end{aligned}$

Since \(f(1)=1=1\cdot w(1)\), induction yields

$f(n)=n\,w(n).$

Therefore the required sum is simply

$S(N)=\sum_{n=1}^{N} n^2 w(n)^2.$

Step 2: Precompute moments over all \(m\)-bit suffixes

For each \(m\ge 0\), let \(X_m=\{0,1,\dots,2^m-1\}\). For \(a,b\in\{0,1,2\}\), define

$T_m^{a,b}=\sum_{x\in X_m} x^a w(x)^b.$

These nine moments are exactly what we need, because any expansion of a square in the value and a square in the popcount can only produce powers \(x^a w(x)^b\) with \(a,b\le 2\).

The target quantity over a complete \(m\)-bit block is the special case

$T_m^{2,2}=\sum_{x=0}^{2^m-1} x^2 w(x)^2.$

The base case is immediate:

$T_0^{0,0}=1,\qquad T_0^{a,b}=0\ \text{for}\ (a,b)\ne(0,0),$

because the only \(0\)-bit suffix is \(x=0\).

Step 3: Extend from \(m\) bits to \(m+1\) bits

Every \((m+1)\)-bit number is either \(2x\) or \(2x+1\) with \(x\in X_m\). Their popcounts satisfy

$w(2x)=w(x),\qquad w(2x+1)=w(x)+1.$

So each new moment obeys

$T_{m+1}^{a,b}=\sum_{x\in X_m}(2x)^a w(x)^b+\sum_{x\in X_m}(2x+1)^a (w(x)+1)^b.$

Because \(a,b\le 2\), the right-hand side always reduces to a linear combination of the same nine moments. The smallest updates are

$T_{m+1}^{0,0}=2T_m^{0,0},\qquad T_{m+1}^{0,1}=2T_m^{0,1}+T_m^{0,0},$

while the target second-order moment expands to

$\begin{aligned} T_{m+1}^{2,2}={}&8T_m^{2,2}+8T_m^{2,1}+4T_m^{1,2}+8T_m^{1,1}\\ &+4T_m^{2,0}+4T_m^{1,0}+T_m^{0,2}+2T_m^{0,1}+T_m^{0,0}. \end{aligned}$

This is the reason the implementations maintain exactly nine aggregate tables and update them level by level.

Step 4: Evaluate a full suffix block in one formula

During the final summation, the higher bits of a number are fixed first. Suppose the already chosen prefix contributes a value \(p\) and contains \(r\) ones. If \(m\) lower bits are still free, every number in that block has the form

$n=p+x,\qquad 0\le x\lt 2^m,$

and its popcount is

$w(n)=r+w(x),$

because the suffix occupies positions strictly below the fixed prefix.

Hence the whole block contributes

$B(p,r,m)=\sum_{x=0}^{2^m-1}(p+x)^2(r+w(x))^2.$

Expanding both squares gives a fixed linear combination of the precomputed moments:

$\begin{aligned} B(p,r,m)={}&p^2r^2T_m^{0,0}+2p^2rT_m^{0,1}+p^2T_m^{0,2}\\ &+2pr^2T_m^{1,0}+4prT_m^{1,1}+2pT_m^{1,2}\\ &+r^2T_m^{2,0}+2rT_m^{2,1}+T_m^{2,2}. \end{aligned}$

So once the nine moments are known, an entire interval of length \(2^m\) is summed in constant time.

Step 5: Scan \(N\) bit by bit

Write \(N\) in binary and scan from the most significant bit down to the least significant bit. Whenever the current bit of \(N\) is \(1\), there is a complete block of smaller numbers obtained by keeping all earlier bits equal to the current prefix, setting the current bit to \(0\), and letting all lower bits vary freely. That full block is added with \(B(p,r,m)\).

After adding the block, the scan turns the current bit on in the prefix and increases the prefix popcount by one. When the loop ends, the prefix itself equals \(N\), so we add the single endpoint term

$N^2 w(N)^2.$

Every integer from \(0\) to \(N\) appears exactly once in this decomposition, and the term for \(0\) is automatically zero.

Step 6: Worked example with \(N=13\)

Since \(13=1101_2\), the scan splits the range into four parts:

$[0,7],\qquad [8,11],\qquad [12,12],\qquad \{13\}.$

These four pieces correspond to

$B(0,0,3),\qquad B(8,1,2),\qquad B(12,2,0),\qquad 13^2\cdot 3^2.$

Their values are

$742,\qquad 1877,\qquad 576,\qquad 1521,$

so

$S(13)=742+1877+576+1521=4716.$

As another small checkpoint, the same closed form gives

$S(10)=\sum_{n=1}^{10} n^2 w(n)^2=1530.$

How the Code Works

The C++, Python, and Java implementations follow the same pipeline. First they precompute the nine tables \(T_m^{a,b}\) for every bit length up to 64. In parallel, they precompute the powers of two modulo \(10^9+7\), so a fixed prefix value can be updated immediately when the scan accepts a new set bit.

Next they scan the bits of \(N\) from high to low. At each set bit, the implementation adds the complete suffix block determined by the current prefix, using the expanded formula for \(B(p,r,m)\). Then it extends the prefix by turning that bit on and increasing the prefix popcount. After all bits have been processed, it adds the single endpoint term \(N^2 w(N)^2\).

All arithmetic is performed modulo \(10^9+7\), and the algorithm never iterates through the integers \(1,2,\dots,N\) one by one.

Complexity Analysis

Let \(B=\lfloor\log_2 N\rfloor+1\). Building the nine moment tables up to bit length \(B\) takes \(O(B)\) time and \(O(B)\) memory, because each new level is obtained from the previous one by a constant amount of arithmetic. The prefix scan over \(N\) also costs \(O(B)\) time. For the actual input size here, the implementations use a fixed 64-bit bound, so both time and memory are effectively \(O(64)\).

Footnotes and References

1. Problem page: https://projecteuler.net/problem=759
2. Hamming weight / popcount: Wikipedia — Hamming weight
3. Binary number system: Wikipedia — Binary number
4. Moment (mathematics): Wikipedia — Moment (mathematics)
5. Dynamic programming: Wikipedia — Dynamic programming