Problem 1011: Modified Euclidean Algorithm

Project Euler Problem 1011 Solution

Reduce the sum of modified Euclidean values to weighted quotient pairs and evaluate them with self-similar lookup tables in O(N log N) time. Exact implementations are available in C++, Python and Java.

Detailed mathematical approach

Problem Summary

For positive integers $a,b$, one step of the modified Euclidean algorithm divides the larger number by the smaller one and replaces the larger number by the integer quotient, not by the remainder. The process stops as soon as one entry equals $1$, and $f(a,b)$ is the other entry. We must evaluate

$$E(N)=\sum_{1\le a,b\lt N}f(a,b),\qquad N=3\,000\,000.$$

The statement gives $f(123,456)=3$, $E(10)=343$ and $E(100)=269288$. ^[1]

Mathematical Approach

For $N\ge2$, the idea is to count many starting pairs together. First separate the boundary and diagonal, then group the remaining pairs by the quotient of their first division. Finally, compute the repeated values in those groups from a small set of initial values. The weights count how often each value occurs; the table avoids computing that value again.

1. Basic properties

The classical algorithm keeps the remainder; this variant keeps the quotient. ^[2] A step maps $(a,b)$ to $(a,\lfloor b/a\rfloor)$ when $a\le b$, and to $(\lfloor a/b\rfloor,b)$ when $a\gt b$. Four facts follow directly from this rule.

Termination. While both entries are at least $2$, the divided entry at least halves, because $\lfloor b/a\rfloor\le b/2$ for $a\ge2$. Every quotient is at least $1$, so no entry reaches $0$. The product of the two entries at least halves at each step, so at most $\log_2(ab)$ steps occur.
Symmetry. For unequal entries, swapping them commutes with the division rule. If the entries become equal, either choice of entry to divide leaves the same surviving value. Thus $f(a,b)=f(b,a)$.
Boundary values. If an entry is already $1$, nothing is divided: $f(1,b)=b$ and $f(a,1)=a$. For $a\ge2$, one step turns $(a,a)$ into $(a,1)$, so $f(a,a)=a$.
No growth. Entries never increase, since a quotient never exceeds the number divided. Hence $f(a,b)\le\max(a,b)$.

For example, $(123,456)\to(123,3)\to(41,3)\to(13,3)\to(4,3)\to(1,3)$, so $f(123,456)=3$.

2. The first division

Let $T=\sum_{c=1}^{N-1}c=N(N-1)/2$. The pairs containing a $1$ contribute $2T-1$: the row $a=1$ and the column $b=1$ each sum to $T$, and $(1,1)$, whose value is $1$, was counted twice. The diagonal $2\le a=b\lt N$ contributes $T-1$. By symmetry, the remaining pairs split into two equal halves, and we keep the half with $a\lt b$:

$$E(N)=3T-2+2\sum_{2\le a\lt b\lt N}f(a,b).$$

For $2\le a\lt b$, the first step replaces $b$ by $q=\lfloor b/a\rfloor\ge1$, so $f(a,b)=f(a,q)$. Two cases arise.

If $q=1$, that is $a\lt b\lt 2a$, the pair becomes $(a,1)$ and $f(a,b)=a$. For fixed $a$ there are $\min(2a-1,N-1)-a$ such values of $b$.
If $q\ge2$, division with remainder shows that exactly the integers $b$ with $aq\le b\le aq+a-1$ have quotient $q$. ^[3] Intersecting this block with $b\le N-1$ leaves $\min(a,N-aq)$ values whenever $aq\lt N$, and all of them exceed $a$.

Therefore

$$E(N)=3T-2+2\bigl(Q_1+S\bigr),\qquad Q_1=\sum_{a=2}^{N-1}a\bigl(\min(2a-1,N-1)-a\bigr),$$

$$S=\sum_{\substack{a,q\ge2\\ aq\lt N}}\min(a,N-aq)\,f(a,q).$$

The quotient-one sum $Q_1$ costs one pass over $a$. The whole difficulty now lies in $S$, which has about $N\ln N$ terms.

3. A self-similar table for the smaller entry

In $S$, put $m=\min(a,q)$ and $x=\max(a,q)$. Then $2\le m\le x$ and $mx\lt N$, so $m^2\le N-1$ and $m\le r=\lfloor\sqrt{N-1}\rfloor$. For fixed $m$, define $v_m(x)=f(x,m)$. If $x\gt m$, the first step divides $x$ by $m$. If $x=m$, the value is $m=f(1,m)$. In both cases

$$v_m(x)=v_m\bigl(\lfloor x/m\rfloor\bigr)\quad(x\ge m),\qquad v_m(1)=m.$$

The values $v_m(x)$ with $2\le x\lt m$ are obtained by running the algorithm. All indices with the same quotient $k=\lfloor x/m\rfloor$ form the block $km\le x\le km+m-1$, and the whole block copies the single earlier value $v_m(k)$. Since $k\le x/2$, that value is already known when the block is reached.

There is also a digit interpretation of this recurrence. Repeatedly replacing $x$ by $\lfloor x/m\rfloor$ removes its last base-$m$ digit until only its leading digit $d\in\{1,\ldots,m-1\}$ remains. Consequently $v_m(x)=f(d,m)$. For instance, the table below starts from $v_3(1)=3$ and $v_3(2)=2$. The intervals $3\ldots5$ and $9\ldots17$ start with digit $1$ in base $3$, whereas $6\ldots8$ and $18\ldots26$ start with digit $2$. This explains the repeated blocks without treating their values as a coincidence.

$x$	$v_3(x)$
$1$	$3$
$2$	$2$
$3\le x\le5$	$3$
$6\le x\le8$	$2$
$9\le x\le17$	$3$
$18\le x\le26$	$2$

A pair $\{m,x\}$ with $x\gt m$ comes from the two ordered pairs $(a,q)=(m,x)$ and $(a,q)=(x,m)$, with multiplicities $\min(m,N-mx)$ and $\min(x,N-mx)$. By symmetry both have the value $v_m(x)$. For $x=m$ there is only one ordered pair. Hence

$$S=\sum_{m=2}^{r}\ \sum_{x=m}^{\lfloor (N-1)/m\rfloor}w_m(x)\,v_m(x),\qquad w_m(x)=\min(m,N-mx)+[x\ne m]\,\min(x,N-mx).$$

Here $[x\ne m]$ equals $1$ when $x\ne m$ and $0$ otherwise. Every term of $S$ is now read from a table that is filled by copying, not by running the algorithm again.

4. Worked example: $E(10)$

For $N=10$ we have $T=45$, so $3T-2=133$. The quotient-one sum is

$$Q_1=2\cdot1+3\cdot2+4\cdot3+5\cdot4+6\cdot3+7\cdot2+8\cdot1+9\cdot0=80.$$

Here $r=3$. For $m=2$, the indices $x=2,3,4$ all have $v_2(x)=2$, with weights $w_2(2)=2$, $w_2(3)=2+3=5$ and $w_2(4)=2+2=4$, so they contribute $2\cdot11=22$. For $m=3$, only $x=3$ occurs, with $w_3(3)=\min(3,1)=1$ and $v_3(3)=3$. Thus $S=25$ and

$$E(10)=133+2(80+25)=343.$$

5. Why every contribution is counted correctly

The boundary, the remaining diagonal, the quotient-one pairs and the pairs contributing to $S$ are disjoint and exhaust the original square. Every pair in $S$ has exactly one representation with $m\le x$; the weight combines its two orientations only when they are distinct. Finally, the table is correct by induction on $x$: its initial values are evaluated directly, and every later value refers to $\lfloor x/m\rfloor<x$, which has already been computed. Together these facts prove the formula for $E(N)$, independently of the numerical checks.

How the Code Works

modified_euclid simulates the definition directly. It fills the short range $2\le x\lt m$ of every table and drives the brute-force checks.

solve first computes $T$ and $Q_1$ in one pass and finds $r$ as the largest integer with $r^2\lt N$. Worker threads then share the values $m=2,\ldots,r$. Each worker claims the next $m$ with an atomic fetch_add, so every $m$ is processed exactly once. ^[4] A small $m$ carries far more work, roughly $N/m$ table entries, so this dynamic assignment balances the threads better than fixed ranges would.

Each worker owns one table value and enlarges it only when the current range $\lfloor(N-1)/m\rfloor$ does not fit. A worker receives increasing values of $m$, so its first allocation is its largest. For each $m$, the worker sets value[1] to $m$, fills $2\le x\lt m$ directly, and sweeps the blocks $k=1,2,\ldots$. In block $k$ it copies value[k] into every index $x$ with $\lfloor x/m\rfloor=k$, up to $\lfloor(N-1)/m\rfloor$, and adds $w_m(x)\,v_m(x)$ to its sum.

A crude bound shows why wide integers are used. From $f(a,b)\le N-1$ we only know $E(N)\le(N-1)^3$, and $(N-1)^3$ exceeds $2^{64}$ for the target $N$. The C++ therefore keeps its partial sums in unsigned 128-bit integers, a compiler extension supported by GCC and Clang. ^[5] Each single product $w_m(x)\,v_m(x)$ is below $2N^2$ and fits in 64 bits. The result is printed by a short decimal conversion routine, and the elapsed time goes to standard error.

Integer bounds and the three implementations

A sharper bound explains why the Java port can use long. If both starting entries are at least $2$, their minimum cannot increase before the last division. That last division creates a $1$ and leaves the previous minimum as the result, so $f(a,b)\le\min(a,b)$. Put $M=N-1$. The sum of $\min(a,b)$ over the full $M\times M$ square is $\sum_{j=1}^{M}j^2$: for each level $k$, exactly $(M-k+1)^2$ pairs have both entries at least $k$. Correcting the row and column containing $1$ adds $M(M-1)$, giving

$$E(N)\le\frac{M(M+1)(2M+1)}6+M(M-1),\qquad M=N-1.$$

At the target $N$, this upper bound is below $2^{63}-1$, so the nonnegative accumulated contributions fit in signed 64-bit integers. The C++ code keeps its conservative 128-bit accumulators. Python uses arbitrary-precision integers and processes the values of $m$ sequentially; Java uses worker threads and Math.addExact/Math.multiplyExact for the large sums. The thread-specific checks listed below describe the C++ and Java versions; the Python version performs the corresponding serial checks. ^[8]

Complexity and Verification

Computing $T$ and $Q_1$ takes $O(N)$ time. The direct part runs the algorithm on $m-2$ pairs for each $m$, which is $(r-1)(r-2)/2=O(N)$ runs of $O(\log N)$ steps. The block sweep makes $\lfloor(N-1)/m\rfloor-m+1$ updates for each $m$, so its total is

$$\sum_{m=2}^{r}\Bigl(\Bigl\lfloor\frac{N-1}{m}\Bigr\rfloor-m+1\Bigr)=\tfrac12N\ln N+\Bigl(\gamma-\tfrac32\Bigr)N+O\bigl(\sqrt N\bigr),$$

because the harmonic numbers satisfy $H_r=\ln r+\gamma+O(1/r)$ and $r^2=N+O(\sqrt N)$. ^[6] The running time is therefore $O(N\log N)$. For $N=3\,000\,000$ we get $r=1732$: the sweep makes $19\,603\,698$ updates and the direct part $1\,497\,315$ runs. The largest table has $\lfloor(N-1)/2\rfloor+1=1\,500\,000$ entries of 32 bits, about 6 MB. The result does not depend on how the values of $m$ are distributed among the threads, because every $m$ is processed once and integer addition is exact.

The 6 MB figure describes the largest single C++/Java table, not the whole process: every active worker retains its own table. With a fixed worker count, total table storage is $O(N)$. Python uses a list of object references, so its memory use is not four bytes per entry. These representation differences change resource use, but not the recurrence or the exact sum.

Before the main computation, the program runs these exact checks: ^[7]

$f(123,456)=f(456,123)=3$.
$E(10)=343$ and $E(100)=269288$, each by brute force and by the fast method.
The single-threaded fast method equals a brute-force double loop for every $2\le N\le300$.
The multithreaded fast method equals brute force for $N=1000$, $2023$ and $4096$.
One thread and the full thread count give the same value for $N=200003$.

A separate full-size computation checked the regrouping independently. It evaluated $S$ directly over all $39\,205\,665$ pairs $(a,q)$ with $a,q\ge2$ and $aq\lt N$, running the modified algorithm on each pair, without the grouping by $m$ and without the block tables. It agreed with the fast method at $N=3\,000\,000$, as did the Python and Java ports. This comparison is additional validation, separate from the packaged checks.

Footnotes and References

The references below supply the background facts used above. The reduction to weighted quotient pairs, the self-similar tables and the operation counts are derived explicitly in this article.

Project Euler 1011 — Modified Euclidean Algorithm. The official statement defines the quotient-replacing step and gives the checkpoints for f(123, 456), E(10) and E(100). These examples test the implementation without disclosing the requested final value.
Euclid's Elements, Book VII, Proposition 2. D. E. Joyce's online edition, Clark University. Euclid finds the greatest common measure by repeatedly subtracting the smaller number from the larger, which amounts to keeping the remainder. The variant in this problem keeps the quotient instead; Section 1 shows why it still terminates.
Concrete Mathematics. R. L. Graham, D. E. Knuth and O. Patashnik, Concrete Mathematics, second edition, Chapter 3, §3.1 (floors and ceilings) and §3.4 (the binary operation mod). Writing b = a⌊b/a⌋ + (b mod a) with 0 ≤ b mod a < a identifies the block of integers that share one quotient, as used in Sections 2 and 3.
cppreference — std::atomic<T>::fetch_add. The operation atomically adds to the stored value and returns the previous value. Each worker uses it to claim the next m, so no value is skipped or processed twice.
GCC manual — 128-bit Integers. The unsigned __int128 extension holds the partial sums, since the a priori bound (N−1)³ does not fit in 64 bits.
NIST Digital Library of Mathematical Functions — §5.4, equation 5.4.14. It expresses the harmonic number H_n through the digamma function as ψ(n+1) + γ, and the asymptotic expansion 5.11.2 of ψ then gives H_n = ln n + γ + O(1/n). This yields the operation count in the complexity section.
C++ — Euler1011.cpp. The linked immutable C++ revision contains the weighted quotient-pair sum, the worker threads and the checks that run before the main computation. It is the implementation record for this article. The separate full-size comparison described above is additional validation, not a test packaged in that source file.
Java Math — exact integer arithmetic documents the overflow checks used by the Java port. Python — numeric types documents unlimited integer precision. The bound above proves that the target computation fits in Java's signed 64-bit range; the checked operations additionally detect overflow if the computation is changed.

Problem 1011 source code

C++

#include <algorithm>
#include <atomic>
#include <chrono>
#include <cstdint>
#include <cstdlib>
#include <iostream>
#include <pthread.h>
#include <stdexcept>
#include <string>
#include <thread>
#include <vector>

namespace {

using u32 = std::uint32_t;
using u64 = std::uint64_t;
using u128 = unsigned __int128;

constexpr u32 TARGET = 3'000'000;

u32 modified_euclid(u32 a, u32 b) {
    while (a != 1 && b != 1) {
        if (a <= b) {
            b /= a;
        } else {
            a /= b;
        }
    }
    return a == 1 ? b : a;
}

u128 brute_force(const u32 n) {
    u128 total = 0;
    for (u32 a = 1; a < n; ++a) {
        for (u32 b = 1; b < n; ++b) total += modified_euclid(a, b);
    }
    return total;
}

struct Task {
    u32 n = 0;
    u32 root = 0;
    std::atomic<u32>* next = nullptr;
    u128 total = 0;
};

// Pairs a,q >= 2 with aq < n carry weight #{b < n : floor(b/a) = q} = min(a, n-aq);
// for m = min(a,q) and x >= m, f(x,m) = f(floor(x/m),m).
void* pair_worker(void* argument) {
    Task& task = *static_cast<Task*>(argument);
    const u32 n = task.n;
    std::vector<u32> value;
    for (u32 m = task.next->fetch_add(1); m <= task.root; m = task.next->fetch_add(1)) {
        const u32 limit = (n - 1) / m;
        if (value.size() <= limit) value.resize(limit + 1);
        value[1] = m;
        for (u32 x = 2; x < m; ++x) value[x] = modified_euclid(x, m);

        u128 sum = 0;
        for (u32 k = 1, x = m; x <= limit; ++k) {
            const u32 current = value[k];
            for (const u32 end = std::min(x + m - 1, limit); x <= end; ++x) {
                value[x] = current;
                const u32 rest = n - x * m;
                u64 weight = std::min(m, rest);
                if (x != m) weight += std::min(x, rest);
                sum += weight * current;
            }
        }
        task.total += sum;
    }
    return nullptr;
}

void require(const bool condition, const std::string& description) {
    if (!condition) throw std::runtime_error("Check failed: " + description);
}

u128 solve(const u32 n, const unsigned thread_count) {
    const u128 triangle = static_cast<u128>(n) * (n - 1) / 2;
    u128 quotient_one = 0;
    for (u64 a = 2; a < n; ++a) quotient_one += a * (std::min<u64>(2 * a - 1, n - 1) - a);

    u32 root = 1;
    while (static_cast<u64>(root + 1) * (root + 1) < n) ++root;

    std::atomic<u32> next{2};
    std::vector<Task> tasks(thread_count);
    std::vector<pthread_t> threads(thread_count);
    for (unsigned t = 0; t < thread_count; ++t) {
        tasks[t].n = n;
        tasks[t].root = root;
        tasks[t].next = &next;
        require(pthread_create(&threads[t], nullptr, pair_worker, &tasks[t]) == 0, "pthread_create");
    }
    u128 paired = 0;
    for (unsigned t = 0; t < thread_count; ++t) {
        require(pthread_join(threads[t], nullptr) == 0, "pthread_join");
        paired += tasks[t].total;
    }
    return 3 * triangle - 2 + 2 * (quotient_one + paired);
}

std::string to_string(u128 value) {
    std::string digits;
    do {
        digits.push_back(static_cast<char>('0' + value % 10));
        value /= 10;
    } while (value != 0);
    std::reverse(digits.begin(), digits.end());
    return digits;
}

void run_tests(const unsigned thread_count) {
    require(modified_euclid(123, 456) == 3 && modified_euclid(456, 123) == 3, "f(123,456) = 3");
    require(brute_force(10) == 343 && solve(10, 1) == 343, "E(10) = 343");
    require(brute_force(100) == 269288 && solve(100, thread_count) == 269288, "E(100) = 269288");
    for (u32 n = 2; n <= 300; ++n) {
        require(solve(n, 1) == brute_force(n), "brute force n=" + std::to_string(n));
    }
    for (const u32 n : {1000U, 2023U, 4096U}) {
        require(solve(n, thread_count) == brute_force(n), "brute force n=" + std::to_string(n));
    }
    require(solve(200'003, 1) == solve(200'003, thread_count), "thread consistency");
    std::cout << "All checks passed.\n";
}

}  // namespace

int main() {
    try {
        const unsigned thread_count = std::min(16U, std::max(1U, std::thread::hardware_concurrency()));
        run_tests(thread_count);
        const auto start = std::chrono::steady_clock::now();
        const u128 answer = solve(TARGET, thread_count);
        const std::chrono::duration<double> elapsed = std::chrono::steady_clock::now() - start;
        std::cout << to_string(answer) << '\n';
        std::cerr << "Computed in " << elapsed.count() << "s with " << thread_count << " threads.\n";
    } catch (const std::exception& error) {
        std::cerr << error.what() << '\n';
        return EXIT_FAILURE;
    }
    return EXIT_SUCCESS;
}

Python

#!/usr/bin/env python3
"""Project Euler Problem 1011 - Modified Euclidean Algorithm."""

import sys
import time

TARGET = 3_000_000


def modified_euclid(a, b):
    while a != 1 and b != 1:
        if a <= b:
            b //= a
        else:
            a //= b
    return b if a == 1 else a


def brute_force(n):
    total = 0
    for a in range(1, n):
        for b in range(1, n):
            total += modified_euclid(a, b)
    return total


def pair_sum(n, m, value):
    # Pairs a,q >= 2 with aq < n carry weight #{b < n : floor(b/a) = q} = min(a, n-aq);
    # for m = min(a,q) and x >= m, f(x,m) = f(floor(x/m),m).
    limit = (n - 1) // m
    value[1] = m
    for x in range(2, m):
        value[x] = modified_euclid(x, m)

    total = 0
    k = 1
    x = m
    while x <= limit:
        current = value[k]
        end = min(x + m - 1, limit)
        while x <= end:
            value[x] = current
            rest = n - x * m
            weight = min(m, rest)
            if x != m:
                weight += min(x, rest)
            total += weight * current
            x += 1
        k += 1
    return total


def solve(n):
    triangle = n * (n - 1) // 2
    quotient_one = 0
    for a in range(2, n):
        quotient_one += a * (min(2 * a - 1, n - 1) - a)

    root = 1
    while (root + 1) * (root + 1) < n:
        root += 1

    # The C++ version shares m = 2..root among threads; this port runs them in order.
    value = [0] * ((n - 1) // 2 + 1)
    paired = 0
    for m in range(2, root + 1):
        paired += pair_sum(n, m, value)
    return 3 * triangle - 2 + 2 * (quotient_one + paired)


def require(condition, description):
    if not condition:
        raise AssertionError("Check failed: " + description)


def run_tests():
    require(modified_euclid(123, 456) == 3 and modified_euclid(456, 123) == 3, "f(123,456) = 3")
    require(brute_force(10) == 343 and solve(10) == 343, "E(10) = 343")
    require(brute_force(100) == 269_288 and solve(100) == 269_288, "E(100) = 269288")
    for n in range(2, 301):
        require(solve(n) == brute_force(n), f"brute force n={n}")
    for n in (1000, 2023, 4096):
        require(solve(n) == brute_force(n), f"brute force n={n}")
    print("All checks passed.")


def main():
    run_tests()
    start = time.perf_counter()
    answer = solve(TARGET)
    elapsed = time.perf_counter() - start
    print(answer)
    print(f"Computed in {elapsed:.3f}s.", file=sys.stderr)


if __name__ == "__main__":
    try:
        main()
    except AssertionError as error:
        print(error, file=sys.stderr)
        sys.exit(1)

Java

import java.util.ArrayList;
import java.util.List;
import java.util.Locale;
import java.util.concurrent.ExecutionException;
import java.util.concurrent.ExecutorService;
import java.util.concurrent.Executors;
import java.util.concurrent.Future;
import java.util.concurrent.atomic.AtomicInteger;

public class Euler1011 {
    private static final int TARGET = 3_000_000;

    private static int modifiedEuclid(int a, int b) {
        while (a != 1 && b != 1) {
            if (a <= b) {
                b /= a;
            } else {
                a /= b;
            }
        }
        return a == 1 ? b : a;
    }

    private static long bruteForce(int n) {
        long total = 0;
        for (int a = 1; a < n; ++a) {
            for (int b = 1; b < n; ++b) total += modifiedEuclid(a, b);
        }
        return total;
    }

    // Pairs a,q >= 2 with aq < n carry weight #{b < n : floor(b/a) = q} = min(a, n-aq);
    // for m = min(a,q) and x >= m, f(x,m) = f(floor(x/m),m).
    private static long pairWorker(int n, int root, AtomicInteger next) {
        int[] value = new int[0];
        long total = 0;
        for (int m = next.getAndIncrement(); m <= root; m = next.getAndIncrement()) {
            int limit = (n - 1) / m;
            if (value.length <= limit) value = new int[limit + 1];
            value[1] = m;
            for (int x = 2; x < m; ++x) value[x] = modifiedEuclid(x, m);

            long sum = 0;
            for (int k = 1, x = m; x <= limit; ++k) {
                int current = value[k];
                for (int end = Math.min(x + m - 1, limit); x <= end; ++x) {
                    value[x] = current;
                    int rest = n - x * m;
                    long weight = Math.min(m, rest);
                    if (x != m) weight += Math.min(x, rest);
                    sum = Math.addExact(sum, weight * current);
                }
            }
            total = Math.addExact(total, sum);
        }
        return total;
    }

    private static void require(boolean condition, String description) {
        if (!condition) throw new IllegalStateException("Check failed: " + description);
    }

    // Checked arithmetic replaces the C++ 128-bit accumulators: an overflow throws instead of wrapping.
    private static long solve(int n, int threadCount) throws InterruptedException, ExecutionException {
        long triangle = (long) n * (n - 1) / 2;
        long quotientOne = 0;
        for (long a = 2; a < n; ++a) {
            quotientOne = Math.addExact(quotientOne, a * (Math.min(2 * a - 1, n - 1) - a));
        }

        int root = 1;
        while ((long) (root + 1) * (root + 1) < n) ++root;
        final int lastM = root;

        AtomicInteger next = new AtomicInteger(2);
        ExecutorService pool = Executors.newFixedThreadPool(threadCount);
        long paired = 0;
        try {
            List<Future<Long>> parts = new ArrayList<>();
            for (int t = 0; t < threadCount; ++t) {
                parts.add(pool.submit(() -> pairWorker(n, lastM, next)));
            }
            for (Future<Long> part : parts) paired = Math.addExact(paired, part.get());
        } finally {
            pool.shutdown();
        }
        long doubled = Math.multiplyExact(2L, Math.addExact(quotientOne, paired));
        return Math.addExact(Math.multiplyExact(3L, triangle) - 2, doubled);
    }

    private static void runTests(int threadCount) throws InterruptedException, ExecutionException {
        require(modifiedEuclid(123, 456) == 3 && modifiedEuclid(456, 123) == 3, "f(123,456) = 3");
        require(bruteForce(10) == 343 && solve(10, 1) == 343, "E(10) = 343");
        require(bruteForce(100) == 269_288 && solve(100, threadCount) == 269_288, "E(100) = 269288");
        for (int n = 2; n <= 300; ++n) {
            require(solve(n, 1) == bruteForce(n), "brute force n=" + n);
        }
        for (int n : new int[] {1000, 2023, 4096}) {
            require(solve(n, threadCount) == bruteForce(n), "brute force n=" + n);
        }
        require(solve(200_003, 1) == solve(200_003, threadCount), "thread consistency");
        System.out.println("All checks passed.");
    }

    public static void main(String[] args) {
        try {
            int threadCount = Math.min(16, Math.max(1, Runtime.getRuntime().availableProcessors()));
            runTests(threadCount);
            long start = System.nanoTime();
            long answer = solve(TARGET, threadCount);
            double elapsed = (System.nanoTime() - start) / 1e9;
            System.out.println(answer);
            System.err.println(String.format(Locale.ROOT, "Computed in %.3fs with %d threads.", elapsed, threadCount));
        } catch (IllegalStateException | ArithmeticException | InterruptedException | ExecutionException error) {
            System.err.println(error.getMessage());
            System.exit(1);
        }
    }
}

\(x\)	\(v_3(x)\)
\(1\)	\(3\)
\(2\)	\(2\)
\(3\le x\le5\)	\(3\)
\(6\le x\le8\)	\(2\)
\(9\le x\le17\)	\(3\)
\(18\le x\le26\)	\(2\)