Problem 488: Unbalanced Nim
View on Project EulerProject Euler Problem 488 Solution
EulerSolve provides an optimized solution for Project Euler Problem 488, Unbalanced Nim, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary We consider three distinct heap sizes \(a\lt b\lt c\). A legal move decreases exactly one heap to a smaller nonnegative size, after which the three heap sizes are reordered and must still be distinct. Let \(\mathcal{P}\) denote the set of losing positions for this game. Define \(\mathcal{P}_N=\{(a,b,c)\in\mathcal{P}: 0\lt a\lt b\lt c\lt N\}\). The quantity to compute is $$F(N)=\sum_{(a,b,c)\in\mathcal{P}_N}(a+b+c)\pmod{10^9}.$$ A direct search over all triples below \(N\) is only feasible for tiny bounds. The implementation avoids that cubic search space by using a structural parametrization of all losing positions and then summing each block in closed form. Mathematical Approach The main game-specific fact used by the implementation is that every losing position belongs to a unique dyadic block indexed by $$p=2^k-1,\qquad k\ge 1,$$ and inside that block the losing triples are exactly $$c=p+n,\qquad b=p+t,\qquad a=(n\oplus t)-1,$$ with $$1\le n\le p,\qquad 0\le t\lt n.$$ Here \(\oplus\) is bitwise XOR. Because \(c=p+n\) lies in the interval \(p+1\le c\le 2p\), the largest heap chooses the block uniquely, so different values of \(k\) never overlap. Step 1: Restrict the Block to the Bound \(N\) To contribute to \(F(N)\), the largest heap must satisfy \(c\lt N\)....
Detailed mathematical approach
Problem Summary
We consider three distinct heap sizes \(a\lt b\lt c\). A legal move decreases exactly one heap to a smaller nonnegative size, after which the three heap sizes are reordered and must still be distinct. Let \(\mathcal{P}\) denote the set of losing positions for this game.
Define \(\mathcal{P}_N=\{(a,b,c)\in\mathcal{P}: 0\lt a\lt b\lt c\lt N\}\).
The quantity to compute is
$$F(N)=\sum_{(a,b,c)\in\mathcal{P}_N}(a+b+c)\pmod{10^9}.$$
A direct search over all triples below \(N\) is only feasible for tiny bounds. The implementation avoids that cubic search space by using a structural parametrization of all losing positions and then summing each block in closed form.
Mathematical Approach
The main game-specific fact used by the implementation is that every losing position belongs to a unique dyadic block indexed by
$$p=2^k-1,\qquad k\ge 1,$$
and inside that block the losing triples are exactly
$$c=p+n,\qquad b=p+t,\qquad a=(n\oplus t)-1,$$
with
$$1\le n\le p,\qquad 0\le t\lt n.$$
Here \(\oplus\) is bitwise XOR. Because \(c=p+n\) lies in the interval \(p+1\le c\le 2p\), the largest heap chooses the block uniquely, so different values of \(k\) never overlap.
Step 1: Restrict the Block to the Bound \(N\)
To contribute to \(F(N)\), the largest heap must satisfy \(c\lt N\). Since \(c=p+n\), this becomes
$$p+n\lt N\iff n\le N-1-p.$$
Therefore the relevant range of \(n\) is
$$1\le n\le m,\qquad m=\min\{p,\ N-1-p\}.$$
For each such \(n\), the parameter \(t\) runs through \(0,1,\dots,n-1\), which automatically ensures \(b\lt c\).
Step 2: Sum the Contribution of One Block
Fix a block parameter \(p\). For a fixed \(n\), summing over all allowed \(t\) gives
$$a+b+c=(n\oplus t)-1+(p+t)+(p+n).$$
Define the XOR row sum
$$X(n)=\sum_{t=0}^{n-1}(n\oplus t).$$
Then
$$\sum_{t=0}^{n-1}(a+b+c)=X(n)+\sum_{t=0}^{n-1}t+n^2+n(2p-1).$$
Because \(\sum_{t=0}^{n-1}t=\frac{n(n-1)}{2}\), the raw contribution of the block is
$$\sum_{n=1}^{m}\left(X(n)+\frac{n(n-1)}{2}+n^2+(2p-1)n\right).$$
So the whole problem reduces to four standard sums:
$$\sum_{n=1}^{m}X(n),\qquad \sum_{n=1}^{m}n,\qquad \sum_{n=1}^{m}n^2,\qquad \sum_{n=1}^{m}\frac{n(n-1)}{2}.$$
Step 3: Closed Forms for the Polynomial Terms
The ordinary arithmetic sums are classical:
$$\sum_{n=1}^{m}n=\frac{m(m+1)}{2},$$
$$\sum_{n=1}^{m}n^2=\frac{m(m+1)(2m+1)}{6},$$
$$\sum_{n=1}^{m}\frac{n(n-1)}{2}=\frac{m(m+1)(m-1)}{6}.$$
The implementation evaluates these formulas exactly by dividing by \(2\) and \(3\) before modular multiplication, so everything remains valid modulo \(10^9\).
Step 4: Fast Evaluation of the XOR Prefix Sum
Set \(X(0)=0\) and define
$$S(n)=\sum_{i=0}^{n-1}X(i).$$
Then the needed XOR contribution is simply \(\sum_{n=1}^{m}X(n)=S(m+1)\).
To evaluate \(S\) quickly, write
$$n=2^m+r,\qquad 0\le r\lt 2^m,$$
and let \(p=2^m\). For \(0\le r\lt p\), split the definition of \(X(p+r)\) into two ranges:
$$X(p+r)=\sum_{t=0}^{p-1}((p+r)\oplus t)+\sum_{u=0}^{r-1}((p+r)\oplus(p+u)).$$
In the first sum the top bit stays set, so \((p+r)\oplus t=p+(r\oplus t)\). As \(t\) runs through \(0,\dots,p-1\), the values \(r\oplus t\) are a permutation of \(0,\dots,p-1\). Hence
$$\sum_{t=0}^{p-1}((p+r)\oplus t)=p^2+\sum_{v=0}^{p-1}v=p^2+\frac{p(p-1)}{2}=\frac{p(3p-1)}{2}.$$
In the second sum the top bit cancels, so \((p+r)\oplus(p+u)=r\oplus u\). Therefore
$$X(p+r)=C_m+X(r),\qquad C_m=\frac{2^m(3\cdot 2^m-1)}{2}.$$
Summing this identity over \(r\) gives the key recursion
$$\boxed{S(2^m+r)=S(2^m)+r\,C_m+S(r).}$$
For powers of two we get
$$S(2^{m+1})=2S(2^m)+2^mC_m,$$
so the implementation precomputes \(S(2^m)\) and \(C_m\) once and then evaluates any \(S(n)\) in \(O(\log n)\) recursive steps.
Step 5: Correct the Auxiliary States with \(a=0\)
The parametrization above also produces states whose smallest heap is \(a=0\). These states appear naturally in the game graph, but they are not part of \(F(N)\), which only sums \(0\lt a\lt b\lt c\).
We have \(a=0\) exactly when
$$n\oplus t=1,$$
so the candidate is \(t=n\oplus 1\). This is inside the legal range \(0\le t\lt n\) if and only if \(n\) is odd, in which case \(t=n-1\). Thus each odd \(n\) contributes exactly one invalid triple.
For odd \(n\), the excluded sum equals
$$0+(p+n-1)+(p+n)=2p+2n-1.$$
If
$$q=\left\lfloor\frac{m+1}{2}\right\rfloor,$$
then the odd integers in \([1,m]\) are \(1,3,\dots,2q-1\), so
$$\sum_{\substack{1\le n\le m\\ 2\nmid n}}1=q,\qquad \sum_{\substack{1\le n\le m\\ 2\nmid n}}n=q^2.$$
Therefore the correction term for one block is
$$\boxed{(2p-1)q+2q^2.}$$
Step 6: Final Block Formula
Combining the raw block sum with the odd-\(n\) correction, the contribution of block \(p=2^k-1\) is
$$\begin{aligned} B(p,m)=&\ S(m+1)+\frac{m(m+1)(m-1)}{6}+\frac{m(m+1)(2m+1)}{6}\\ &+(2p-1)\frac{m(m+1)}{2}-\left((2p-1)q+2q^2\right), \end{aligned}$$
where
$$m=\min\{p,\ N-1-p\},\qquad q=\left\lfloor\frac{m+1}{2}\right\rfloor.$$
Hence the whole answer is
$$\boxed{F(N)=\sum_{\substack{k\ge 1\\ 2^k\lt N}} B(2^k-1,\ \min\{2^k-1,\ N-2^k\})\pmod{10^9}.}$$
After the structural theorem about losing positions, everything else is pure arithmetic.
Worked Example: \(F(8)=42\)
For \(N=8\), the only possible blocks are \(p=1\) and \(p=3\).
For \(p=1\), we have \(m=1\). The unique parametrized triple is \((0,1,2)\), so the entire block disappears after the \(a=0\) correction.
For \(p=3\), we have \(m=3\). First compute
$$X(1)=1,\qquad X(2)=2+3=5,\qquad X(3)=3+2+1=6,$$
hence
$$\sum_{n=1}^{3}X(n)=12,\qquad \sum_{n=1}^{3}n=6,\qquad \sum_{n=1}^{3}n^2=14,\qquad \sum_{n=1}^{3}\frac{n(n-1)}{2}=4.$$
The raw block total is
$$12+14+4+(2\cdot 3-1)\cdot 6=12+14+4+30=60.$$
Now \(q=\left\lfloor\frac{3+1}{2}\right\rfloor=2\), so the correction is
$$(2\cdot 3-1)\cdot 2+2\cdot 2^2=10+8=18.$$
Therefore
$$F(8)=60-18=42,$$
which matches the checkpoint used by the implementation. A second checkpoint is \(F(128)=496062\).
How the Implementation Works
The C++, Python, and Java implementations all follow the same plan. They enumerate the dyadic blocks \(p=2^k-1\) with \(p+1\lt N\), compute the cutoff \(m=\min\{p,N-1-p\}\), and add the closed-form block contribution modulo \(10^9\).
The only nontrivial subroutine is the evaluator for \(S(n)\). It uses the recursion above together with precomputed values at powers of two, so no block ever loops over all \(t\) or all \(n\). Because different blocks are independent, the outer accumulation can also be split across worker threads or processes.
Complexity Analysis
There are \(O(\log N)\) dyadic blocks because \(p=2^k-1\). Each block requires only a constant number of modular polynomial sums plus one evaluation of \(S(m+1)\), and the recursion for \(S\) has depth \(O(\log N)\). Therefore the total running time is
$$O((\log N)^2),$$
and the memory usage is \(O(\log N)\) for the precomputed tables of powers of two. This is exponentially smaller than any direct exploration of all triples below \(N\).
Footnotes and References
- Problem page: https://projecteuler.net/problem=488
- Nim and losing positions: Wikipedia — Nim
- XOR and bitwise algebra: Wikipedia — XOR
- Powers of two and binary intervals: Wikipedia — Power of two
- Berlekamp, Conway, Guy. Winning Ways for your Mathematical Plays.
Problem 488 source code
C++
#include <array>
#include <cstdint>
#include <cstdlib>
#include <iomanip>
#include <iostream>
#include <thread>
#include <vector>
using namespace std;
namespace {
constexpr uint64_t MOD = 1'000'000'000ULL;
inline uint64_t add_mod(uint64_t a, uint64_t b) {
a += b;
if (a >= MOD) a -= MOD;
return a;
}
inline uint64_t sub_mod(uint64_t a, uint64_t b) {
return (a >= b) ? (a - b) : (a + MOD - b);
}
inline uint64_t mul_mod(uint64_t a, uint64_t b) {
return static_cast<uint64_t>((__int128)(a % MOD) * (b % MOD) % MOD);
}
uint64_t sum_linear_mod(uint64_t m) {
// 1 + 2 + ... + m
if (m == 0) return 0;
uint64_t a = m, b = m + 1;
if ((a & 1ULL) == 0) {
a >>= 1;
} else {
b >>= 1;
}
return mul_mod(a, b);
}
uint64_t sum_square_mod(uint64_t m) {
// 1^2 + 2^2 + ... + m^2 = m(m+1)(2m+1)/6
if (m == 0) return 0;
uint64_t a = m, b = m + 1, c = 2 * m + 1;
// divide by 2
if ((a & 1ULL) == 0) {
a >>= 1;
} else if ((b & 1ULL) == 0) {
b >>= 1;
} else {
c >>= 1;
}
// divide by 3
if (a % 3ULL == 0) {
a /= 3ULL;
} else if (b % 3ULL == 0) {
b /= 3ULL;
} else {
c /= 3ULL;
}
return mul_mod(mul_mod(a, b), c);
}
uint64_t sum_choose2_mod(uint64_t m) {
// sum_{n=1..m} n(n-1)/2 = m(m+1)(m-1)/6
if (m <= 1) return 0;
uint64_t a = m, b = m + 1, c = m - 1;
// divide by 2
if ((a & 1ULL) == 0) {
a >>= 1;
} else if ((b & 1ULL) == 0) {
b >>= 1;
} else {
c >>= 1;
}
// divide by 3
if (a % 3ULL == 0) {
a /= 3ULL;
} else if (b % 3ULL == 0) {
b /= 3ULL;
} else {
c /= 3ULL;
}
return mul_mod(mul_mod(a, b), c);
}
struct Block {
uint64_t p;
uint64_t m;
};
class Solver488 {
public:
Solver488() {
precompute_prefix_xor_sums();
}
uint64_t solve_mod(uint64_t N, bool allow_multithreading = true) const {
vector<Block> blocks;
blocks.reserve(64);
for (int k = 1; k <= 62; ++k) {
uint64_t p = (1ULL << k) - 1ULL;
if (p + 1ULL >= N) break;
uint64_t m = min<uint64_t>(p, N - 1ULL - p);
blocks.push_back({p, m});
}
if (blocks.empty()) return 0;
unsigned workers = 1;
if (allow_multithreading) {
unsigned hw = thread::hardware_concurrency();
if (hw == 0) hw = 1;
workers = min<unsigned>(hw, static_cast<unsigned>(blocks.size()));
if (blocks.size() < 8) workers = 1;
}
if (workers == 1) {
uint64_t ans = 0;
for (const auto &blk : blocks) {
ans = add_mod(ans, block_sum_mod(blk.p, blk.m));
}
return ans;
}
vector<uint64_t> partial(workers, 0);
vector<thread> ts;
ts.reserve(workers);
for (unsigned tid = 0; tid < workers; ++tid) {
ts.emplace_back([&, tid]() {
size_t L = blocks.size() * tid / workers;
size_t R = blocks.size() * (tid + 1) / workers;
uint64_t local = 0;
for (size_t i = L; i < R; ++i) {
local = add_mod(local, block_sum_mod(blocks[i].p, blocks[i].m));
}
partial[tid] = local;
});
}
for (auto &th : ts) th.join();
uint64_t ans = 0;
for (uint64_t x : partial) ans = add_mod(ans, x);
return ans;
}
private:
// S(n) = sum_{i=0..n-1} X(i), X(i)=sum_{t=0..i-1}(i xor t)
array<uint64_t, 63> s_pow2_mod{};
array<uint64_t, 63> c_mod{};
void precompute_prefix_xor_sums() {
s_pow2_mod.fill(0);
c_mod.fill(0);
for (int m = 0; m <= 62; ++m) {
uint64_t p = (1ULL << m);
uint64_t x = p;
uint64_t y = 3ULL * p - 1ULL;
if ((x & 1ULL) == 0) {
x >>= 1;
} else {
y >>= 1;
}
c_mod[m] = mul_mod(x, y); // C_m = p(3p-1)/2
}
for (int m = 0; m < 62; ++m) {
uint64_t p = (1ULL << m) % MOD;
uint64_t term = mul_mod(p, c_mod[m]);
s_pow2_mod[m + 1] = add_mod(add_mod(s_pow2_mod[m], s_pow2_mod[m]), term);
}
}
uint64_t prefix_sum_x_mod(uint64_t n) const {
// returns S(n)=sum_{i=0..n-1}X(i)
if (n <= 1) return 0;
int msb = 63 - __builtin_clzll(n);
uint64_t p = (1ULL << msb);
if (n == p) return s_pow2_mod[msb];
uint64_t r = n - p;
uint64_t term = mul_mod(r % MOD, c_mod[msb]);
return add_mod(add_mod(s_pow2_mod[msb], term), prefix_sum_x_mod(r));
}
uint64_t block_sum_mod(uint64_t p, uint64_t m) const {
// P-positions in block p=2^k-1 are parameterized by:
// c=p+n, b=p+t, a=(n xor t)-1 with 1<=n<=m and 0<=t<n.
// Sum is then reduced to closed forms below.
uint64_t s_x = prefix_sum_x_mod(m + 1); // sum_{n=1..m} X(n)
uint64_t s1 = sum_linear_mod(m); // sum n
uint64_t s2 = sum_square_mod(m); // sum n^2
uint64_t sC2 = sum_choose2_mod(m); // sum n(n-1)/2
uint64_t two_p_minus_one = (2ULL * p - 1ULL) % MOD;
uint64_t total = 0;
total = add_mod(total, s_x);
total = add_mod(total, s2);
total = add_mod(total, sC2);
total = add_mod(total, mul_mod(two_p_minus_one, s1));
// Remove n odd term where t=n xor 1 = n-1 -> a=0 (not counted in F(N)).
uint64_t odd_cnt = (m + 1ULL) >> 1;
uint64_t odd_sq = mul_mod(odd_cnt % MOD, odd_cnt % MOD);
uint64_t correction = add_mod(mul_mod(odd_cnt % MOD, two_p_minus_one),
mul_mod(2ULL, odd_sq));
return sub_mod(total, correction);
}
};
uint64_t brute_F(uint32_t N) {
auto idx = [N](uint32_t a, uint32_t b, uint32_t c) {
return (a * N + b) * N + c;
};
vector<uint8_t> win(static_cast<size_t>(N) * N * N, 0);
uint64_t ans = 0;
for (uint32_t c = 2; c < N; ++c) {
for (uint32_t b = 1; b < c; ++b) {
for (uint32_t a = 0; a < b; ++a) {
bool w = false;
// Decrease a.
for (uint32_t x = 0; x < a && !w; ++x) {
if (!win[idx(x, b, c)]) w = true;
}
// Decrease b.
for (uint32_t x = 0; x < b && !w; ++x) {
if (x == a || x == c) continue;
uint32_t u, v, z;
if (x < a) {
u = x; v = a; z = c;
} else {
u = a; v = x; z = c;
}
if (!win[idx(u, v, z)]) w = true;
}
// Decrease c.
for (uint32_t x = 0; x < c && !w; ++x) {
if (x == a || x == b) continue;
uint32_t u, v, z;
if (x < b) {
if (x < a) {
u = x; v = a; z = b;
} else {
u = a; v = x; z = b;
}
} else {
u = a; v = b; z = x;
}
if (!win[idx(u, v, z)]) w = true;
}
win[idx(a, b, c)] = static_cast<uint8_t>(w);
if (!w && a > 0) {
ans += static_cast<uint64_t>(a) + b + c;
}
}
}
}
return ans;
}
void validate() {
Solver488 solver;
const uint64_t f8 = solver.solve_mod(8, false);
if (f8 != 42ULL) {
cerr << "Validation failed: F(8)=" << f8 << " (expected 42)\n";
exit(1);
}
const uint64_t f128 = solver.solve_mod(128, false);
if (f128 != 496062ULL) {
cerr << "Validation failed: F(128)=" << f128 << " (expected 496062)\n";
exit(1);
}
// Extra structural checks against brute force.
for (uint32_t n : {16u, 24u, 32u, 40u}) {
uint64_t exact_brute = brute_F(n);
uint64_t fast = solver.solve_mod(n, false);
if (exact_brute != fast) {
cerr << "Brute-force mismatch at N=" << n
<< ": brute=" << exact_brute
<< ", fast=" << fast << "\n";
exit(1);
}
}
}
} // namespace
int main() {
ios::sync_with_stdio(false);
cin.tie(nullptr);
validate();
Solver488 solver;
const uint64_t ans = solver.solve_mod(1'000'000'000'000'000'000ULL, true);
cout << setw(9) << setfill('0') << ans << "\n";
return 0;
}
Python
import multiprocessing
MOD = 1000000000
def add_mod(a, b):
return (a + b) % MOD
def sub_mod(a, b):
return (a - b) % MOD if a >= b else (a + MOD - b) % MOD
def mul_mod(a, b):
return (a * b) % MOD
def sum_linear_mod(m):
if m == 0: return 0
a, b = m, m + 1
if a % 2 == 0:
a //= 2
else:
b //= 2
return mul_mod(a, b)
def sum_square_mod(m):
if m == 0: return 0
a, b, c = m, m + 1, 2 * m + 1
if a % 2 == 0: a //= 2
elif b % 2 == 0: b //= 2
else: c //= 2
if a % 3 == 0: a //= 3
elif b % 3 == 0: b //= 3
else: c //= 3
return mul_mod(mul_mod(a, b), c)
def sum_choose2_mod(m):
if m <= 1: return 0
a, b, c = m, m + 1, m - 1
if a % 2 == 0: a //= 2
elif b % 2 == 0: b //= 2
else: c //= 2
if a % 3 == 0: a //= 3
elif b % 3 == 0: b //= 3
else: c //= 3
return mul_mod(mul_mod(a, b), c)
class Solver488:
def __init__(self):
self.s_pow2_mod = [0] * 63
self.c_mod = [0] * 63
self.precompute_prefix_xor_sums()
def precompute_prefix_xor_sums(self):
for m in range(63):
p = 1 << m
x, y = p, 3 * p - 1
if x % 2 == 0: x //= 2
else: y //= 2
self.c_mod[m] = mul_mod(x % MOD, y % MOD)
for m in range(62):
p = (1 << m) % MOD
term = mul_mod(p, self.c_mod[m])
self.s_pow2_mod[m + 1] = add_mod(add_mod(self.s_pow2_mod[m], self.s_pow2_mod[m]), term)
def prefix_sum_x_mod(self, n):
if n <= 1: return 0
msb = n.bit_length() - 1
p = 1 << msb
if n == p: return self.s_pow2_mod[msb]
r = n - p
term = mul_mod(r % MOD, self.c_mod[msb])
return add_mod(add_mod(self.s_pow2_mod[msb], term), self.prefix_sum_x_mod(r))
def block_sum_mod(self, p, m):
s_x = self.prefix_sum_x_mod(m + 1)
s1 = sum_linear_mod(m)
s2 = sum_square_mod(m)
sC2 = sum_choose2_mod(m)
two_p_minus_one = (2 * (p % MOD) - 1) % MOD
total = add_mod(add_mod(add_mod(s_x, s2), sC2), mul_mod(two_p_minus_one, s1))
odd_cnt = (m + 1) // 2
odd_sq = mul_mod(odd_cnt % MOD, odd_cnt % MOD)
correction = add_mod(mul_mod(odd_cnt % MOD, two_p_minus_one), mul_mod(2, odd_sq))
return sub_mod(total, correction)
def worker_func(start_idx, end_idx, blocks, solver):
local_sum = 0
for i in range(start_idx, end_idx):
local_sum = add_mod(local_sum, solver.block_sum_mod(blocks[i][0], blocks[i][1]))
return local_sum
def solve_multiprocessing(N):
solver = Solver488()
blocks = []
for k in range(1, 63):
p = (1 << k) - 1
if p + 1 >= N: break
m = min(p, N - 1 - p)
blocks.append((p, m))
if not blocks:
return 0
threads = multiprocessing.cpu_count()
if threads == 0: threads = 1
threads = min(threads, len(blocks))
if len(blocks) < 8: threads = 1
if threads == 1:
ans = 0
for p, m in blocks:
ans = add_mod(ans, solver.block_sum_mod(p, m))
return ans
chunk_size = (len(blocks) + threads - 1) // threads
tasks = []
for i in range(threads):
start = i * chunk_size
end = min(len(blocks), start + chunk_size)
tasks.append((start, end, blocks, solver))
with multiprocessing.Pool(threads) as pool:
results = pool.starmap(worker_func, tasks)
ans = 0
for res in results:
ans = add_mod(ans, res)
return ans
def solve():
N = 1000000000000000000
ans = solve_multiprocessing(N)
return f"{ans:09d}"
if __name__ == '__main__':
print(solve())
Java
import java.util.ArrayList;
import java.util.List;
import java.util.concurrent.Callable;
import java.util.concurrent.ExecutionException;
import java.util.concurrent.ExecutorService;
import java.util.concurrent.Executors;
import java.util.concurrent.Future;
public class Euler488 {
private static final long MOD = 1_000_000_000L;
private static long addMod(long a, long b) {
a += b;
if (a >= MOD)
a -= MOD;
return a;
}
private static long subMod(long a, long b) {
return (a >= b) ? (a - b) : (a + MOD - b);
}
private static long mulMod(long a, long b) {
return ((a % MOD) * (b % MOD)) % MOD;
}
private static long sumLinearMod(long m) {
if (m == 0)
return 0;
long a = m, b = m + 1;
if ((a & 1) == 0)
a >>= 1;
else
b >>= 1;
return mulMod(a, b);
}
private static long sumSquareMod(long m) {
if (m == 0)
return 0;
long a = m, b = m + 1, c = 2 * m + 1;
if ((a & 1) == 0)
a >>= 1;
else if ((b & 1) == 0)
b >>= 1;
else
c >>= 1;
if (a % 3 == 0)
a /= 3;
else if (b % 3 == 0)
b /= 3;
else
c /= 3;
return mulMod(mulMod(a, b), c);
}
private static long sumChoose2Mod(long m) {
if (m <= 1)
return 0;
long a = m, b = m + 1, c = m - 1;
if ((a & 1) == 0)
a >>= 1;
else if ((b & 1) == 0)
b >>= 1;
else
c >>= 1;
if (a % 3 == 0)
a /= 3;
else if (b % 3 == 0)
b /= 3;
else
c /= 3;
return mulMod(mulMod(a, b), c);
}
static class Block {
long p, m;
Block(long p, long m) {
this.p = p;
this.m = m;
}
}
static class Solver {
long[] sPow2Mod = new long[63];
long[] cMod = new long[63];
Solver() {
precomputePrefixXorSums();
}
private void precomputePrefixXorSums() {
for (int m = 0; m <= 62; ++m) {
long p = 1L << m;
long x = p;
long y = 3L * p - 1L;
if ((x & 1L) == 0)
x >>= 1;
else
y >>= 1;
cMod[m] = mulMod(x, y);
}
for (int m = 0; m < 62; ++m) {
long p = (1L << m) % MOD;
long term = mulMod(p, cMod[m]);
sPow2Mod[m + 1] = addMod(addMod(sPow2Mod[m], sPow2Mod[m]), term);
}
}
long prefixSumXMod(long n) {
if (n <= 1)
return 0;
int msb = 63 - Long.numberOfLeadingZeros(n);
long p = 1L << msb;
if (n == p)
return sPow2Mod[msb];
long r = n - p;
long term = mulMod(r % MOD, cMod[msb]);
return addMod(addMod(sPow2Mod[msb], term), prefixSumXMod(r));
}
long blockSumMod(long p, long m) {
long sX = prefixSumXMod(m + 1);
long s1 = sumLinearMod(m);
long s2 = sumSquareMod(m);
long sC2 = sumChoose2Mod(m);
long twoPMinusOne = (2L * (p % MOD) - 1L) % MOD;
if (twoPMinusOne < 0)
twoPMinusOne += MOD;
long total = 0;
total = addMod(total, sX);
total = addMod(total, s2);
total = addMod(total, sC2);
total = addMod(total, mulMod(twoPMinusOne, s1));
long oddCnt = (m + 1L) >> 1;
long oddSq = mulMod(oddCnt % MOD, oddCnt % MOD);
long correction = addMod(mulMod(oddCnt % MOD, twoPMinusOne), mulMod(2L, oddSq));
return subMod(total, correction);
}
long solveMod(long N, boolean allowMultithreading) throws InterruptedException, ExecutionException {
List<Block> blocks = new ArrayList<>();
for (int k = 1; k <= 62; ++k) {
long p = (1L << k) - 1L;
if (p + 1L >= N)
break;
long m = Math.min(p, N - 1L - p);
blocks.add(new Block(p, m));
}
if (blocks.isEmpty())
return 0;
int workers = 1;
if (allowMultithreading) {
int hw = Runtime.getRuntime().availableProcessors();
if (hw > 0)
workers = Math.min(hw, blocks.size());
if (blocks.size() < 8)
workers = 1;
}
if (workers == 1) {
long ans = 0;
for (Block blk : blocks) {
ans = addMod(ans, blockSumMod(blk.p, blk.m));
}
return ans;
}
ExecutorService executor = Executors.newFixedThreadPool(workers);
List<Future<Long>> futures = new ArrayList<>();
for (int tid = 0; tid < workers; ++tid) {
final int currentTid = tid;
final int finalWorkers = workers;
futures.add(executor.submit(() -> {
int L = blocks.size() * currentTid / finalWorkers;
int R = blocks.size() * (currentTid + 1) / finalWorkers;
long local = 0;
for (int i = L; i < R; ++i) {
local = addMod(local, blockSumMod(blocks.get(i).p, blocks.get(i).m));
}
return local;
}));
}
long ans = 0;
for (Future<Long> f : futures) {
ans = addMod(ans, f.get());
}
executor.shutdown();
return ans;
}
}
public static void main(String[] args) throws Exception {
Solver solver = new Solver();
long ans = solver.solveMod(1_000_000_000_000_000_000L, true);
System.out.printf("%09d\n", ans);
}
}