Problem 986: Another Infinite Game
View on Project EulerProject Euler Problem 986 Solution
EulerSolve provides an optimized solution for Project Euler Problem 986, Another Infinite Game, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary For each pair \(1 \le c,d \le 160\), let \(s=c+d\). We place nonnegative integers on a cycle of length \(s\), start from a single pulse \(k\) at one position, and then repeatedly perform a cyclic in-place averaging update with a floor. The key quantity is \(K(c,d)\), the largest starting value that still eventually dies out to the all-zero state. Once \(K(c,d)\) is known, the problem asks for \[ G(c,d)=2K(c,d)+1, \qquad \sum_{c=1}^{160}\sum_{d=1}^{160} G(c,d). \] A brute-force search over all starting values is impossible because the thresholds become large. The implementations succeed by exploiting the monotone structure of the update rule, a reduction by \(\gcd(c,d)\), and a further collapse of the reduced problem to one-dimensional threshold tables. Mathematical Approach The whole problem is governed by a discrete dynamical system on a cycle. The useful mathematics is not an abstract template; it is the exact recurrence, the exact invariant regions, and the exact threshold structure used by the search. The cyclic averaging recurrence Index the cycle by \(0,1,\dots,s-1\), and let the state after some number of elementary updates be \(x=(x_0,\dots,x_{s-1})\). The initial state is \[ x_{s-1}=k, \qquad x_i=0 \text{ for } i\ne s-1. \] At each elementary step, only one coordinate changes....
Detailed mathematical approach
Problem Summary
For each pair \(1 \le c,d \le 160\), let \(s=c+d\). We place nonnegative integers on a cycle of length \(s\), start from a single pulse \(k\) at one position, and then repeatedly perform a cyclic in-place averaging update with a floor. The key quantity is \(K(c,d)\), the largest starting value that still eventually dies out to the all-zero state.
Once \(K(c,d)\) is known, the problem asks for
\[ G(c,d)=2K(c,d)+1, \qquad \sum_{c=1}^{160}\sum_{d=1}^{160} G(c,d). \]
A brute-force search over all starting values is impossible because the thresholds become large. The implementations succeed by exploiting the monotone structure of the update rule, a reduction by \(\gcd(c,d)\), and a further collapse of the reduced problem to one-dimensional threshold tables.
Mathematical Approach
The whole problem is governed by a discrete dynamical system on a cycle. The useful mathematics is not an abstract template; it is the exact recurrence, the exact invariant regions, and the exact threshold structure used by the search.
The cyclic averaging recurrence
Index the cycle by \(0,1,\dots,s-1\), and let the state after some number of elementary updates be \(x=(x_0,\dots,x_{s-1})\). The initial state is
\[ x_{s-1}=k, \qquad x_i=0 \text{ for } i\ne s-1. \]
At each elementary step, only one coordinate changes. If the current head position is \(h\), then the update is
\[ x_h \leftarrow \left\lfloor \frac{x_h + x_{h-d \bmod s}}{2} \right\rfloor, \]
and then the head advances to \(h+1 \bmod s\). Because \(s=c+d\), the same rule can also be written as
\[ x_h \leftarrow \left\lfloor \frac{x_h + x_{h+c \bmod s}}{2} \right\rfloor. \]
This is the exact recurrence implemented in all three languages. The process is asynchronous: values written earlier in the sweep are immediately visible to later updates in the same sweep.
Order preservation and absorbing regions
Two facts drive the search.
First, the local map \((a,b)\mapsto \lfloor(a+b)/2\rfloor\) is monotone in both inputs. So if two initial pulses satisfy \(k_1 \le k_2\), and we run the same update schedule on both, then the componentwise order is preserved forever. Larger initial data can never produce a smaller trajectory than smaller initial data.
Second, the dynamics has two obvious forward-invariant regions:
\[ x\equiv 0 \implies \text{all future states are } 0, \]
and
\[ x_i \ge 1 \text{ for every } i \implies \left\lfloor \frac{x_i+x_j}{2} \right\rfloor \ge 1, \]
so a state with no zeros can never reach the all-zero state later. This is why the implementations can stop as soon as either every entry is zero or every entry is positive.
There is also a simple boundedness invariant:
\[ 0 \le x_i \le k \qquad \text{for every coordinate and every time.} \]
No update can create a negative value or exceed the initial pulse.
The extinction threshold \(K(c,d)\)
Define the extinction predicate
\[ \mathcal{E}_{c,d}(k)= \begin{cases} 1, & \text{if the trajectory started from } k \text{ reaches } 0,\\ 0, & \text{otherwise.} \end{cases} \]
By monotonicity, \(\mathcal{E}_{c,d}(k)\) is non-increasing in \(k\). Therefore the set of extinguishing initial values is an initial interval
\[ \{k\ge 0 : \mathcal{E}_{c,d}(k)=1\}=\{0,1,\dots,K(c,d)\}. \]
This is the mathematical reason binary search is valid: once some \(k\) fails to die out, every larger \(k\) also fails.
Reduction by the greatest common divisor
Let
\[ g=\gcd(c,d), \qquad r_c=\frac{c}{g}, \qquad r_d=\frac{d}{g}. \]
Since the update for position \(h\) reads from \(h+c \bmod s\), every dependency preserves the residue class modulo \(g\). So the cycle of length \(s=c+d\) splits into \(g\) independent smaller cycles, each of length
\[ \frac{s}{g}=r_c+r_d. \]
Only one of those residue classes contains the initial pulse. The other \(g-1\) classes start at zero and stay zero forever. Therefore the extinction threshold depends only on the reduced pair:
\[ K(c,d)=K(r_c,r_d). \]
This reduction is exact, not heuristic. It explains why the final summation always begins by dividing \(c\) and \(d\) by their gcd.
The one-parameter collapse used by the solver
After the gcd reduction, the implementations do not build a full two-dimensional table of reduced thresholds. Instead they use a stronger structural relation and reconstruct every reduced pair from one-dimensional families.
Write
\[ K_1(n)=K(1,n), \qquad K_{d=1}(m)=K(m,1). \]
The reduced threshold used by the solver is then
\[ K(r_c,r_d)= \begin{cases} K_{d=1}(r_c), & r_d=1,\\ K_1\!\left(r_d+\left\lfloor\frac{r_c-1}{2}\right\rfloor\right), & r_d>1. \end{cases} \]
So the two-parameter family is collapsed to the line \(K(1,n)\), plus a small exceptional strip for \(K(m,1)\). In particular, once \(K_1(1),K_1(2),\dots,K_1(239)\) are known, almost every reduced pair needed by the problem is already determined.
Worked example: the threshold for \((c,d)=(2,1)\)
Here \(s=3\), so the initial state is \((0,0,k)\). The update repeatedly averages a position with the previous one on the 3-cycle.
For \(k=3\), one possible compressed trace is
\[ (0,0,3)\to(1,0,3)\to(1,0,1)\to(1,0,0)\to(0,0,0). \]
So \(k=3\) does become extinct.
For \(k=4\), the process behaves differently:
\[ (0,0,4)\to(2,0,4)\to(2,1,4)\to(2,1,2)\to(2,1,1)\to(1,1,1). \]
Once the state is \((1,1,1)\), every coordinate is positive, so extinction is impossible from then on. Hence \(K(2,1)=3\), and therefore
\[ G(2,1)=2\cdot 3+1=7. \]
This small example captures the general pattern: monotonicity turns the dynamical question into a threshold search.
How the Code Works
Exact extinction tests
The C++, Python, and Java implementations simulate the in-place recurrence exactly for a fixed triple \((c,d,k)\). They maintain the cyclic state, the current head, and a counter of how many entries are zero. That counter gives two exact stopping conditions: all zero means extinction has been reached, while no zero means the state has entered the forward-invariant positive region and cannot die out anymore.
Precomputing reduced thresholds
Each exact threshold is found by doubling an upper bound until extinction fails, then binary-searching between the last extinguishing value and the first non-extinguishing one. For the special family \(K_1(n)=K(1,n)\), the implementations use an eight-way cubic predictor indexed by \(n \bmod 8\) to place the initial search window near the true answer. That predictor changes only the speed of the search, not its correctness.
The one-dimensional table \(K_1(1),\dots,K_1(239)\) is the main precomputation. The C++ and Java implementations parallelize it across threads; the Python implementation uses multiple processes when available and falls back to serial evaluation if needed. For the branch \(K(m,1)\), a few small values are computed directly, and the larger ones are reconstructed from the same collapsed family with direct spot checks for safety.
Final reconstruction of the sum
For every pair \((c,d)\), the implementation first reduces it to \((r_c,r_d)\) by dividing by \(\gcd(c,d)\). It then recovers the needed threshold from the precomputed one-dimensional tables, forms
\[ G(c,d)=2K(r_c,r_d)+1, \]
and accumulates the result over all \(160^2\) pairs. Because the final answer is far larger than a machine word, the total is stored with arbitrary-precision integer arithmetic.
Complexity Analysis
Let \(\tau(c,d,k)\) be the number of elementary updates performed before the extinction test for \((c,d,k)\) reaches one of its two exact stopping conditions. Then one threshold query costs
\[ O\!\big(\tau(c,d,\cdot)\log K(c,d)\big), \]
because doubling plus binary search makes only logarithmically many extinction tests.
The dominant work is the precomputation of the table \(K_1(1),\dots,K_1(239)\). After that, the final double sum over all \(c,d\le 160\) is cheap: each pair needs one gcd, one table lookup, and a few arithmetic operations. Memory usage is small. A single simulation stores at most \(c+d\le 320\) integers, and the precomputed tables themselves have only a few hundred entries.
Footnotes and References
- Problem page: Project Euler 986
- Greatest common divisor: Wikipedia - Greatest common divisor
- Binary search algorithm: Wikipedia - Binary search algorithm
- Discrete dynamical system: Wikipedia - Discrete dynamical system
- Asynchronous iteration: Wikipedia - Asynchronous iteration
- Floor and ceiling functions: Wikipedia - Floor and ceiling functions
Problem 986 source code
C++
#include <algorithm>
#include <array>
#include <atomic>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <vector>
#include <pthread.h>
#include <unistd.h>
#include <boost/multiprecision/cpp_int.hpp>
namespace {
using u64 = std::uint64_t;
using bigint = boost::multiprecision::cpp_int;
bool extinct_for_k(int c, int d, u64 k) {
if (k == 0) return true;
const int s = c + d;
std::array<u64, 321> window{};
window[static_cast<std::size_t>(s - 1)] = k;
int zero_count = s - 1;
int head = 0;
int tap = s - d;
u64 steps = 0;
constexpr u64 CHECK_MASK = 63ULL;
for (;;) {
const u64 old = window[static_cast<std::size_t>(head)];
const u64 next = (window[static_cast<std::size_t>(head)] + window[static_cast<std::size_t>(tap)]) >> 1U;
window[static_cast<std::size_t>(head)] = next;
zero_count += static_cast<int>(old != 0 && next == 0) - static_cast<int>(old == 0 && next != 0);
if (++head == s) head = 0;
if (++tap == s) tap = 0;
++steps;
if ((steps & CHECK_MASK) == 0) {
if (zero_count == s) return true;
if (zero_count == 0) return false;
}
}
}
bool extinct_for_k_c1(int n, u64 k) {
if (k == 0) return true;
const int s = n + 1;
std::array<u64, 240> window{};
window[static_cast<std::size_t>(s - 1)] = k;
int zero_count = s - 1;
for (;;) {
for (int i = 0; i < s; ++i) {
const int j = (i + 1 == s) ? 0 : (i + 1);
const u64 old = window[static_cast<std::size_t>(i)];
const u64 next = (window[static_cast<std::size_t>(i)] + window[static_cast<std::size_t>(j)]) >> 1U;
window[static_cast<std::size_t>(i)] = next;
zero_count += static_cast<int>(old != 0 && next == 0) - static_cast<int>(old == 0 && next != 0);
}
if (zero_count == s) return true;
if (zero_count == 0) return false;
}
}
u64 compute_k(int c, int d) {
u64 lo = 0;
u64 hi = 1;
while (extinct_for_k(c, d, hi)) {
lo = hi;
hi <<= 1U;
}
while (lo + 1 < hi) {
const u64 mid = lo + ((hi - lo) >> 1U);
if (extinct_for_k(c, d, mid)) {
lo = mid;
} else {
hi = mid;
}
}
return lo;
}
u64 estimate_k1(int n) {
static constexpr std::array<std::array<long double, 4>, 8> coef{{
{{0.17656501777829126L, 0.35981185819773187L, -9.093343254260652L, 128.0591554042931L}},
{{0.17791093179771922L, -0.02526258663290543L, 17.778074509864577L, -228.639325004742L}},
{{0.17688473006194433L, 0.30195121097992167L, -8.649989398001244L, 159.32792700127638L}},
{{0.1774911030303856L, 0.09144288721929997L, 9.76665768770627L, -161.46727161611716L}},
{{0.17665123701842447L, 0.37639947984877686L, -14.382073443985455L, 264.68150347640244L}},
{{0.17652781023196232L, 0.37645445228054497L, -10.389054848768398L, 121.34965708864206L}},
{{0.17725060801487988L, 0.15363937288855242L, 7.49278183636027L, -167.17442628344003L}},
{{0.17658285799047665L, 0.3729423126656912L, -11.25711787747357L, 179.7243004469314L}},
}};
const auto& c = coef[static_cast<std::size_t>(n & 7)];
const long double x = static_cast<long double>(n);
const long double est = ((c[0] * x + c[1]) * x + c[2]) * x + c[3];
return est > 0.0L ? static_cast<u64>(std::llround(est)) : 0ULL;
}
u64 compute_k1_from_estimate(int n) {
const u64 guess = estimate_k1(n);
u64 lo = (guess > 2048ULL) ? (guess - 2048ULL) : 0ULL;
u64 hi = guess + 2048ULL;
while (lo > 0 && !extinct_for_k_c1(n, lo)) {
hi = lo;
lo >>= 1U;
}
while (extinct_for_k_c1(n, hi)) {
lo = hi;
hi <<= 1U;
}
while (lo + 1 < hi) {
const u64 mid = lo + ((hi - lo) >> 1U);
if (extinct_for_k_c1(n, mid)) {
lo = mid;
} else {
hi = mid;
}
}
return lo;
}
struct K1WorkerCtx {
std::atomic<int>* next_n = nullptr;
std::array<u64, 240>* k1 = nullptr;
};
void* k1_worker_main(void* raw_ctx) {
auto* ctx = static_cast<K1WorkerCtx*>(raw_ctx);
for (;;) {
const int n = ctx->next_n->fetch_add(1, std::memory_order_relaxed);
if (n > 239) break;
(*ctx->k1)[static_cast<std::size_t>(n)] = compute_k1_from_estimate(n);
}
return nullptr;
}
} // namespace
int main() {
std::ios::sync_with_stdio(false);
std::cin.tie(nullptr);
std::array<u64, 240> k1{};
std::array<u64, 161> k_d1{};
std::atomic<int> next_n(1);
K1WorkerCtx worker_ctx{&next_n, &k1};
long cpu_count = ::sysconf(_SC_NPROCESSORS_ONLN);
if (cpu_count < 1) cpu_count = 1;
const int thread_count = static_cast<int>(std::min<long>(64, cpu_count * 2));
std::vector<pthread_t> threads(static_cast<std::size_t>(thread_count));
for (int t = 0; t < thread_count; ++t) {
const int rc = pthread_create(&threads[static_cast<std::size_t>(t)], nullptr, k1_worker_main, &worker_ctx);
assert(rc == 0);
}
for (int t = 0; t < thread_count; ++t) {
const int rc = pthread_join(threads[static_cast<std::size_t>(t)], nullptr);
assert(rc == 0);
}
for (int c = 1; c <= 10; ++c) {
k_d1[static_cast<std::size_t>(c)] = compute_k(c, 1);
}
for (int c = 11; c <= 160; ++c) {
k_d1[static_cast<std::size_t>(c)] = k1[static_cast<std::size_t>((c + 1) / 2)];
}
auto k_for_reduced = [&](int rc, int rd) -> u64 {
if (rd == 1) return k_d1[static_cast<std::size_t>(rc)];
return k1[static_cast<std::size_t>(rd + (rc - 1) / 2)];
};
auto G = [&](int c, int d) -> u64 {
const int g = std::gcd(c, d);
const int rc = c / g;
const int rd = d / g;
return 2 * k_for_reduced(rc, rd) + 1;
};
assert(G(2, 1) == 7);
assert(G(1, 2) == 7);
assert(G(3, 1) == 11);
assert(G(2, 2) == 3);
assert(G(1, 3) == 15);
assert(k_d1[11] == compute_k(11, 1));
assert(k_d1[57] == compute_k(57, 1));
assert(k_d1[160] == compute_k(160, 1));
bigint sum = 0;
for (int c = 1; c <= 160; ++c) {
for (int d = 1; d <= 160; ++d) {
sum += G(c, d);
}
}
std::cout << sum << '\n';
return 0;
}
Python
import math
import os
from concurrent.futures import ProcessPoolExecutor
def extinct_for_k(c, d, k):
if k == 0:
return True
s = c + d
window = [0] * s
window[s - 1] = k
zero_count = s - 1
head = 0
tap = s - d
steps = 0
while True:
old = window[head]
nxt = (old + window[tap]) >> 1
window[head] = nxt
if old != 0:
if nxt == 0:
zero_count += 1
elif nxt != 0:
zero_count -= 1
head += 1
if head == s:
head = 0
tap += 1
if tap == s:
tap = 0
steps += 1
if (steps & 63) == 0:
if zero_count == s:
return True
if zero_count == 0:
return False
def extinct_for_k_c1(n, k):
if k == 0:
return True
s = n + 1
window = [0] * s
window[s - 1] = k
zero_count = s - 1
while True:
for i in range(s):
j = i + 1
if j == s:
j = 0
old = window[i]
nxt = (old + window[j]) >> 1
window[i] = nxt
if old != 0:
if nxt == 0:
zero_count += 1
elif nxt != 0:
zero_count -= 1
if zero_count == s:
return True
if zero_count == 0:
return False
def compute_k(c, d):
lo = 0
hi = 1
while extinct_for_k(c, d, hi):
lo = hi
hi <<= 1
while lo + 1 < hi:
mid = lo + ((hi - lo) >> 1)
if extinct_for_k(c, d, mid):
lo = mid
else:
hi = mid
return lo
COEF = [
(0.17656501777829126, 0.35981185819773187, -9.093343254260652, 128.0591554042931),
(0.17791093179771922, -0.02526258663290543, 17.778074509864577, -228.639325004742),
(0.17688473006194433, 0.30195121097992167, -8.649989398001244, 159.32792700127638),
(0.1774911030303856, 0.09144288721929997, 9.76665768770627, -161.46727161611716),
(0.17665123701842447, 0.37639947984877686, -14.382073443985455, 264.68150347640244),
(0.17652781023196232, 0.37645445228054497, -10.389054848768398, 121.34965708864206),
(0.17725060801487988, 0.15363937288855242, 7.49278183636027, -167.17442628344003),
(0.17658285799047665, 0.3729423126656912, -11.25711787747357, 179.7243004469314),
]
def estimate_k1(n):
c0, c1, c2, c3 = COEF[n & 7]
x = float(n)
est = ((c0 * x + c1) * x + c2) * x + c3
if est <= 0.0:
return 0
return int(round(est))
def compute_k1_from_estimate(n):
guess = estimate_k1(n)
lo = guess - 2048
if lo < 0:
lo = 0
hi = guess + 2048
while lo > 0 and not extinct_for_k_c1(n, lo):
hi = lo
lo >>= 1
while extinct_for_k_c1(n, hi):
lo = hi
hi <<= 1
while lo + 1 < hi:
mid = lo + ((hi - lo) >> 1)
if extinct_for_k_c1(n, mid):
lo = mid
else:
hi = mid
return lo
def compute_all_k1():
k1 = [0] * 240
n_values = list(range(1, 240))
worker_count = min(64, max(1, (os.cpu_count() or 1) * 2))
if worker_count == 1:
for n in n_values:
k1[n] = compute_k1_from_estimate(n)
return k1
try:
with ProcessPoolExecutor(max_workers=worker_count) as executor:
for n, value in zip(n_values, executor.map(compute_k1_from_estimate, n_values, chunksize=1)):
k1[n] = value
except Exception:
for n in n_values:
k1[n] = compute_k1_from_estimate(n)
return k1
def solve():
k1 = compute_all_k1()
k_d1 = [0] * 161
for c in range(1, 11):
k_d1[c] = compute_k(c, 1)
for c in range(11, 161):
k_d1[c] = k1[(c + 1) // 2]
def k_for_reduced(rc, rd):
if rd == 1:
return k_d1[rc]
return k1[rd + (rc - 1) // 2]
def g_value(c, d):
g = math.gcd(c, d)
rc = c // g
rd = d // g
return 2 * k_for_reduced(rc, rd) + 1
assert g_value(2, 1) == 7
assert g_value(1, 2) == 7
assert g_value(3, 1) == 11
assert g_value(2, 2) == 3
assert g_value(1, 3) == 15
assert k_d1[11] == compute_k(11, 1)
assert k_d1[57] == compute_k(57, 1)
assert k_d1[160] == compute_k(160, 1)
total = 0
for c in range(1, 161):
for d in range(1, 161):
total += g_value(c, d)
return str(total)
if __name__ == "__main__":
print(solve())
Java
import java.math.BigInteger;
import java.util.concurrent.atomic.AtomicInteger;
public class Euler986 {
private static final double[][] COEF = {
{0.17656501777829126, 0.35981185819773187, -9.093343254260652, 128.0591554042931},
{0.17791093179771922, -0.02526258663290543, 17.778074509864577, -228.639325004742},
{0.17688473006194433, 0.30195121097992167, -8.649989398001244, 159.32792700127638},
{0.1774911030303856, 0.09144288721929997, 9.76665768770627, -161.46727161611716},
{0.17665123701842447, 0.37639947984877686, -14.382073443985455, 264.68150347640244},
{0.17652781023196232, 0.37645445228054497, -10.389054848768398, 121.34965708864206},
{0.17725060801487988, 0.15363937288855242, 7.49278183636027, -167.17442628344003},
{0.17658285799047665, 0.3729423126656912, -11.25711787747357, 179.7243004469314},
};
private static boolean extinctForK(int c, int d, long k) {
if (k == 0) {
return true;
}
int s = c + d;
long[] window = new long[s];
window[s - 1] = k;
int zeroCount = s - 1;
int head = 0;
int tap = s - d;
long steps = 0;
while (true) {
long old = window[head];
long next = (window[head] + window[tap]) >>> 1;
window[head] = next;
if (old != 0) {
if (next == 0) {
zeroCount++;
}
} else if (next != 0) {
zeroCount--;
}
head++;
if (head == s) {
head = 0;
}
tap++;
if (tap == s) {
tap = 0;
}
steps++;
if ((steps & 63L) == 0L) {
if (zeroCount == s) {
return true;
}
if (zeroCount == 0) {
return false;
}
}
}
}
private static boolean extinctForKC1(int n, long k) {
if (k == 0) {
return true;
}
int s = n + 1;
long[] window = new long[s];
window[s - 1] = k;
int zeroCount = s - 1;
while (true) {
for (int i = 0; i < s; i++) {
int j = i + 1;
if (j == s) {
j = 0;
}
long old = window[i];
long next = (window[i] + window[j]) >>> 1;
window[i] = next;
if (old != 0) {
if (next == 0) {
zeroCount++;
}
} else if (next != 0) {
zeroCount--;
}
}
if (zeroCount == s) {
return true;
}
if (zeroCount == 0) {
return false;
}
}
}
private static long computeK(int c, int d) {
long lo = 0;
long hi = 1;
while (extinctForK(c, d, hi)) {
lo = hi;
hi <<= 1;
}
while (lo + 1 < hi) {
long mid = lo + ((hi - lo) >>> 1);
if (extinctForK(c, d, mid)) {
lo = mid;
} else {
hi = mid;
}
}
return lo;
}
private static long estimateK1(int n) {
double[] c = COEF[n & 7];
double x = n;
double est = ((c[0] * x + c[1]) * x + c[2]) * x + c[3];
if (est <= 0.0) {
return 0;
}
return Math.round(est);
}
private static long computeK1FromEstimate(int n) {
long guess = estimateK1(n);
long lo = guess > 2048L ? guess - 2048L : 0L;
long hi = guess + 2048L;
while (lo > 0 && !extinctForKC1(n, lo)) {
hi = lo;
lo >>>= 1;
}
while (extinctForKC1(n, hi)) {
lo = hi;
hi <<= 1;
}
while (lo + 1 < hi) {
long mid = lo + ((hi - lo) >>> 1);
if (extinctForKC1(n, mid)) {
lo = mid;
} else {
hi = mid;
}
}
return lo;
}
private static long gcd(long a, long b) {
while (b != 0) {
long t = a % b;
a = b;
b = t;
}
return a;
}
private static long gValue(int c, int d, long[] kD1, long[] k1) {
long g = gcd(c, d);
int rc = (int) (c / g);
int rd = (int) (d / g);
long k;
if (rd == 1) {
k = kD1[rc];
} else {
k = k1[rd + (rc - 1) / 2];
}
return 2L * k + 1L;
}
private static void require(boolean condition, String message) {
if (!condition) {
throw new IllegalStateException(message);
}
}
public static String solve() {
long[] k1 = new long[240];
long[] kD1 = new long[161];
AtomicInteger nextN = new AtomicInteger(1);
int cpuCount = Runtime.getRuntime().availableProcessors();
if (cpuCount < 1) {
cpuCount = 1;
}
int threadCount = Math.min(64, cpuCount * 2);
Thread[] threads = new Thread[threadCount];
for (int t = 0; t < threadCount; t++) {
threads[t] = new Thread(() -> {
while (true) {
int n = nextN.getAndIncrement();
if (n > 239) {
break;
}
k1[n] = computeK1FromEstimate(n);
}
});
threads[t].start();
}
for (Thread thread : threads) {
try {
thread.join();
} catch (InterruptedException e) {
Thread.currentThread().interrupt();
throw new RuntimeException("Thread interrupted", e);
}
}
for (int c = 1; c <= 10; c++) {
kD1[c] = computeK(c, 1);
}
for (int c = 11; c <= 160; c++) {
kD1[c] = k1[(c + 1) / 2];
}
require(gValue(2, 1, kD1, k1) == 7, "Check G(2,1)");
require(gValue(1, 2, kD1, k1) == 7, "Check G(1,2)");
require(gValue(3, 1, kD1, k1) == 11, "Check G(3,1)");
require(gValue(2, 2, kD1, k1) == 3, "Check G(2,2)");
require(gValue(1, 3, kD1, k1) == 15, "Check G(1,3)");
require(kD1[11] == computeK(11, 1), "Check k_d1[11]");
require(kD1[57] == computeK(57, 1), "Check k_d1[57]");
require(kD1[160] == computeK(160, 1), "Check k_d1[160]");
BigInteger sum = BigInteger.ZERO;
for (int c = 1; c <= 160; c++) {
for (int d = 1; d <= 160; d++) {
sum = sum.add(BigInteger.valueOf(gValue(c, d, kD1, k1)));
}
}
return sum.toString();
}
public static void main(String[] args) {
System.out.println(solve());
}
}