Problem 759: A Squared Recurrence Relation

View on Project Euler

Project Euler Problem 759 Solution

EulerSolve provides an optimized solution for Project Euler Problem 759, A Squared Recurrence Relation, with C++, Python, Java, and a step-by-step mathematical explanation.

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\}\)....

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

Problem 759 source code

C++

#include <array>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <vector>

namespace {

using i64 = long long;
using i128 = __int128_t;
using u64 = std::uint64_t;

constexpr i64 MOD = 1'000'000'007LL;
constexpr int MAX_BITS = 64;

struct Aggregates {
    std::array<i64, MAX_BITS + 1> cnt{};
    std::array<i64, MAX_BITS + 1> sx{};
    std::array<i64, MAX_BITS + 1> sx2{};
    std::array<i64, MAX_BITS + 1> sc{};
    std::array<i64, MAX_BITS + 1> sc2{};
    std::array<i64, MAX_BITS + 1> sxc{};
    std::array<i64, MAX_BITS + 1> sxc2{};
    std::array<i64, MAX_BITS + 1> sx2c{};
    std::array<i64, MAX_BITS + 1> sx2c2{};
};

i64 add_mod(i64 a, i64 b) {
    i64 r = a + b;
    if (r >= MOD) {
        r -= MOD;
    }
    return r;
}

i64 sub_mod(i64 a, i64 b) {
    i64 r = a - b;
    if (r < 0) {
        r += MOD;
    }
    return r;
}

i64 mul_mod(i64 a, i64 b) {
    return static_cast<i64>((static_cast<i128>(a) * b) % MOD);
}

Aggregates build_aggregates() {
    Aggregates ag;
    ag.cnt[0] = 1;

    for (int m = 0; m < MAX_BITS; ++m) {
        const i64 n = ag.cnt[m];
        const i64 sx = ag.sx[m];
        const i64 sx2 = ag.sx2[m];
        const i64 sc = ag.sc[m];
        const i64 sc2 = ag.sc2[m];
        const i64 sxc = ag.sxc[m];
        const i64 sxc2 = ag.sxc2[m];
        const i64 sx2c = ag.sx2c[m];
        const i64 sx2c2 = ag.sx2c2[m];

        const i64 two_sx = mul_mod(2, sx);
        const i64 four_sx = mul_mod(4, sx);
        const i64 four_sx2 = mul_mod(4, sx2);
        const i64 two_sc = mul_mod(2, sc);

        const i64 e_cnt = n;
        const i64 e_sx = two_sx;
        const i64 e_sx2 = four_sx2;
        const i64 e_sc = sc;
        const i64 e_sc2 = sc2;
        const i64 e_sxc = mul_mod(2, sxc);
        const i64 e_sxc2 = mul_mod(2, sxc2);
        const i64 e_sx2c = mul_mod(4, sx2c);
        const i64 e_sx2c2 = mul_mod(4, sx2c2);

        const i64 o_cnt = n;
        const i64 o_sx = add_mod(two_sx, n);
        const i64 o_sx2 = add_mod(add_mod(four_sx2, four_sx), n);
        const i64 o_sc = add_mod(sc, n);
        const i64 o_sc2 = add_mod(add_mod(sc2, two_sc), n);

        i64 o_sxc = 0;
        o_sxc = add_mod(o_sxc, mul_mod(2, sxc));
        o_sxc = add_mod(o_sxc, sc);
        o_sxc = add_mod(o_sxc, two_sx);
        o_sxc = add_mod(o_sxc, n);

        i64 o_sxc2 = 0;
        o_sxc2 = add_mod(o_sxc2, mul_mod(2, sxc2));
        o_sxc2 = add_mod(o_sxc2, sc2);
        o_sxc2 = add_mod(o_sxc2, mul_mod(4, sxc));
        o_sxc2 = add_mod(o_sxc2, mul_mod(2, sc));
        o_sxc2 = add_mod(o_sxc2, two_sx);
        o_sxc2 = add_mod(o_sxc2, n);

        i64 o_sx2c = 0;
        o_sx2c = add_mod(o_sx2c, mul_mod(4, sx2c));
        o_sx2c = add_mod(o_sx2c, mul_mod(4, sxc));
        o_sx2c = add_mod(o_sx2c, sc);
        o_sx2c = add_mod(o_sx2c, four_sx2);
        o_sx2c = add_mod(o_sx2c, four_sx);
        o_sx2c = add_mod(o_sx2c, n);

        i64 o_sx2c2 = 0;
        o_sx2c2 = add_mod(o_sx2c2, mul_mod(4, sx2c2));
        o_sx2c2 = add_mod(o_sx2c2, mul_mod(4, sxc2));
        o_sx2c2 = add_mod(o_sx2c2, sc2);
        o_sx2c2 = add_mod(o_sx2c2, mul_mod(8, sx2c));
        o_sx2c2 = add_mod(o_sx2c2, mul_mod(8, sxc));
        o_sx2c2 = add_mod(o_sx2c2, mul_mod(2, sc));
        o_sx2c2 = add_mod(o_sx2c2, four_sx2);
        o_sx2c2 = add_mod(o_sx2c2, four_sx);
        o_sx2c2 = add_mod(o_sx2c2, n);

        ag.cnt[m + 1] = add_mod(e_cnt, o_cnt);
        ag.sx[m + 1] = add_mod(e_sx, o_sx);
        ag.sx2[m + 1] = add_mod(e_sx2, o_sx2);
        ag.sc[m + 1] = add_mod(e_sc, o_sc);
        ag.sc2[m + 1] = add_mod(e_sc2, o_sc2);
        ag.sxc[m + 1] = add_mod(e_sxc, o_sxc);
        ag.sxc2[m + 1] = add_mod(e_sxc2, o_sxc2);
        ag.sx2c[m + 1] = add_mod(e_sx2c, o_sx2c);
        ag.sx2c2[m + 1] = add_mod(e_sx2c2, o_sx2c2);
    }

    return ag;
}

i64 block_sum(i64 prefix_value_mod, i64 prefix_popcount, int lower_bits, const Aggregates& ag) {
    const i64 p = prefix_value_mod;
    const i64 c = prefix_popcount % MOD;
    const i64 p2 = mul_mod(p, p);
    const i64 c2 = mul_mod(c, c);

    i64 ans = 0;
    ans = add_mod(ans, mul_mod(mul_mod(ag.cnt[lower_bits], p2), c2));
    ans = add_mod(ans, mul_mod(ag.sc[lower_bits], mul_mod(mul_mod(2, p2), c)));
    ans = add_mod(ans, mul_mod(ag.sc2[lower_bits], p2));
    ans = add_mod(ans, mul_mod(ag.sx[lower_bits], mul_mod(mul_mod(2, p), c2)));
    ans = add_mod(ans, mul_mod(ag.sxc[lower_bits], mul_mod(mul_mod(4, p), c)));
    ans = add_mod(ans, mul_mod(ag.sxc2[lower_bits], mul_mod(2, p)));
    ans = add_mod(ans, mul_mod(ag.sx2[lower_bits], c2));
    ans = add_mod(ans, mul_mod(ag.sx2c[lower_bits], mul_mod(2, c)));
    ans = add_mod(ans, ag.sx2c2[lower_bits]);
    return ans;
}

i64 S(u64 n, const Aggregates& ag, const std::array<i64, MAX_BITS + 1>& pow2_mod) {
    i64 ans = 0;
    i64 prefix_value_mod = 0;
    i64 prefix_popcount = 0;

    for (int bit = MAX_BITS - 1; bit >= 0; --bit) {
        if (((n >> bit) & 1ULL) == 0ULL) {
            continue;
        }
        ans = add_mod(ans, block_sum(prefix_value_mod, prefix_popcount, bit, ag));
        prefix_value_mod = add_mod(prefix_value_mod, pow2_mod[bit]);
        ++prefix_popcount;
    }

    const i64 p2 = mul_mod(prefix_value_mod, prefix_value_mod);
    const i64 c = prefix_popcount % MOD;
    const i64 c2 = mul_mod(c, c);
    ans = add_mod(ans, mul_mod(p2, c2));
    return ans;
}

u64 brute_S(int n) {
    std::vector<u64> f(static_cast<std::size_t>(n + 1), 0ULL);
    f[1] = 1ULL;
    for (int i = 2; i <= n; ++i) {
        if ((i & 1) == 0) {
            f[i] = 2ULL * f[i / 2];
        } else {
            int m = (i - 1) / 2;
            f[i] = static_cast<u64>(2 * m + 1) + 2ULL * f[m] + f[m] / static_cast<u64>(m);
        }
    }
    u64 sum = 0ULL;
    for (int i = 1; i <= n; ++i) {
        sum += f[i] * f[i];
    }
    return sum;
}

}  // namespace

int main() {
    const Aggregates ag = build_aggregates();

    std::array<i64, MAX_BITS + 1> pow2_mod{};
    pow2_mod[0] = 1;
    for (int i = 1; i <= MAX_BITS; ++i) {
        pow2_mod[i] = mul_mod(2, pow2_mod[i - 1]);
    }

    assert(brute_S(10) == 1530ULL);
    assert(brute_S(100) == 4'798'445ULL);
    assert(S(10ULL, ag, pow2_mod) == 1530LL);
    assert(S(100ULL, ag, pow2_mod) == 4'798'445LL);

    std::cout << S(10'000'000'000'000'000ULL, ag, pow2_mod) << '\n';
    return 0;
}

Python

MOD = 1000000007
MAX_BITS = 64

class Aggregates:
    def __init__(self):
        self.cnt = [0] * (MAX_BITS + 1)
        self.sx = [0] * (MAX_BITS + 1)
        self.sx2 = [0] * (MAX_BITS + 1)
        self.sc = [0] * (MAX_BITS + 1)
        self.sc2 = [0] * (MAX_BITS + 1)
        self.sxc = [0] * (MAX_BITS + 1)
        self.sxc2 = [0] * (MAX_BITS + 1)
        self.sx2c = [0] * (MAX_BITS + 1)
        self.sx2c2 = [0] * (MAX_BITS + 1)

def build_aggregates():
    ag = Aggregates()
    ag.cnt[0] = 1
    
    for m in range(MAX_BITS):
        n = ag.cnt[m]
        sx = ag.sx[m]
        sx2 = ag.sx2[m]
        sc = ag.sc[m]
        sc2 = ag.sc2[m]
        sxc = ag.sxc[m]
        sxc2 = ag.sxc2[m]
        sx2c = ag.sx2c[m]
        sx2c2 = ag.sx2c2[m]
        
        two_sx = (2 * sx) % MOD
        four_sx = (4 * sx) % MOD
        four_sx2 = (4 * sx2) % MOD
        two_sc = (2 * sc) % MOD
        
        e_cnt = n
        e_sx = two_sx
        e_sx2 = four_sx2
        e_sc = sc
        e_sc2 = sc2
        e_sxc = (2 * sxc) % MOD
        e_sxc2 = (2 * sxc2) % MOD
        e_sx2c = (4 * sx2c) % MOD
        e_sx2c2 = (4 * sx2c2) % MOD
        
        o_cnt = n
        o_sx = (two_sx + n) % MOD
        o_sx2 = (four_sx2 + four_sx + n) % MOD
        o_sc = (sc + n) % MOD
        o_sc2 = (sc2 + two_sc + n) % MOD
        
        o_sxc = (2 * sxc + sc + two_sx + n) % MOD
        o_sxc2 = (2 * sxc2 + sc2 + 4 * sxc + 2 * sc + two_sx + n) % MOD
        o_sx2c = (4 * sx2c + 4 * sxc + sc + four_sx2 + four_sx + n) % MOD
        o_sx2c2 = (4 * sx2c2 + 4 * sxc2 + sc2 + 8 * sx2c + 8 * sxc + 2 * sc + four_sx2 + four_sx + n) % MOD
        
        ag.cnt[m + 1] = (e_cnt + o_cnt) % MOD
        ag.sx[m + 1] = (e_sx + o_sx) % MOD
        ag.sx2[m + 1] = (e_sx2 + o_sx2) % MOD
        ag.sc[m + 1] = (e_sc + o_sc) % MOD
        ag.sc2[m + 1] = (e_sc2 + o_sc2) % MOD
        ag.sxc[m + 1] = (e_sxc + o_sxc) % MOD
        ag.sxc2[m + 1] = (e_sxc2 + o_sxc2) % MOD
        ag.sx2c[m + 1] = (e_sx2c + o_sx2c) % MOD
        ag.sx2c2[m + 1] = (e_sx2c2 + o_sx2c2) % MOD
        
    return ag

def block_sum(prefix_value_mod, prefix_popcount, lower_bits, ag):
    p = prefix_value_mod % MOD
    c = prefix_popcount % MOD
    p2 = (p * p) % MOD
    c2 = (c * c) % MOD
    
    ans = 0
    ans = (ans + ag.cnt[lower_bits] * p2 % MOD * c2) % MOD
    ans = (ans + ag.sc[lower_bits] * p2 % MOD * 2 % MOD * c) % MOD
    ans = (ans + ag.sc2[lower_bits] * p2) % MOD
    ans = (ans + ag.sx[lower_bits] * p % MOD * 2 % MOD * c2) % MOD
    ans = (ans + ag.sxc[lower_bits] * p % MOD * 4 % MOD * c) % MOD
    ans = (ans + ag.sxc2[lower_bits] * p % MOD * 2) % MOD
    ans = (ans + ag.sx2[lower_bits] * c2) % MOD
    ans = (ans + ag.sx2c[lower_bits] * 2 % MOD * c) % MOD
    ans = (ans + ag.sx2c2[lower_bits]) % MOD
    
    return ans

def S(n, ag, pow2_mod):
    ans = 0
    prefix_value_mod = 0
    prefix_popcount = 0
    
    for bit in range(MAX_BITS - 1, -1, -1):
        if ((n >> bit) & 1) == 0:
            continue
            
        ans = (ans + block_sum(prefix_value_mod, prefix_popcount, bit, ag)) % MOD
        prefix_value_mod = (prefix_value_mod + pow2_mod[bit]) % MOD
        prefix_popcount += 1
        
    p2 = (prefix_value_mod * prefix_value_mod) % MOD
    c = prefix_popcount % MOD
    c2 = (c * c) % MOD
    ans = (ans + p2 * c2) % MOD
    
    return ans

def solve():
    ag = build_aggregates()
    pow2_mod = [0] * (MAX_BITS + 1)
    pow2_mod[0] = 1
    for i in range(1, MAX_BITS + 1):
        pow2_mod[i] = (pow2_mod[i-1] * 2) % MOD
        
    return str(S(10000000000000000, ag, pow2_mod))

if __name__ == "__main__":
    print(solve())

Java

public class Euler759 {
    static final long MOD = 1000000007L;
    static final int MAX_BITS = 64;

    static class Aggregates {
        long[] cnt = new long[MAX_BITS + 1];
        long[] sx = new long[MAX_BITS + 1];
        long[] sx2 = new long[MAX_BITS + 1];
        long[] sc = new long[MAX_BITS + 1];
        long[] sc2 = new long[MAX_BITS + 1];
        long[] sxc = new long[MAX_BITS + 1];
        long[] sxc2 = new long[MAX_BITS + 1];
        long[] sx2c = new long[MAX_BITS + 1];
        long[] sx2c2 = new long[MAX_BITS + 1];
    }

    static long addMod(long a, long b) {
        long r = a + b;
        return (r >= MOD) ? (r - MOD) : r;
    }

    static long mulMod(long a, long b) {
        return (a * b) % MOD;
    }

    static Aggregates buildAggregates() {
        Aggregates ag = new Aggregates();
        ag.cnt[0] = 1;

        for (int m = 0; m < MAX_BITS; ++m) {
            long n = ag.cnt[m];
            long sx = ag.sx[m];
            long sx2 = ag.sx2[m];
            long sc = ag.sc[m];
            long sc2 = ag.sc2[m];
            long sxc = ag.sxc[m];
            long sxc2 = ag.sxc2[m];
            long sx2c = ag.sx2c[m];
            long sx2c2 = ag.sx2c2[m];

            long twoSx = mulMod(2, sx);
            long fourSx = mulMod(4, sx);
            long fourSx2 = mulMod(4, sx2);
            long twoSc = mulMod(2, sc);

            long eCnt = n;
            long eSx = twoSx;
            long eSx2 = fourSx2;
            long eSc = sc;
            long eSc2 = sc2;
            long eSxc = mulMod(2, sxc);
            long eSxc2 = mulMod(2, sxc2);
            long eSx2c = mulMod(4, sx2c);
            long eSx2c2 = mulMod(4, sx2c2);

            long oCnt = n;
            long oSx = addMod(twoSx, n);
            long oSx2 = addMod(addMod(fourSx2, fourSx), n);
            long oSc = addMod(sc, n);
            long oSc2 = addMod(addMod(sc2, twoSc), n);

            long oSxc = 0;
            oSxc = addMod(oSxc, mulMod(2, sxc));
            oSxc = addMod(oSxc, sc);
            oSxc = addMod(oSxc, twoSx);
            oSxc = addMod(oSxc, n);

            long oSxc2 = 0;
            oSxc2 = addMod(oSxc2, mulMod(2, sxc2));
            oSxc2 = addMod(oSxc2, sc2);
            oSxc2 = addMod(oSxc2, mulMod(4, sxc));
            oSxc2 = addMod(oSxc2, mulMod(2, sc));
            oSxc2 = addMod(oSxc2, twoSx);
            oSxc2 = addMod(oSxc2, n);

            long oSx2c = 0;
            oSx2c = addMod(oSx2c, mulMod(4, sx2c));
            oSx2c = addMod(oSx2c, mulMod(4, sxc));
            oSx2c = addMod(oSx2c, sc);
            oSx2c = addMod(oSx2c, fourSx2);
            oSx2c = addMod(oSx2c, fourSx);
            oSx2c = addMod(oSx2c, n);

            long oSx2c2 = 0;
            oSx2c2 = addMod(oSx2c2, mulMod(4, sx2c2));
            oSx2c2 = addMod(oSx2c2, mulMod(4, sxc2));
            oSx2c2 = addMod(oSx2c2, sc2);
            oSx2c2 = addMod(oSx2c2, mulMod(8, sx2c));
            oSx2c2 = addMod(oSx2c2, mulMod(8, sxc));
            oSx2c2 = addMod(oSx2c2, mulMod(2, sc));
            oSx2c2 = addMod(oSx2c2, fourSx2);
            oSx2c2 = addMod(oSx2c2, fourSx);
            oSx2c2 = addMod(oSx2c2, n);

            ag.cnt[m + 1] = addMod(eCnt, oCnt);
            ag.sx[m + 1] = addMod(eSx, oSx);
            ag.sx2[m + 1] = addMod(eSx2, oSx2);
            ag.sc[m + 1] = addMod(eSc, oSc);
            ag.sc2[m + 1] = addMod(eSc2, oSc2);
            ag.sxc[m + 1] = addMod(eSxc, oSxc);
            ag.sxc2[m + 1] = addMod(eSxc2, oSxc2);
            ag.sx2c[m + 1] = addMod(eSx2c, oSx2c);
            ag.sx2c2[m + 1] = addMod(eSx2c2, oSx2c2);
        }

        return ag;
    }

    static long blockSum(long prefixValueMod, long prefixPopcount, int lowerBits, Aggregates ag) {
        long p = prefixValueMod;
        long c = prefixPopcount % MOD;
        long p2 = mulMod(p, p);
        long c2 = mulMod(c, c);

        long ans = 0;
        ans = addMod(ans, mulMod(mulMod(ag.cnt[lowerBits], p2), c2));
        ans = addMod(ans, mulMod(ag.sc[lowerBits], mulMod(mulMod(2, p2), c)));
        ans = addMod(ans, mulMod(ag.sc2[lowerBits], p2));
        ans = addMod(ans, mulMod(ag.sx[lowerBits], mulMod(mulMod(2, p), c2)));
        ans = addMod(ans, mulMod(ag.sxc[lowerBits], mulMod(mulMod(4, p), c)));
        ans = addMod(ans, mulMod(ag.sxc2[lowerBits], mulMod(2, p)));
        ans = addMod(ans, mulMod(ag.sx2[lowerBits], c2));
        ans = addMod(ans, mulMod(ag.sx2c[lowerBits], mulMod(2, c)));
        ans = addMod(ans, ag.sx2c2[lowerBits]);
        return ans;
    }

    static long S(long n, Aggregates ag, long[] pow2Mod) {
        long ans = 0;
        long prefixValueMod = 0;
        long prefixPopcount = 0;

        for (int bit = MAX_BITS - 1; bit >= 0; --bit) {
            if (((n >> bit) & 1L) == 0L) {
                continue;
            }
            ans = addMod(ans, blockSum(prefixValueMod, prefixPopcount, bit, ag));
            prefixValueMod = addMod(prefixValueMod, pow2Mod[bit]);
            ++prefixPopcount;
        }

        long p2 = mulMod(prefixValueMod, prefixValueMod);
        long c = prefixPopcount % MOD;
        long c2 = mulMod(c, c);
        ans = addMod(ans, mulMod(p2, c2));
        return ans;
    }

    public static String solve() {
        Aggregates ag = buildAggregates();
        long[] pow2Mod = new long[MAX_BITS + 1];
        pow2Mod[0] = 1;
        for (int i = 1; i <= MAX_BITS; ++i) {
            pow2Mod[i] = mulMod(2, pow2Mod[i - 1]);
        }

        return Long.toString(S(10000000000000000L, ag, pow2Mod));
    }

    public static void main(String[] args) {
        System.out.println(solve());
    }
}