Problem 433: Steps in Euclid's Algorithm
View on Project EulerProject Euler Problem 433 Solution
EulerSolve provides an optimized solution for Project Euler Problem 433, Steps in Euclid's Algorithm, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary For positive integers \(x\) and \(y\), let \(E(x,y)\) denote the number of divisions performed by Euclid's algorithm until the remainder becomes \(0\). The problem asks for $$S(N)=\sum_{1\le x,y\le N} E(x,y)$$ at the very large value \(N=5\cdot 10^6\). A direct double loop is hopeless: there are \(N^2\) ordered pairs, and each pair would still require a separate Euclidean run. The solution therefore reorganizes the sum so that most of the work is done by arithmetic formulas, Möbius inversion, and floor-sum recurrences. Mathematical Approach 1. Reduce the Full Sum to the Lower Triangle When \(x>y\), starting from \((y,x)\) performs one swap before entering exactly the same remainder chain, so $$E(y,x)=E(x,y)+1 \qquad (x>y).$$ Also, \(E(x,x)=1\) for every diagonal pair. If we define $$L(N)=\sum_{1\le y<x\le N} E(x,y),$$ then the whole ordered sum becomes $$\begin{aligned} S(N) &=\sum_{x=1}^{N}E(x,x)+\sum_{1\le y<x\le N}\bigl(E(x,y)+E(y,x)\bigr) \\ &=N+\sum_{1\le y<x\le N}\bigl(2E(x,y)+1\bigr) \\ &=2L(N)+\frac{N(N+1)}{2}. \end{aligned}$$ So it is enough to evaluate the lower-triangular sum \(L(N)\). 2. Peel Off the Universal First Steps Every pair with \(x>y\) needs at least one Euclidean step....
Detailed mathematical approach
Problem Summary
For positive integers \(x\) and \(y\), let \(E(x,y)\) denote the number of divisions performed by Euclid's algorithm until the remainder becomes \(0\). The problem asks for
$$S(N)=\sum_{1\le x,y\le N} E(x,y)$$
at the very large value \(N=5\cdot 10^6\). A direct double loop is hopeless: there are \(N^2\) ordered pairs, and each pair would still require a separate Euclidean run. The solution therefore reorganizes the sum so that most of the work is done by arithmetic formulas, Möbius inversion, and floor-sum recurrences.
Mathematical Approach
1. Reduce the Full Sum to the Lower Triangle
When \(x>y\), starting from \((y,x)\) performs one swap before entering exactly the same remainder chain, so
$$E(y,x)=E(x,y)+1 \qquad (x>y).$$
Also, \(E(x,x)=1\) for every diagonal pair. If we define
$$L(N)=\sum_{1\le y<x\le N} E(x,y),$$
then the whole ordered sum becomes
$$\begin{aligned} S(N) &=\sum_{x=1}^{N}E(x,x)+\sum_{1\le y<x\le N}\bigl(E(x,y)+E(y,x)\bigr) \\ &=N+\sum_{1\le y<x\le N}\bigl(2E(x,y)+1\bigr) \\ &=2L(N)+\frac{N(N+1)}{2}. \end{aligned}$$
So it is enough to evaluate the lower-triangular sum \(L(N)\).
2. Peel Off the Universal First Steps
Every pair with \(x>y\) needs at least one Euclidean step. There is also an easily recognized second step: if \(x<2y\), then the first quotient is \(1\), so after one subtraction we are left with \((y,x-y)\), which is still a non-terminal pair. Hence
$$E(x,y)=1+\mathbf{1}_{x<2y}+R(x,y),\qquad R(x,y)\ge 0.$$
The easy contribution is therefore
$$B(N)=\sum_{1\le y<x\le N}\left(1+\mathbf{1}_{x<2y}\right).$$
The first part is just the number of pairs below the diagonal:
$$\sum_{1\le y<x\le N}1=\frac{N(N-1)}{2}.$$
For the indicator term, count the pairs with \(y<x<2y\). Writing \(N=2m\) or \(N=2m+1\), one gets
$$\sum_{1\le y<x\le 2m}\mathbf{1}_{x<2y}=m(m-1),$$
$$\sum_{1\le y<x\le 2m+1}\mathbf{1}_{x<2y}=m^2.$$
So the closed base term is
$$B(N)= \begin{cases} (3m-2)m, & N=2m, \\ 3m^2+m, & N=2m+1. \end{cases}$$
This is exactly the closed form used by the implementations before any expensive computation begins.
3. Residual Contribution and the GCD Structure
The remaining term \(R(x,y)\) measures everything beyond those universal first one or two steps. Because Euclid's algorithm is unchanged by scaling both arguments by the same factor,
$$R(dx,dy)=R(x,y)\qquad (d\ge 1).$$
That makes the gcd decomposition natural: the residual depends only on the primitive core of the pair.
Another useful observation is that \(R(x,y)=0\) unless the larger coordinate is at least \(5\). Indeed, the first primitive pairs with more than the peeled-off base behavior are \((5,2)\) and \((5,3)\). This is why the Möbius layer only needs to run up to \(\lfloor N/5\rfloor\).
4. Möbius Inversion
Let \(\mu\) be the Möbius function and let
$$M(t)=\sum_{d\le t}\mu(d)$$
be its prefix sum. The fast method isolates primitive pairs by Möbius inversion and rewrites the residual part of the lower triangle in the form
$$L(N)=B(N)+2\sum_{d=1}^{\lfloor N/5\rfloor}\mu(d)\,A\!\left(\left\lfloor\frac{N}{d}\right\rfloor\right),$$
where \(A(V)\) is an unrestricted auxiliary counting function. The role of \(\mu(d)\) is to remove contributions coming from non-coprime multiples. This is why all three implementations build Möbius values and their prefix sums before the main accumulation loop.
5. Reverse Euclid and the Auxiliary Function
To compute \(A(V)\), the implementation reverses one step of Euclid's algorithm. If a reduced tail pair is \((u,v)\) with \(u>v\ge 1\), then every predecessor that reaches it after one division has the shape
$$\bigl(ku+v,\;u\bigr),\qquad k\ge 1.$$
Imposing the upper bound \(ku+v\le V\) turns the counting problem into lattice-point counting under linear inequalities. After summing all admissible predecessors of all reduced tails, the required quantities reduce to floor sums of the form
$$\sum_t \left\lfloor\frac{\alpha t+\beta}{\gamma}\right\rfloor,$$
together with a second weighted variant of the same type. The implementation evaluates both of these by a recursive Euclidean-style transform, which is the key reason the method stays fast.
6. Why the Loop Stops at \(\sqrt{V}\)
The predecessor lattice is symmetric between two coordinate regions. Because of that symmetry, the auxiliary computation only needs to enumerate reduced tail pairs with larger coordinate at most \(\lfloor\sqrt{V}\rfloor\). One accumulated sum covers the full asymmetric region, another covers the square overlap, and the final combination is a doubled total minus that overlap correction. This explains the integer-square-root bound seen in the code.
7. Harmonic Interval Grouping
The Möbius sum is not evaluated one divisor at a time. The quotient
$$q=\left\lfloor\frac{N}{d}\right\rfloor$$
stays constant on intervals \(l\le d\le r\), where
$$r=\left\lfloor\frac{N}{q}\right\rfloor.$$
So an entire block contributes
$$\sum_{d=l}^{r}\mu(d)\,A(q)=\bigl(M(r)-M(l-1)\bigr)A(q).$$
Only \(O(\sqrt{N})\) distinct quotients occur, so the expensive auxiliary function is called much fewer times than in a linear scan.
Worked Check: \(N=10\)
For \(N=10\), we have \(m=5\), so the base part is
$$B(10)=(3\cdot 5-2)\cdot 5=65.$$
The residual term computed by the fast Möbius-floor-sum pipeline is \(18\), hence
$$L(10)=65+18=83.$$
Returning to the full ordered sum gives
$$S(10)=2L(10)+\frac{10\cdot 11}{2}=2\cdot 83+55=221,$$
which matches the small checkpoint used to validate the implementation.
How the Code Works
The C++, Python, and Java implementations all follow the same plan. They first build the Möbius table and its prefix sums up to \(\lfloor N/5\rfloor\). They then evaluate the closed base term \(B(N)\), generate the harmonic intervals where \(\lfloor N/d\rfloor\) is constant, and for each distinct quotient compute the reverse-Euclid auxiliary kernel via recursive floor sums. Each block is weighted by the corresponding Möbius prefix difference, added to the base term, and finally converted back to the full ordered total through
$$S(N)=2L(N)+\frac{N(N+1)}{2}.$$
Complexity Analysis
A naive method would inspect all \(N^2\) ordered pairs and run Euclid's algorithm separately, which is entirely infeasible for \(N=5\cdot 10^6\). The optimized method uses a linear-style Möbius sieve up to \(\lfloor N/5\rfloor\), only \(O(\sqrt{N})\) distinct quotient blocks in the outer summation, and a recursive floor-sum kernel whose depth is logarithmic in the relevant parameters. In practice this turns a quadratic enumeration into a heavily compressed arithmetic computation. The main memory cost is the Möbius and prefix tables, so the space usage is \(O(N/5)\).
References
- Problem page: https://projecteuler.net/problem=433
- Euclidean algorithm: Wikipedia — Euclidean algorithm
- Möbius function and inversion: Wikipedia — Möbius inversion formula
- Continued fractions and Euclid's algorithm: Wikipedia — Continued fraction
- Graham, Knuth, Patashnik, Concrete Mathematics, 2nd ed., sections on floor sums and Möbius inversion.
Problem 433 source code
C++
#include <algorithm>
#include <atomic>
#include <chrono>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <limits>
#include <string>
#include <thread>
#include <vector>
namespace {
using i32 = std::int32_t;
using i64 = std::int64_t;
using i8 = std::int8_t;
using u32 = std::uint32_t;
using u64 = std::uint64_t;
using i128 = __int128_t;
using u128 = unsigned __int128;
constexpr u32 kDefaultN = 5'000'000;
struct Options {
u32 n = kDefaultN;
bool allow_multithreading = true;
unsigned requested_threads = 0;
bool run_validation = true;
u32 frontier_split_x = 128;
u32 chunk_divisor = 256;
};
struct FloorPair {
i128 sum1 = 0;
i128 sum2 = 0;
};
struct IntervalTask {
u32 q = 0;
i32 mu_sum = 0;
};
unsigned choose_thread_count(bool allow_multithreading,
unsigned requested_threads,
std::size_t workload) {
if (!allow_multithreading || workload < 50'000) {
return 1;
}
unsigned threads = requested_threads;
if (threads == 0) {
threads = std::thread::hardware_concurrency();
if (threads == 0) {
threads = 1;
}
}
return std::max(1u, std::min<unsigned>(threads, static_cast<unsigned>(workload)));
}
bool parse_unsigned_after_prefix(const std::string& arg,
const char* prefix,
u32& value_out) {
const std::string p(prefix);
if (arg.rfind(p, 0) != 0) {
return false;
}
const std::string tail = arg.substr(p.size());
if (tail.empty()) {
return false;
}
std::uint64_t parsed = 0;
for (const char c : tail) {
if (c < '0' || c > '9') {
return false;
}
parsed = parsed * 10 + static_cast<std::uint64_t>(c - '0');
if (parsed > static_cast<std::uint64_t>(std::numeric_limits<u32>::max())) {
return false;
}
}
value_out = static_cast<u32>(parsed);
return true;
}
bool parse_options(int argc, char** argv, Options& options) {
for (int i = 1; i < argc; ++i) {
const std::string arg(argv[i]);
if (arg == "--single-thread") {
options.allow_multithreading = false;
continue;
}
if (arg == "--no-validate") {
options.run_validation = false;
continue;
}
u32 parsed = 0;
if (parse_unsigned_after_prefix(arg, "--n=", parsed)) {
if (parsed == 0) {
std::cerr << "Invalid --n value: must be >= 1\n";
return false;
}
options.n = parsed;
continue;
}
if (parse_unsigned_after_prefix(arg, "--threads=", parsed)) {
options.requested_threads = static_cast<unsigned>(parsed);
continue;
}
if (parse_unsigned_after_prefix(arg, "--split=", parsed)) {
if (parsed == 0) {
std::cerr << "Invalid --split value: must be >= 1\n";
return false;
}
options.frontier_split_x = parsed;
continue;
}
if (parse_unsigned_after_prefix(arg, "--chunk-div=", parsed)) {
if (parsed == 0) {
std::cerr << "Invalid --chunk-div value: must be >= 1\n";
return false;
}
options.chunk_divisor = parsed;
continue;
}
std::cerr << "Unknown argument: " << arg << '\n';
return false;
}
return true;
}
inline i64 floor_div(i64 a, i64 b) {
return (a >= 0) ? (a / b) : -(((-a) + b - 1) / b);
}
u32 euclid_steps(u32 x, u32 y) {
u32 steps = 0;
while (y != 0) {
const u32 r = x % y;
x = y;
y = r;
++steps;
}
return steps;
}
u64 brute_force_s(u32 n) {
u64 total = 0;
for (u32 x = 1; x <= n; ++x) {
for (u32 y = 1; y <= n; ++y) {
total += euclid_steps(x, y);
}
}
return total;
}
FloorPair floor_sum_pair(i64 p, i64 q, i64 m, i64 s1, i64 s2) {
FloorPair result{};
i64 sign = 1;
while (true) {
i64 t = floor_div(m, q);
m -= t * q;
result.sum1 += static_cast<i128>(sign) * t * s1;
result.sum2 += static_cast<i128>(sign) * t * s2;
t = floor_div(p, q);
p -= t * q;
result.sum1 += static_cast<i128>(sign) * t * s1 * (s1 + 1) / 2;
result.sum2 += static_cast<i128>(sign) * t * s2 * (s2 + 1) / 2;
if (p == 0) {
return result;
}
t = (p * s1 + m) / q;
result.sum1 += static_cast<i128>(sign) * s1 * t;
s1 = t;
t = (p * s2 + m) / q;
result.sum2 += static_cast<i128>(sign) * s2 * t;
s2 = t;
std::swap(p, q);
m = -m - 1;
sign = -sign;
}
}
u32 isqrt_u32(u32 n) {
u32 r = static_cast<u32>(std::sqrt(static_cast<long double>(n)));
while (static_cast<u64>(r + 1) * static_cast<u64>(r + 1) <= n) {
++r;
}
while (static_cast<u64>(r) * static_cast<u64>(r) > n) {
--r;
}
return r;
}
i128 calc_smart(u32 val) {
const u32 max_xy = isqrt_u32(val);
i128 result1 = 0;
i128 result2 = 0;
for (u32 x = 2; x <= max_xy; ++x) {
for (u32 y = 1; y < x; ++y) {
const i64 x64 = x;
const i64 y64 = y;
const i64 v64 = val;
const i64 split = v64 / (x64 + y64);
result1 += static_cast<i128>(split) * (split - 1) / 2;
const i64 s1 = v64 / x64 - split;
const i64 s2 = (split <= max_xy) ? (static_cast<i64>(max_xy) - split) : 0;
const FloorPair block = floor_sum_pair(-x64, y64, v64 - split * x64, s1, s2);
result1 += block.sum1;
if (split <= max_xy) {
result2 += static_cast<i128>(split) * (split - 1) / 2;
result2 += block.sum2;
} else {
result2 += static_cast<i128>(max_xy) * (max_xy - 1) / 2;
}
}
}
return 2 * result1 - result2;
}
void build_mobius(u32 limit, std::vector<i8>& mu, std::vector<i32>& prefix) {
std::vector<u32> primes;
std::vector<u32> lp(static_cast<std::size_t>(limit) + 1, 0);
mu.assign(static_cast<std::size_t>(limit) + 1, 0);
prefix.assign(static_cast<std::size_t>(limit) + 1, 0);
mu[1] = 1;
for (u32 i = 2; i <= limit; ++i) {
if (lp[i] == 0) {
lp[i] = i;
primes.push_back(i);
mu[i] = -1;
}
for (const u32 p : primes) {
const u64 v = static_cast<u64>(i) * p;
if (v > limit) {
break;
}
lp[static_cast<u32>(v)] = p;
if (p == lp[i]) {
mu[static_cast<u32>(v)] = 0;
break;
}
mu[static_cast<u32>(v)] = static_cast<i8>(-mu[i]);
}
}
for (u32 i = 1; i <= limit; ++i) {
prefix[i] = prefix[i - 1] + static_cast<i32>(mu[i]);
}
}
i128 base_term(u32 n) {
const i64 m = n / 2;
if ((n & 1U) == 0U) {
return static_cast<i128>(3 * m - 2) * m;
}
return static_cast<i128>(3) * m * m + m;
}
u64 compute_s(u32 n,
unsigned threads,
u32 /*frontier_split_x*/,
u32 /*chunk_divisor*/) {
const u32 limit = n / 5;
std::vector<i8> mu;
std::vector<i32> prefix_mu;
build_mobius(limit, mu, prefix_mu);
std::vector<IntervalTask> tasks;
tasks.reserve(static_cast<std::size_t>(2 * isqrt_u32(n) + 32));
for (u32 l = 1; l <= limit;) {
const u32 q = n / l;
const u32 r = std::min<u32>(limit, n / q);
const i32 mu_sum = prefix_mu[r] - prefix_mu[l - 1];
if (mu_sum != 0 && q > 1) {
tasks.push_back(IntervalTask{q, mu_sum});
}
l = r + 1;
}
i128 result = base_term(n);
if (threads <= 1 || tasks.size() < 2) {
for (const IntervalTask& task : tasks) {
result += static_cast<i128>(2) * task.mu_sum * calc_smart(task.q);
}
} else {
std::atomic<std::size_t> next_task(0);
std::vector<i128> partial(static_cast<std::size_t>(threads), 0);
std::vector<std::thread> pool;
pool.reserve(static_cast<std::size_t>(threads));
auto worker = [&](unsigned tid) {
i128 local = 0;
while (true) {
const std::size_t idx = next_task.fetch_add(1, std::memory_order_relaxed);
if (idx >= tasks.size()) {
break;
}
const IntervalTask& task = tasks[idx];
local += static_cast<i128>(2) * task.mu_sum * calc_smart(task.q);
}
partial[static_cast<std::size_t>(tid)] = local;
};
for (unsigned t = 0; t < threads; ++t) {
pool.emplace_back(worker, t);
}
for (auto& th : pool) {
th.join();
}
for (const i128 v : partial) {
result += v;
}
}
const i128 total = static_cast<i128>(2) * result +
static_cast<i128>(n) * static_cast<i128>(n + 1) / 2;
if (total < 0 || total > static_cast<i128>(std::numeric_limits<u64>::max())) {
throw std::overflow_error("Result does not fit in u64");
}
return static_cast<u64>(total);
}
bool run_validation() {
struct Check {
u32 n;
u64 expected;
};
const std::vector<Check> checks = {
{1, 1},
{10, 221},
{100, 39'826},
{1'000, 5'893'024},
};
bool ok = true;
for (const Check& check : checks) {
const u64 got = compute_s(check.n, 1, 128, 256);
std::cout << "Validation S(" << check.n << ") = " << got;
if (got == check.expected) {
std::cout << " [PASS]";
} else {
std::cout << " [FAIL] expected " << check.expected;
ok = false;
}
std::cout << '\n';
}
constexpr u32 kBruteN = 100;
const u64 brute = brute_force_s(kBruteN);
const u64 fast = compute_s(kBruteN, 1, 128, 256);
std::cout << "Validation brute-vs-fast S(" << kBruteN << "): fast=" << fast
<< ", brute=" << brute;
if (fast == brute) {
std::cout << " [PASS]";
} else {
std::cout << " [FAIL]";
ok = false;
}
std::cout << '\n';
return ok;
}
} // namespace
int main(int argc, char** argv) {
Options options;
if (!parse_options(argc, argv, options)) {
std::cerr << "Usage: ./Euler433 [--n=<positive-int>] [--threads=<positive-int>] "
"[--single-thread] [--no-validate] [--split=<positive-int>] "
"[--chunk-div=<positive-int>]\n";
return 1;
}
#if !defined(__OPTIMIZE__)
std::cout << "Build warning: compiled without optimization flags. "
"Use -O3 for significantly faster runtime.\n";
#endif
if (options.run_validation && !run_validation()) {
return 1;
}
const unsigned threads = choose_thread_count(
options.allow_multithreading, options.requested_threads, options.n / 2);
std::cout << "Using " << threads << " thread(s)\n";
if (threads > 1) {
// Kept for CLI compatibility with previous scheduler-tuned versions.
std::cout << "Scheduler: split_x=" << options.frontier_split_x
<< ", chunk_div=" << options.chunk_divisor << '\n';
}
const auto start = std::chrono::steady_clock::now();
const u64 result = compute_s(
options.n, threads, options.frontier_split_x, options.chunk_divisor);
const auto end = std::chrono::steady_clock::now();
const std::chrono::duration<double> elapsed = end - start;
std::cout << "S(" << options.n << ") = " << result << '\n';
std::cout << "Elapsed: " << elapsed.count() << " seconds\n";
return 0;
}
Python
import math
def solve():
n = 5000000
limit = n // 5
def isqrt(x):
r = int(math.isqrt(x))
while (r+1)*(r+1) <= x: r += 1
while r*r > x: r -= 1
return r
# Build Mobius
mu = [0]*(limit+1); mu[1] = 1
lp = [0]*(limit+1); primes = []
for i in range(2, limit+1):
if lp[i] == 0:
lp[i] = i; primes.append(i); mu[i] = -1
for p in primes:
if i*p > limit: break
lp[i*p] = p
if p == lp[i]: mu[i*p] = 0; break
mu[i*p] = -mu[i]
prefix_mu = [0]*(limit+1)
for i in range(1, limit+1): prefix_mu[i] = prefix_mu[i-1] + mu[i]
def floor_div(a, b):
return a // b if a >= 0 else -((-a + b - 1) // b)
def floor_sum_pair(p, q, m, s1, s2):
r1 = r2 = 0; sign = 1
while True:
t = floor_div(m, q); m -= t*q
r1 += sign*t*s1; r2 += sign*t*s2
t = floor_div(p, q); p -= t*q
r1 += sign*t*s1*(s1+1)//2; r2 += sign*t*s2*(s2+1)//2
if p == 0: return r1, r2
t = (p*s1 + m) // q; r1 += sign*s1*t; s1 = t
t = (p*s2 + m) // q; r2 += sign*s2*t; s2 = t
p, q = q, p; m = -m - 1; sign = -sign
def calc_smart(val):
mxy = isqrt(val); r1 = r2 = 0
for x in range(2, mxy+1):
for y in range(1, x):
sp = val // (x + y)
r1 += sp*(sp-1)//2
s1 = val//x - sp
s2 = max(0, mxy - sp) if sp <= mxy else 0
b1, b2 = floor_sum_pair(-x, y, val - sp*x, s1, s2)
r1 += b1
if sp <= mxy:
r2 += sp*(sp-1)//2; r2 += b2
else:
r2 += mxy*(mxy-1)//2
return 2*r1 - r2
def base_term(nn):
m = nn // 2
return (3*m - 2)*m if nn % 2 == 0 else 3*m*m + m
# Build tasks
tasks = []; l = 1
while l <= limit:
q = n // l; r = min(limit, n // q)
ms = prefix_mu[r] - prefix_mu[l-1]
if ms != 0 and q > 1: tasks.append((q, ms))
l = r + 1
result = base_term(n)
for q, ms in tasks: result += 2 * ms * calc_smart(q)
total = 2 * result + n * (n + 1) // 2
return str(total)
if __name__ == '__main__':
print(solve())
Java
import java.util.ArrayList;
import java.util.List;
import java.util.stream.IntStream;
import java.util.stream.LongStream;
public class Euler433 {
static long floorDiv(long a, long b) {
return (a >= 0) ? (a / b) : -(((-a) + b - 1) / b);
}
static class FloorPair {
long sum1;
long sum2;
}
static FloorPair floorSumPair(long p, long q, long m, long s1, long s2) {
FloorPair result = new FloorPair();
long sign = 1;
while (true) {
long t = floorDiv(m, q);
m -= t * q;
result.sum1 += sign * t * s1;
result.sum2 += sign * t * s2;
t = floorDiv(p, q);
p -= t * q;
result.sum1 += sign * t * s1 * (s1 + 1) / 2;
result.sum2 += sign * t * s2 * (s2 + 1) / 2;
if (p == 0)
return result;
t = floorDiv(p * s1 + m, q);
result.sum1 += sign * s1 * t;
s1 = t;
t = floorDiv(p * s2 + m, q);
result.sum2 += sign * s2 * t;
s2 = t;
long tempP = p;
p = q;
q = tempP;
m = -m - 1;
sign = -sign;
}
}
static long calcSmart(int val) {
int maxXY = (int) Math.sqrt(val);
long result1 = 0;
long result2 = 0;
for (int x = 2; x <= maxXY; x++) {
long vOverX = val / x;
for (int y = 1; y < x; y++) {
int split = val / (x + y);
result1 += (long) split * (split - 1) / 2;
long s1 = vOverX - split;
long s2 = (split <= maxXY) ? (maxXY - split) : 0;
FloorPair block = floorSumPair(-x, y, val - (long) split * x, s1, s2);
result1 += block.sum1;
if (split <= maxXY) {
result2 += (long) split * (split - 1) / 2;
result2 += block.sum2;
} else {
result2 += (long) maxXY * (maxXY - 1) / 2;
}
}
}
return 2 * result1 - result2;
}
static class MobiusData {
byte[] mu;
int[] prefix;
}
static MobiusData buildMobius(int limit) {
MobiusData data = new MobiusData();
data.mu = new byte[limit + 1];
data.prefix = new int[limit + 1];
int[] lp = new int[limit + 1];
List<Integer> primes = new ArrayList<>();
data.mu[1] = 1;
for (int i = 2; i <= limit; i++) {
if (lp[i] == 0) {
lp[i] = i;
primes.add(i);
data.mu[i] = -1;
}
for (int p : primes) {
long v = (long) i * p;
if (v > limit)
break;
lp[(int) v] = p;
if (p == lp[i]) {
data.mu[(int) v] = 0;
break;
}
data.mu[(int) v] = (byte) -data.mu[i];
}
}
for (int i = 1; i <= limit; i++) {
data.prefix[i] = data.prefix[i - 1] + data.mu[i];
}
return data;
}
static long baseTerm(int n) {
long m = n / 2;
if ((n & 1) == 0) {
return (3 * m - 2) * m;
}
return 3 * m * m + m;
}
static class IntervalTask {
int q;
int muSum;
IntervalTask(int q, int muSum) {
this.q = q;
this.muSum = muSum;
}
}
static long computeS(int n) {
int limit = n / 5;
MobiusData data = buildMobius(limit);
List<IntervalTask> tasks = new ArrayList<>();
int l = 1;
while (l <= limit) {
int q = n / l;
int r = Math.min(limit, n / q);
int muSum = data.prefix[r] - data.prefix[l - 1];
if (muSum != 0 && q > 1) {
tasks.add(new IntervalTask(q, muSum));
}
l = r + 1;
}
long result = baseTerm(n);
long parallelSum = tasks.parallelStream().mapToLong(task -> {
return 2L * task.muSum * calcSmart(task.q);
}).sum();
result += parallelSum;
return 2 * result + (long) n * (n + 1) / 2;
}
public static String solve() {
return Long.toString(computeS(5000000));
}
public static void main(String[] args) {
System.out.println(solve());
}
}