Problem 1003: Lonely Singles

Project Euler Problem 1003 Solution

EulerSolve provides an optimized solution for Project Euler Problem 1003, Lonely Singles, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The process starts with $n$ stones at position $0$. Positions are then processed in increasing order. If the current position $i$ contains $m$ stones, one singleton is left behind exactly when $m$ is odd; in all cases $\lfloor m/2\rfloor$ stones are sent to $i+1$ and another $\lfloor m/2\rfloor$ stones are sent to $i+3$. Only the positions where the current pile was odd matter for this problem. A singleton is lonely when no other singleton is within distance $2$, so a sad integer is one whose singleton positions are pairwise separated by at least $3$. The function $S(k)$ asks for the sum of all sad starting values whose complete singleton set is contained in $0\le i\lt k$. The target is $S(80)$, and the problem gives the checks $S(14)=159$ and $S(30)=33438$. The main challenge is that many starting values do not terminate. Even when the process keeps moving stones forever, it can still leave only finitely many singleton positions. So the solution cannot rely on simulating until all stones disappear; it must characterize the finite parity trace left by the process. Mathematical Approach Viewing the process as a carrying rule At each position the process removes pairs of stones and turns each pair at $i$ into one stone at $i+1$ and one stone at $i+3$. The possible leftover digit at $i$ is therefore only $0$ or $1$....

Detailed mathematical approach

Problem Summary

The process starts with $n$ stones at position $0$. Positions are then processed in increasing order. If the current position $i$ contains $m$ stones, one singleton is left behind exactly when $m$ is odd; in all cases $\lfloor m/2\rfloor$ stones are sent to $i+1$ and another $\lfloor m/2\rfloor$ stones are sent to $i+3$.

Only the positions where the current pile was odd matter for this problem. A singleton is lonely when no other singleton is within distance $2$, so a sad integer is one whose singleton positions are pairwise separated by at least $3$. The function $S(k)$ asks for the sum of all sad starting values whose complete singleton set is contained in $0\le i\lt k$. The target is $S(80)$, and the problem gives the checks $S(14)=159$ and $S(30)=33438$.

The main challenge is that many starting values do not terminate. Even when the process keeps moving stones forever, it can still leave only finitely many singleton positions. So the solution cannot rely on simulating until all stones disappear; it must characterize the finite parity trace left by the process.

Mathematical Approach

Viewing the process as a carrying rule

At each position the process removes pairs of stones and turns each pair at $i$ into one stone at $i+1$ and one stone at $i+3$. The possible leftover digit at $i$ is therefore only $0$ or $1$. This is analogous to carrying in a positional numeral system, except the carry rule is unusual:

$$2\cdot \text{position }i \quad\longrightarrow\quad \text{position }(i+1)+\text{position }(i+3).$$

If $s_i\in\{0,1\}$ records whether a singleton was left at $i$, then the process deterministically normalizes the starting pile $n$ into a digit sequence $s_0,s_1,s_2,\dots$. The sadness condition is then a condition on the support of this digit sequence.

An algebraic invariant for the carry rule

Attach a formal weight $x^i$ to position $i$. For the carry rule to preserve weight, one pair at position $i$ must have the same value as the two carried stones:

$$2x^i=x^{i+1}+x^{i+3}.$$

After division by $x^i$, this becomes

$$x^3+x-2=0.$$

The polynomial factors as

$$x^3+x-2=(x-1)(x^2+x+2).$$

The factor $x-1$ only tracks total mass, which is not enough here because the process may continue forever while leaving finitely many singletons. The useful finite trace is captured by the quadratic factor, so the implementation works in the ring where

$$x^2+x+2=0,\qquad\text{equivalently}\qquad x^2=-x-2.$$

Reducing each position to two coefficients

Modulo $x^2+x+2$, every power has the form

$$x^i=a_i+b_i x.$$

The first two values are $x^0=1$ and $x^1=x$. If $x^{i-1}=a+b x$, then

$$x^i=x(a+b x)=a x+b x^2=a x+b(-x-2)=-2b+(a-b)x.$$

Thus the code stores only two integers per position and advances by

$$a_i=-2b_{i-1},\qquad b_i=a_{i-1}-b_{i-1}.$$

For a finite singleton set $P$, define

$$V(P)=\sum_{i\in P}x^i=A(P)+B(P)x.$$

If $P$ is the complete finite singleton set for some starting value $n$, then its value in this quadratic quotient must equal the initial value $n$, which has no $x$-part. Therefore the required algebraic test is

$$B(P)=0,\qquad n=A(P),\qquad A(P)\gt0.$$

Why this enumerates the right objects

The left-to-right process is deterministic: once $n$ is fixed, the parity digit $s_i$ and the outgoing carry from position $i$ are fixed. Conversely, a finite digit set $P$ whose quadratic value is a positive integer is a valid finite trace for that integer under this carrying rule. The program uses this as a bijection between sad starting values counted by $S(k)$ and finite subsets $P\subseteq\{0,\dots,k-1\}$ satisfying the algebraic test and the loneliness gap.

This is also why the final answer is a sum of values $A(P)$, not a count of sets. Each accepted set $P$ represents one starting value $n=A(P)$, and $S(k)$ asks for the sum of those starting values.

The loneliness condition as an independent-set condition

A finite set $P$ of singleton positions is lonely exactly when

$$|i-j|\ge 3\qquad\text{for all distinct }i,j\in P.$$

Equivalently, in the path graph on positions $0,1,\dots,k-1$, where vertices at distance $1$ or $2$ conflict, $P$ is an independent set. The problem has now become a finite combinatorial enumeration: list gap-$3$ subsets, reduce their algebraic value to $A+Bx$, and retain only positive integer values with $B=0$.

The examples fall out directly. $\{0\}$ gives $1$. The set $\{2,5,8,13\}$ gives $68$, and $\{1,13\}$ gives $90$. These are the sad contributions below $k=14$, so $S(14)=1+68+90=159$.

Why a direct subset search is too large

Let $F(t)$ be the number of subsets of $\{0,\dots,t-1\}$ with mutual distance at least $3$. Looking at the last position gives

$$F(t)=F(t-1)+F(t-3).$$

The term $F(t-1)$ covers subsets not using the last position; the term $F(t-3)$ covers subsets using it, because the previous two positions are then forbidden. For $t=80$, this recurrence already gives $25058735850088$ candidate subsets. That is far beyond direct enumeration.

Splitting at $m=\lfloor k/2\rfloor$ changes the scale. For $k=80$, each half has length $40$, and the number of valid half-subsets is only $F(40)=5736961$. Millions of states are manageable; tens of trillions are not.

Meet-in-the-middle matching

Write the left-half value as $A_L+B_Lx$ and the right-half value as $A_R+B_Rx$. The union represents a positive integer exactly when

$$B_L+B_R=0,\qquad A_L+A_R\gt0.$$

So the right half only needs to be searchable by coefficient $B_R$. For each left subset, the required right coefficient is $-B_L$, and every compatible right subset with that coefficient contributes the starting value $A_L+A_R$.

The C++ implementation stores right entries sorted by coefficient and keeps prefix sums of their constants. That allows it to support even the general case where many right subsets have the same coefficient: binary search locates the matching block, and a second binary search discards entries whose total constant is not positive. The Python and Java ports use a lighter coefficient-to-constant table because, for the $k\le80$ instances used here, the right coefficients are unique inside each boundary class.

Handling the split boundary

The only missing condition after splitting is the distance rule across the cut. A left subset ending near $m$ can conflict with a right subset beginning near $m$. The code records exactly the information needed for this decision.

For the left half, the category is

$$0:\text{ no selected point at }m-2\text{ or }m-1,\qquad 1:\text{ last}=m-2,\qquad 2:\text{ last}=m-1.$$

For the right half, the category is

$$0:\text{ no selected point at }m\text{ or }m+1,\qquad 1:\text{ first}=m+1,\qquad 2:\text{ first}=m.$$

If the left category is $0$, there is no left-side boundary conflict. If it is $1$, only right category $2$ is too close. If it is $2$, both $m$ and $m+1$ are too close, so only right category $0$ is allowed. This three-by-three compatibility rule is exactly the cross-boundary form of $|i-j|\ge3$.

How the Code Works

build_powers constructs the pairs $(a_i,b_i)$ for $x^i=a_i+b_i x$. evaluate_positions sums those pairs for a proposed singleton set. singleton_positions directly simulates the stone process for small checks, which verifies both the algebraic representation and the interpretation of the singleton trace.

enumerate_right recursively enumerates all right-half gap-$3$ subsets and places them into the three boundary buckets. Each stored entry contains the reduced coefficient and constant. build_left_tasks creates shallow left prefixes; the C++ version uses them to distribute the deeper left search across worker threads.

combine_left_dfs finishes each left subset, computes its boundary category, and queries only compatible right buckets. A match is accepted when the coefficients cancel and the constant sum is positive. That positive constant is added to the answer because it is the represented starting value.

The checkpoint suite verifies the known traces for $n=1$, $n=68$, and $n=90$, then verifies $S(14)=159$ and $S(30)=33438$. Only after those checks does the program evaluate $S(80)$.

Complexity Analysis

The unsplit search would examine $F(k)$ gap-valid subsets, which is not feasible for $k=80$. The meet-in-the-middle method stores $O(F(k/2))$ right-half states and performs $O(F(k/2))$ left-half searches.

In the sorted-table C++ version, right-table construction costs $O(F(k/2)\log F(k/2))$, and each completed left subset performs a constant number of binary-search lookups. The map-based Python and Java ports have expected $O(F(k/2))$ table construction and expected constant-time lookups for the actual $k\le80$ range. Memory usage is $O(F(k/2))$.

Footnotes and References

Problem page: Project Euler 1003 - Lonely Singles
Algebraic number: Wikipedia - Algebraic number
Meet-in-the-middle algorithm: Wikipedia - Meet-in-the-middle
Independent set: Wikipedia - Independent set
Binary search: Wikipedia - Binary search

Problem 1003 source code

C++

#include <algorithm>
#include <array>
#include <cassert>
#include <cstdint>
#include <cstdlib>
#include <iostream>
#include <limits>
#include <string>
#include <thread>
#include <vector>

namespace {

using i64 = long long;
using i128 = __int128_t;

constexpr int kTarget = 80;
constexpr int kTaskPrefix = 24;

struct Options {
    int target_k = kTarget;
    bool run_checkpoints = true;
    unsigned threads = 0;
};

struct Power {
    i64 constant = 0;
    i64 coeff = 0;
};

struct Entry {
    i64 coeff = 0;
    i64 constant = 0;
};

struct RightTable {
    std::vector<Entry> entries;
    std::vector<i128> prefix_constant;
};

struct LeftTask {
    int pos = 0;
    int last = -1000000;
    i64 constant = 0;
    i64 coeff = 0;
};

bool entry_less(const Entry& a, const Entry& b) {
    if (a.coeff != b.coeff) {
        return a.coeff < b.coeff;
    }
    return a.constant < b.constant;
}

std::string to_string_i128(i128 value) {
    if (value == 0) {
        return "0";
    }

    bool negative = false;
    if (value < 0) {
        negative = true;
        value = -value;
    }

    std::string result;
    while (value > 0) {
        const int digit = static_cast<int>(value % 10);
        result.push_back(static_cast<char>('0' + digit));
        value /= 10;
    }
    if (negative) {
        result.push_back('-');
    }
    std::reverse(result.begin(), result.end());
    return result;
}

unsigned resolve_threads(const Options& options) {
    if (options.threads != 0) {
        return std::max(1U, options.threads);
    }
    const unsigned detected = std::thread::hardware_concurrency();
    return detected == 0 ? 1U : detected;
}

std::vector<Power> build_powers(const int k) {
    std::vector<Power> powers(std::max(2, k));
    powers[0] = {1, 0};
    powers[1] = {0, 1};
    for (int i = 2; i < k; ++i) {
        const i64 a = powers[i - 1].constant;
        const i64 b = powers[i - 1].coeff;
        powers[i] = {-2 * b, a - b};
    }
    powers.resize(k);
    return powers;
}

std::pair<i64, i64> evaluate_positions(const std::vector<int>& positions, const std::vector<Power>& powers) {
    i64 constant = 0;
    i64 coeff = 0;
    for (const int pos : positions) {
        constant += powers[pos].constant;
        coeff += powers[pos].coeff;
    }
    return {constant, coeff};
}

std::vector<int> singleton_positions(const i64 n, const int limit) {
    std::vector<i64> stones(static_cast<std::size_t>(limit + 4), 0);
    std::vector<int> positions;
    stones[0] = n;
    for (int i = 0; i < limit; ++i) {
        if ((stones[i] & 1LL) != 0) {
            positions.push_back(i);
        }
        const i64 moved = stones[i] / 2;
        stones[i + 1] += moved;
        stones[i + 3] += moved;
    }
    return positions;
}

int left_category(const int last, const int mid) {
    if (last < 0 || last <= mid - 3) {
        return 0;
    }
    if (last == mid - 2) {
        return 1;
    }
    return 2;
}

int right_category(const int first, const int mid) {
    if (first < 0 || first >= mid + 2) {
        return 0;
    }
    if (first == mid + 1) {
        return 1;
    }
    return 2;
}

bool compatible_boundary(const int left_cat, const int right_cat) {
    if (left_cat == 0) {
        return true;
    }
    if (left_cat == 1) {
        return right_cat != 2;
    }
    return right_cat == 0;
}

void enumerate_right(const int pos,
                     const int end,
                     const int mid,
                     const int last,
                     const int first,
                     const i64 constant,
                     const i64 coeff,
                     const std::vector<Power>& powers,
                     std::array<std::vector<Entry>, 3>& buckets) {
    if (pos == end) {
        buckets[static_cast<std::size_t>(right_category(first, mid))].push_back({coeff, constant});
        return;
    }

    enumerate_right(pos + 1, end, mid, last, first, constant, coeff, powers, buckets);

    if (last < 0 || pos - last >= 3) {
        const int next_first = first < 0 ? pos : first;
        enumerate_right(pos + 1,
                        end,
                        mid,
                        pos,
                        next_first,
                        constant + powers[pos].constant,
                        coeff + powers[pos].coeff,
                        powers,
                        buckets);
    }
}

void build_left_tasks(const int pos,
                      const int stop,
                      const int mid,
                      const int last,
                      const i64 constant,
                      const i64 coeff,
                      const std::vector<Power>& powers,
                      std::vector<LeftTask>& tasks) {
    if (pos == stop) {
        tasks.push_back({pos, last, constant, coeff});
        return;
    }

    build_left_tasks(pos + 1, stop, mid, last, constant, coeff, powers, tasks);

    if (last < 0 || pos - last >= 3) {
        build_left_tasks(pos + 1,
                         stop,
                         mid,
                         pos,
                         constant + powers[pos].constant,
                         coeff + powers[pos].coeff,
                         powers,
                         tasks);
    }
}

void prepare_table(std::vector<Entry>& entries, RightTable& table) {
    std::sort(entries.begin(), entries.end(), entry_less);
    table.entries.swap(entries);
    table.prefix_constant.assign(table.entries.size() + 1, 0);
    for (std::size_t i = 0; i < table.entries.size(); ++i) {
        table.prefix_constant[i + 1] = table.prefix_constant[i] + table.entries[i].constant;
    }
}

i128 query_table(const RightTable& table, const i64 target_coeff, const i64 left_constant) {
    const Entry low{target_coeff, std::numeric_limits<i64>::min()};
    const Entry high{target_coeff, std::numeric_limits<i64>::max()};

    const auto lo = std::lower_bound(table.entries.begin(), table.entries.end(), low, entry_less);
    const auto hi = std::upper_bound(table.entries.begin(), table.entries.end(), high, entry_less);
    if (lo == hi) {
        return 0;
    }

    const Entry threshold{target_coeff, -left_constant};
    const auto first_positive = std::upper_bound(lo, hi, threshold, entry_less);
    const std::size_t first_index = static_cast<std::size_t>(first_positive - table.entries.begin());
    const std::size_t end_index = static_cast<std::size_t>(hi - table.entries.begin());
    const i128 count = static_cast<i128>(end_index - first_index);
    const i128 right_sum = table.prefix_constant[end_index] - table.prefix_constant[first_index];

    return count * static_cast<i128>(left_constant) + right_sum;
}

i128 add_matches(const int left_cat,
                 const i64 left_coeff,
                 const i64 left_constant,
                 const std::array<RightTable, 3>& tables) {
    i128 total = 0;
    const i64 target_coeff = -left_coeff;
    for (int right_cat = 0; right_cat < 3; ++right_cat) {
        if (compatible_boundary(left_cat, right_cat)) {
            total += query_table(tables[static_cast<std::size_t>(right_cat)], target_coeff, left_constant);
        }
    }
    return total;
}

i128 combine_left_dfs(const int pos,
                      const int mid,
                      const int last,
                      const i64 constant,
                      const i64 coeff,
                      const std::vector<Power>& powers,
                      const std::array<RightTable, 3>& tables) {
    if (pos == mid) {
        return add_matches(left_category(last, mid), coeff, constant, tables);
    }

    i128 total = combine_left_dfs(pos + 1, mid, last, constant, coeff, powers, tables);
    if (last < 0 || pos - last >= 3) {
        total += combine_left_dfs(pos + 1,
                                  mid,
                                  pos,
                                  constant + powers[pos].constant,
                                  coeff + powers[pos].coeff,
                                  powers,
                                  tables);
    }
    return total;
}

i128 solve_s(const int k, const unsigned requested_threads) {
    if (k == 0) {
        return 0;
    }

    const int mid = k / 2;
    const std::vector<Power> powers = build_powers(k);

    std::array<std::vector<Entry>, 3> right_buckets;
    enumerate_right(mid, k, mid, -1000000, -1, 0, 0, powers, right_buckets);

    std::array<RightTable, 3> tables;
    for (int i = 0; i < 3; ++i) {
        prepare_table(right_buckets[static_cast<std::size_t>(i)], tables[static_cast<std::size_t>(i)]);
    }

    const int task_stop = std::min(mid, kTaskPrefix);
    std::vector<LeftTask> tasks;
    build_left_tasks(0, task_stop, mid, -1000000, 0, 0, powers, tasks);

    const unsigned workers =
        std::max(1U, std::min<unsigned>(requested_threads, static_cast<unsigned>(std::max<std::size_t>(1, tasks.size()))));

    if (workers == 1) {
        i128 total = 0;
        for (const LeftTask& task : tasks) {
            total += combine_left_dfs(task.pos, mid, task.last, task.constant, task.coeff, powers, tables);
        }
        return total;
    }

    std::vector<i128> partial(workers, 0);
    std::vector<std::thread> pool;
    pool.reserve(workers);
    for (unsigned worker = 0; worker < workers; ++worker) {
        pool.emplace_back([&, worker]() {
            i128 local = 0;
            for (std::size_t i = worker; i < tasks.size(); i += workers) {
                const LeftTask& task = tasks[i];
                local += combine_left_dfs(task.pos, mid, task.last, task.constant, task.coeff, powers, tables);
            }
            partial[worker] = local;
        });
    }
    for (std::thread& thread : pool) {
        thread.join();
    }

    i128 total = 0;
    for (const i128 value : partial) {
        total += value;
    }
    return total;
}

void require_checkpoint(const bool ok, const std::string& message) {
    if (!ok) {
        std::cerr << "Checkpoint failed: " << message << '\n';
        std::exit(EXIT_FAILURE);
    }
}

void run_checkpoints(const unsigned threads) {
    const std::vector<Power> powers = build_powers(kTarget);

    const std::vector<int> p68{2, 5, 8, 13};
    const std::vector<int> p90{1, 13};
    require_checkpoint(evaluate_positions({0}, powers) == std::make_pair<i64, i64>(1, 0), "n=1 polynomial");
    require_checkpoint(evaluate_positions(p68, powers) == std::make_pair<i64, i64>(68, 0), "n=68 polynomial");
    require_checkpoint(evaluate_positions(p90, powers) == std::make_pair<i64, i64>(90, 0), "n=90 polynomial");

    require_checkpoint(singleton_positions(1, 40) == std::vector<int>{0}, "n=1 trace");
    require_checkpoint(singleton_positions(68, 50) == p68, "n=68 trace");
    require_checkpoint(singleton_positions(90, 50) == p90, "n=90 trace");

    require_checkpoint(solve_s(14, 1) == 159, "S(14)");
    require_checkpoint(solve_s(30, 1) == 33438, "S(30)");
    require_checkpoint(solve_s(30, threads) == 33438, "threaded S(30)");

    std::cerr << "Validation checkpoints passed.\n";
}

void usage() {
    std::cerr << "Usage:\n"
              << "  ./Euler1003 [k] [--skip-checkpoints] [--single-thread] [--threads=N]\n";
}

Options parse_options(const int argc, char** argv) {
    Options options;
    for (int i = 1; i < argc; ++i) {
        const std::string arg(argv[i]);
        if (arg == "--skip-checkpoints") {
            options.run_checkpoints = false;
        } else if (arg == "--single-thread") {
            options.threads = 1;
        } else if (arg.rfind("--threads=", 0) == 0) {
            options.threads = static_cast<unsigned>(std::stoul(arg.substr(10)));
            if (options.threads == 0) {
                options.threads = 1;
            }
        } else if (!arg.empty() && arg[0] == '-') {
            usage();
            std::exit(EXIT_FAILURE);
        } else {
            options.target_k = std::stoi(arg);
        }
    }

    if (options.target_k < 0 || options.target_k > kTarget) {
        std::cerr << "k must satisfy 0 <= k <= " << kTarget << ".\n";
        std::exit(EXIT_FAILURE);
    }

    return options;
}

}  // namespace

int main(int argc, char** argv) {
    const Options options = parse_options(argc, argv);
    const unsigned threads = resolve_threads(options);

    if (options.run_checkpoints) {
        run_checkpoints(threads);
    }

    std::cout << to_string_i128(solve_s(options.target_k, threads)) << '\n';
    return 0;
}

Python

#!/usr/bin/env python3
"""Project Euler Problem 1003 - Lonely Singles."""

from __future__ import annotations

import argparse
import os
import sys
from dataclasses import dataclass


TARGET_K = 80
NEG_INF = -10**9


@dataclass(frozen=True)
class Power:
    constant: int
    coeff: int


@dataclass(frozen=True)
class LeftTask:
    pos: int
    last: int
    constant: int
    coeff: int


def build_powers(k: int) -> list[Power]:
    powers = [Power(0, 0) for _ in range(max(2, k))]
    powers[0] = Power(1, 0)
    powers[1] = Power(0, 1)
    for i in range(2, k):
        prev = powers[i - 1]
        powers[i] = Power(-2 * prev.coeff, prev.constant - prev.coeff)
    return powers[:k]


def evaluate_positions(positions: list[int], powers: list[Power]) -> tuple[int, int]:
    constant = 0
    coeff = 0
    for pos in positions:
        constant += powers[pos].constant
        coeff += powers[pos].coeff
    return constant, coeff


def singleton_positions(n: int, limit: int) -> list[int]:
    stones = [0] * (limit + 4)
    positions: list[int] = []
    stones[0] = n
    for i in range(limit):
        if stones[i] & 1:
            positions.append(i)
        moved = stones[i] // 2
        stones[i + 1] += moved
        stones[i + 3] += moved
    return positions


def left_category(last: int, mid: int) -> int:
    if last < 0 or last <= mid - 3:
        return 0
    if last == mid - 2:
        return 1
    return 2


def right_category(first: int, mid: int) -> int:
    if first < 0 or first >= mid + 2:
        return 0
    if first == mid + 1:
        return 1
    return 2


def compatible_boundary(left_cat: int, right_cat: int) -> bool:
    if left_cat == 0:
        return True
    if left_cat == 1:
        return right_cat != 2
    return right_cat == 0


def add_right_entry(table: dict[int, int], coeff: int, constant: int) -> None:
    if coeff in table:
        raise RuntimeError("unexpected duplicate right-side coefficient")
    table[coeff] = constant


def enumerate_right(
    pos: int,
    end: int,
    mid: int,
    last: int,
    first: int,
    constant: int,
    coeff: int,
    powers: list[Power],
    buckets: list[dict[int, int]],
) -> None:
    if pos == end:
        add_right_entry(buckets[right_category(first, mid)], coeff, constant)
        return

    enumerate_right(pos + 1, end, mid, last, first, constant, coeff, powers, buckets)

    if last < 0 or pos - last >= 3:
        power = powers[pos]
        next_first = pos if first < 0 else first
        enumerate_right(
            pos + 1,
            end,
            mid,
            pos,
            next_first,
            constant + power.constant,
            coeff + power.coeff,
            powers,
            buckets,
        )


def build_left_tasks(
    pos: int,
    stop: int,
    last: int,
    constant: int,
    coeff: int,
    powers: list[Power],
    tasks: list[LeftTask],
) -> None:
    if pos == stop:
        tasks.append(LeftTask(pos, last, constant, coeff))
        return

    build_left_tasks(pos + 1, stop, last, constant, coeff, powers, tasks)

    if last < 0 or pos - last >= 3:
        power = powers[pos]
        build_left_tasks(
            pos + 1,
            stop,
            pos,
            constant + power.constant,
            coeff + power.coeff,
            powers,
            tasks,
        )


def query_table(table: dict[int, int], target_coeff: int, left_constant: int) -> int:
    right_constant = table.get(target_coeff)
    if right_constant is None:
        return 0
    value = left_constant + right_constant
    return value if value > 0 else 0


def add_matches(
    left_cat: int,
    left_coeff: int,
    left_constant: int,
    tables: list[dict[int, int]],
) -> int:
    total = 0
    target_coeff = -left_coeff
    for right_cat, table in enumerate(tables):
        if compatible_boundary(left_cat, right_cat):
            total += query_table(table, target_coeff, left_constant)
    return total


def combine_left_dfs(
    pos: int,
    mid: int,
    last: int,
    constant: int,
    coeff: int,
    powers: list[Power],
    tables: list[dict[int, int]],
) -> int:
    if pos == mid:
        return add_matches(left_category(last, mid), coeff, constant, tables)

    total = combine_left_dfs(pos + 1, mid, last, constant, coeff, powers, tables)
    if last < 0 or pos - last >= 3:
        power = powers[pos]
        total += combine_left_dfs(
            pos + 1,
            mid,
            pos,
            constant + power.constant,
            coeff + power.coeff,
            powers,
            tables,
        )
    return total


def solve_s(k: int) -> int:
    if k == 0:
        return 0

    mid = k // 2
    powers = build_powers(k)

    right_buckets: list[dict[int, int]] = [dict(), dict(), dict()]
    enumerate_right(mid, k, mid, NEG_INF, -1, 0, 0, powers, right_buckets)

    task_stop = min(mid, 24)
    tasks: list[LeftTask] = []
    build_left_tasks(0, task_stop, NEG_INF, 0, 0, powers, tasks)

    total = 0
    for task in tasks:
        total += combine_left_dfs(
            task.pos,
            mid,
            task.last,
            task.constant,
            task.coeff,
            powers,
            right_buckets,
        )
    return total


def require_checkpoint(ok: bool, message: str) -> None:
    if not ok:
        raise SystemExit(f"Checkpoint failed: {message}")


def run_checkpoints() -> None:
    powers = build_powers(TARGET_K)

    p68 = [2, 5, 8, 13]
    p90 = [1, 13]
    require_checkpoint(evaluate_positions([0], powers) == (1, 0), "n=1 polynomial")
    require_checkpoint(evaluate_positions(p68, powers) == (68, 0), "n=68 polynomial")
    require_checkpoint(evaluate_positions(p90, powers) == (90, 0), "n=90 polynomial")

    require_checkpoint(singleton_positions(1, 40) == [0], "n=1 trace")
    require_checkpoint(singleton_positions(68, 50) == p68, "n=68 trace")
    require_checkpoint(singleton_positions(90, 50) == p90, "n=90 trace")

    require_checkpoint(solve_s(14) == 159, "S(14)")
    require_checkpoint(solve_s(30) == 33438, "S(30)")

    print("Validation checkpoints passed.", file=sys.stderr)


def parse_args() -> argparse.Namespace:
    parser = argparse.ArgumentParser(description="Solve Project Euler Problem 1003.")
    parser.add_argument("k", nargs="?", type=int, default=TARGET_K)
    parser.add_argument("--skip-checkpoints", action="store_true")
    parser.add_argument("--single-thread", action="store_true", help=argparse.SUPPRESS)
    parser.add_argument("--threads", type=int, default=os.cpu_count() or 1, help=argparse.SUPPRESS)
    args = parser.parse_args()
    if args.k < 0 or args.k > TARGET_K:
        parser.error(f"k must satisfy 0 <= k <= {TARGET_K}.")
    return args


def main() -> int:
    args = parse_args()
    if not args.skip_checkpoints:
        run_checkpoints()
    print(solve_s(args.k))
    return 0


if __name__ == "__main__":
    raise SystemExit(main())

Java

import java.util.ArrayList;
import java.util.Arrays;
import java.util.List;

public class Euler1003 {
    private static final int TARGET_K = 80;
    private static final int TASK_PREFIX = 24;
    private static final int NEG_INF = -1_000_000_000;
    private static final long EMPTY = Long.MIN_VALUE;

    private static final class Options {
        int targetK = TARGET_K;
        boolean runCheckpoints = true;
    }

    private static final class Power {
        final long constant;
        final long coeff;

        Power(long constant, long coeff) {
            this.constant = constant;
            this.coeff = coeff;
        }
    }

    private static final class LeftTask {
        final int pos;
        final int last;
        final long constant;
        final long coeff;

        LeftTask(int pos, int last, long constant, long coeff) {
            this.pos = pos;
            this.last = last;
            this.constant = constant;
            this.coeff = coeff;
        }
    }

    private static final class LongLongMap {
        private final long[] keys;
        private final long[] values;
        private final int mask;

        LongLongMap(int expectedSize) {
            int capacity = 1;
            while (capacity < Math.max(4, expectedSize * 2)) {
                capacity <<= 1;
            }
            keys = new long[capacity];
            values = new long[capacity];
            Arrays.fill(keys, EMPTY);
            mask = capacity - 1;
        }

        private static long mix(long x) {
            x ^= x >>> 33;
            x *= 0xff51afd7ed558ccdL;
            x ^= x >>> 33;
            x *= 0xc4ceb9fe1a85ec53L;
            x ^= x >>> 33;
            return x;
        }

        void put(long key, long value) {
            if (key == EMPTY) {
                throw new IllegalArgumentException("reserved key");
            }
            int slot = (int) mix(key) & mask;
            while (keys[slot] != EMPTY) {
                if (keys[slot] == key) {
                    throw new IllegalStateException("unexpected duplicate right-side coefficient");
                }
                slot = (slot + 1) & mask;
            }
            keys[slot] = key;
            values[slot] = value;
        }

        long get(long key) {
            int slot = (int) mix(key) & mask;
            while (keys[slot] != EMPTY) {
                if (keys[slot] == key) {
                    return values[slot];
                }
                slot = (slot + 1) & mask;
            }
            return EMPTY;
        }
    }

    private static final class RightBucket {
        private long[] coeffs = new long[1024];
        private long[] constants = new long[1024];
        private int size = 0;

        void add(long coeff, long constant) {
            if (size == coeffs.length) {
                int next = coeffs.length << 1;
                coeffs = Arrays.copyOf(coeffs, next);
                constants = Arrays.copyOf(constants, next);
            }
            coeffs[size] = coeff;
            constants[size] = constant;
            size++;
        }

        LongLongMap finish() {
            LongLongMap map = new LongLongMap(size);
            for (int i = 0; i < size; ++i) {
                map.put(coeffs[i], constants[i]);
            }
            coeffs = null;
            constants = null;
            return map;
        }
    }

    private static Power[] buildPowers(int k) {
        Power[] powers = new Power[Math.max(2, k)];
        powers[0] = new Power(1, 0);
        powers[1] = new Power(0, 1);
        for (int i = 2; i < k; ++i) {
            Power prev = powers[i - 1];
            powers[i] = new Power(-2L * prev.coeff, prev.constant - prev.coeff);
        }
        return Arrays.copyOf(powers, k);
    }

    private static long[] evaluatePositions(int[] positions, Power[] powers) {
        long constant = 0;
        long coeff = 0;
        for (int pos : positions) {
            constant += powers[pos].constant;
            coeff += powers[pos].coeff;
        }
        return new long[]{constant, coeff};
    }

    private static List<Integer> singletonPositions(long n, int limit) {
        long[] stones = new long[limit + 4];
        List<Integer> positions = new ArrayList<>();
        stones[0] = n;
        for (int i = 0; i < limit; ++i) {
            if ((stones[i] & 1L) != 0) {
                positions.add(i);
            }
            long moved = stones[i] / 2;
            stones[i + 1] += moved;
            stones[i + 3] += moved;
        }
        return positions;
    }

    private static int leftCategory(int last, int mid) {
        if (last < 0 || last <= mid - 3) {
            return 0;
        }
        if (last == mid - 2) {
            return 1;
        }
        return 2;
    }

    private static int rightCategory(int first, int mid) {
        if (first < 0 || first >= mid + 2) {
            return 0;
        }
        if (first == mid + 1) {
            return 1;
        }
        return 2;
    }

    private static boolean compatibleBoundary(int leftCat, int rightCat) {
        if (leftCat == 0) {
            return true;
        }
        if (leftCat == 1) {
            return rightCat != 2;
        }
        return rightCat == 0;
    }

    private static void enumerateRight(
            int pos,
            int end,
            int mid,
            int last,
            int first,
            long constant,
            long coeff,
            Power[] powers,
            RightBucket[] buckets) {
        if (pos == end) {
            buckets[rightCategory(first, mid)].add(coeff, constant);
            return;
        }

        enumerateRight(pos + 1, end, mid, last, first, constant, coeff, powers, buckets);

        if (last < 0 || pos - last >= 3) {
            Power power = powers[pos];
            int nextFirst = first < 0 ? pos : first;
            enumerateRight(
                    pos + 1,
                    end,
                    mid,
                    pos,
                    nextFirst,
                    constant + power.constant,
                    coeff + power.coeff,
                    powers,
                    buckets);
        }
    }

    private static void buildLeftTasks(
            int pos,
            int stop,
            int last,
            long constant,
            long coeff,
            Power[] powers,
            List<LeftTask> tasks) {
        if (pos == stop) {
            tasks.add(new LeftTask(pos, last, constant, coeff));
            return;
        }

        buildLeftTasks(pos + 1, stop, last, constant, coeff, powers, tasks);

        if (last < 0 || pos - last >= 3) {
            Power power = powers[pos];
            buildLeftTasks(
                    pos + 1,
                    stop,
                    pos,
                    constant + power.constant,
                    coeff + power.coeff,
                    powers,
                    tasks);
        }
    }

    private static long queryTable(LongLongMap table, long targetCoeff, long leftConstant) {
        long rightConstant = table.get(targetCoeff);
        if (rightConstant == EMPTY) {
            return 0;
        }
        long value = leftConstant + rightConstant;
        return value > 0 ? value : 0;
    }

    private static long addMatches(
            int leftCat,
            long leftCoeff,
            long leftConstant,
            LongLongMap[] tables) {
        long total = 0;
        long targetCoeff = -leftCoeff;
        for (int rightCat = 0; rightCat < 3; ++rightCat) {
            if (compatibleBoundary(leftCat, rightCat)) {
                total += queryTable(tables[rightCat], targetCoeff, leftConstant);
            }
        }
        return total;
    }

    private static long combineLeftDfs(
            int pos,
            int mid,
            int last,
            long constant,
            long coeff,
            Power[] powers,
            LongLongMap[] tables) {
        if (pos == mid) {
            return addMatches(leftCategory(last, mid), coeff, constant, tables);
        }

        long total = combineLeftDfs(pos + 1, mid, last, constant, coeff, powers, tables);
        if (last < 0 || pos - last >= 3) {
            Power power = powers[pos];
            total += combineLeftDfs(
                    pos + 1,
                    mid,
                    pos,
                    constant + power.constant,
                    coeff + power.coeff,
                    powers,
                    tables);
        }
        return total;
    }

    private static long solveS(int k) {
        if (k == 0) {
            return 0;
        }

        int mid = k / 2;
        Power[] powers = buildPowers(k);

        RightBucket[] rightBuckets = {new RightBucket(), new RightBucket(), new RightBucket()};
        enumerateRight(mid, k, mid, NEG_INF, -1, 0, 0, powers, rightBuckets);

        LongLongMap[] tables = new LongLongMap[3];
        for (int i = 0; i < 3; ++i) {
            tables[i] = rightBuckets[i].finish();
            rightBuckets[i] = null;
        }

        int taskStop = Math.min(mid, TASK_PREFIX);
        List<LeftTask> tasks = new ArrayList<>();
        buildLeftTasks(0, taskStop, NEG_INF, 0, 0, powers, tasks);

        long total = 0;
        for (LeftTask task : tasks) {
            total += combineLeftDfs(
                    task.pos,
                    mid,
                    task.last,
                    task.constant,
                    task.coeff,
                    powers,
                    tables);
        }
        return total;
    }

    private static void check(boolean ok, String message) {
        if (!ok) {
            throw new IllegalStateException("Checkpoint failed: " + message);
        }
    }

    private static void runCheckpoints() {
        Power[] powers = buildPowers(TARGET_K);

        int[] p68 = {2, 5, 8, 13};
        int[] p90 = {1, 13};
        check(Arrays.equals(evaluatePositions(new int[]{0}, powers), new long[]{1, 0}), "n=1 polynomial");
        check(Arrays.equals(evaluatePositions(p68, powers), new long[]{68, 0}), "n=68 polynomial");
        check(Arrays.equals(evaluatePositions(p90, powers), new long[]{90, 0}), "n=90 polynomial");

        check(singletonPositions(1, 40).equals(List.of(0)), "n=1 trace");
        check(singletonPositions(68, 50).equals(List.of(2, 5, 8, 13)), "n=68 trace");
        check(singletonPositions(90, 50).equals(List.of(1, 13)), "n=90 trace");

        check(solveS(14) == 159, "S(14)");
        check(solveS(30) == 33438, "S(30)");

        System.err.println("Validation checkpoints passed.");
    }

    private static void usage() {
        System.err.println("Usage:");
        System.err.println("  java Euler1003 [k] [--skip-checkpoints] [--single-thread] [--threads=N]");
    }

    private static Options parseOptions(String[] args) {
        Options options = new Options();
        for (String arg : args) {
            if (arg.equals("--skip-checkpoints")) {
                options.runCheckpoints = false;
            } else if (arg.equals("--single-thread") || arg.startsWith("--threads=")) {
                // Accepted for command-line compatibility with the C++ version.
            } else if (arg.startsWith("-")) {
                usage();
                System.exit(1);
            } else {
                options.targetK = Integer.parseInt(arg);
            }
        }
        if (options.targetK < 0 || options.targetK > TARGET_K) {
            throw new IllegalArgumentException("k must satisfy 0 <= k <= " + TARGET_K + ".");
        }
        return options;
    }

    public static void main(String[] args) {
        Options options = parseOptions(args);
        if (options.runCheckpoints) {
            runCheckpoints();
        }
        System.out.println(solveS(options.targetK));
    }
}