Problem 870: Stone Game IV

Project Euler Problem 870 Solution

EulerSolve provides an optimized solution for Project Euler Problem 870, Stone Game IV, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary Stone Game IV studies a one-heap take-away game controlled by a real parameter $r$. On the opening turn, a player may remove any positive number of stones except the whole heap. After a move of size $x$, the next player may remove at most $r x$ stones. For each fixed $r$, some starting heap sizes are losing positions, meaning the player to move has no winning strategy. These losing positions form an increasing sequence, and that sequence changes only when $r$ crosses certain rational thresholds. The problem asks for the $123456$th value in the increasing threshold sequence obtained by starting from $r_1=1$ and repeatedly jumping to the next larger transition. Mathematical Approach The implementation does not search whole game trees. Instead it tracks how the losing-position sequence depends on $r$, then extracts the next threshold where that sequence must change. Step 1: Define the Smallest Winning Opening Move For a fixed $r$, let $w_r(n)$ be the smallest opening move that wins from a heap of size $n$. If no winning opening move exists, define $w_r(n)=n$. With this convention, $n$ is a losing position exactly when $$w_r(n)=n.$$ The recursion comes directly from the move rule. If the current player removes $k$ stones, the opponent faces $n-k$ stones and may remove at most $r k$....

Detailed mathematical approach

Problem Summary

Stone Game IV studies a one-heap take-away game controlled by a real parameter $r$. On the opening turn, a player may remove any positive number of stones except the whole heap. After a move of size $x$, the next player may remove at most $r x$ stones. For each fixed $r$, some starting heap sizes are losing positions, meaning the player to move has no winning strategy. These losing positions form an increasing sequence, and that sequence changes only when $r$ crosses certain rational thresholds. The problem asks for the $123456$th value in the increasing threshold sequence obtained by starting from $r_1=1$ and repeatedly jumping to the next larger transition.

Mathematical Approach

The implementation does not search whole game trees. Instead it tracks how the losing-position sequence depends on $r$, then extracts the next threshold where that sequence must change.

Step 1: Define the Smallest Winning Opening Move

For a fixed $r$, let $w_r(n)$ be the smallest opening move that wins from a heap of size $n$. If no winning opening move exists, define $w_r(n)=n$. With this convention, $n$ is a losing position exactly when

$$w_r(n)=n.$$

The recursion comes directly from the move rule. If the current player removes $k$ stones, the opponent faces $n-k$ stones and may remove at most $r k$. That move is winning exactly when the opponent has no winning reply of size at most $r k$. Since the opponent's smallest winning reply from $n-k$ stones is $w_r(n-k)$, we obtain

$$w_r(1)=1,\qquad w_r(n)=\min\left\{1\le k\le n:\; r k<w_r(n-k)\right\}\quad (n\ge 2).$$

The value $k=n$ is a sentinel: it means every legal opening move $1\le k\le n-1$ fails, so $n$ is losing.

Step 2: Extract the Losing Positions

Let

$$L_1=1<L_2=2<L_3<L_4<\cdots$$

be the increasing sequence of losing starting heaps. If a move of size $k$ leaves $L_j$ stones, then that move is winning precisely when

$$r k<L_j,$$

because $L_j$ itself has no winning reply smaller than $L_j$. Therefore the profitable gaps after $L_j$ are exactly the earlier losing positions whose size is still strictly below $L_j/r$.

Step 3: Recurrence for the Losing Sequence

Suppose $L_1,\dots,L_{i-1}$ are already known. Define $d_i$ as the largest lag such that

$$\frac{L_{i-1}}{L_{i-d_i}}\le r.$$

Equivalently, $L_{i-d_i}$ is the smallest earlier losing position that is at least $L_{i-1}/r$. Then the next losing position is

$$L_i=L_{i-1}+L_{i-d_i}.$$

Why does this work? Every smaller earlier losing position $L_t$ with $L_t<L_{i-1}/r$ gives a winning move from $L_{i-1}+L_t$: remove $L_t$ stones and leave the losing heap $L_{i-1}$, while the reply bound $rL_t$ is still too small to unlock a winning response. The first earlier losing position that fails this inequality is exactly the gap where the winning block ends, so adding that gap produces the next losing heap.

Step 4: Critical Transition Fractions

For a fixed $r$, the lag $d_i$ stays unchanged until one more earlier losing position becomes admissible. The critical value for step $i$ is therefore

$$c_i=\frac{L_{i-1}}{L_{i-(d_i+1)}},$$

whenever $d_i<i-1$. Crossing $c_i$ makes the next smaller denominator available, so the recurrence changes. Hence the next global transition above $r$ is

$$T(r)=\min\left\{c_i:\; d_i<i-1,\ c_i>r\right\}.$$

This is the key observation behind the solver: build the losing sequence for the current threshold, record every candidate fraction just above $r$, and keep the smallest one.

Worked Example: $r=2$

Start with $L_1=1$ and $L_2=2$.

For $L_3$, the largest admissible lag is $d_3=2$ because

$$\frac{L_2}{L_1}=\frac{2}{1}=2\le 2.$$

$$L_3=L_2+L_1=3.$$

For $L_4$, we have

$$\frac{L_3}{L_2}=\frac{3}{2}\le 2,\qquad \frac{L_3}{L_1}=3>2,$$

so $d_4=2$ and

$$L_4=L_3+L_2=5.$$

Next,

$$\frac{L_4}{L_3}=\frac{5}{3}\le 2,\qquad \frac{L_4}{L_2}=\frac{5}{2}>2,$$

which gives

$$L_5=L_4+L_3=8.$$

Continuing in the same way yields

$$1,2,3,5,8,13,\dots,$$

so the losing positions are the Fibonacci numbers. The transition candidates visible in these first steps are

$$c_4=3,\qquad c_5=\frac{5}{2},\qquad c_6=\frac{8}{3},\dots$$

The smallest candidate strictly above $2$ is $5/2$, therefore

$$T(2)=\frac{5}{2}.$$

Likewise, when $r=1$, the recurrence doubles every term and the losing positions become $1,2,4,8,\dots$, so the first transition is $T(1)=2$.

How the Code Works

The implementation stores the current threshold as an exact reduced fraction. For one transition evaluation, it generates the losing-position sequence only up to a finite horizon $M$. A monotone lag pointer is advanced while the inequality $L_{i-1}/L_{i-d}\le r$ remains true, so each new term is produced in amortized constant time. At the same moment, the implementation also forms the critical fraction $L_{i-1}/L_{i-(d+1)}$ and keeps the smallest candidate that is strictly larger than the current threshold.

All fraction comparisons are done by exact cross-multiplication, avoiding floating-point error. Sequence growth is protected by overflow checks before every addition. Because the decisive transition might occur beyond the initial horizon, the implementation repeats the computation with larger and larger horizons until the same best candidate reappears far enough before the end of the computed region. That stabilization test makes the returned fraction reliable without pretending that a short finite prefix is automatically sufficient.

Finally, the transition map is iterated from

$$r_1=1,\qquad r_{m+1}=T(r_m),$$

until $m=123456$. The C++, Python, and Java solutions all expose this same exact recurrence-based computation, and built-in checkpoints verify early values such as $r_2=2$ and $r_{22}=145/23$.

Complexity Analysis

For a fixed horizon $M$, one transition evaluation performs $O(M)$ arithmetic operations and stores $O(M)$ sequence values. The lag pointer only moves forward, so the inner search is amortized linear rather than quadratic. If $M_{\mathrm{eff}}(r)$ is the final stabilized horizon needed for threshold $r$, then a single transition costs $O(M_{\mathrm{eff}}(r))$ time and $O(M_{\mathrm{eff}}(r))$ memory. Reaching the target index therefore costs

$$O\left(\sum_{m=1}^{123455} M_{\mathrm{eff}}(r_m)\right)$$

time and

$$O\left(\max_m M_{\mathrm{eff}}(r_m)\right)$$

memory. In practice, the method is efficient because it follows only the recurrence of losing positions and the exact transition fractions, not the full game tree.

Footnotes and References

Problem page: https://projecteuler.net/problem=870
Combinatorial game theory: Wikipedia - Combinatorial game theory
Fibonacci Nim: Wikipedia - Fibonacci Nim
Rational number: Wikipedia - Rational number
Recurrence relation: Wikipedia - Recurrence relation

Problem 870 source code

C++

#include <algorithm>
#include <cstdint>
#include <cstdlib>
#include <iomanip>
#include <iostream>
#include <limits>
#include <numeric>
#include <string>
#include <vector>

namespace {

struct Fraction {
    std::uint64_t num = 0;
    std::uint64_t den = 1;
};

Fraction reduce(Fraction f) {
    std::uint64_t g = std::gcd(f.num, f.den);
    f.num /= g;
    f.den /= g;
    return f;
}

bool less_frac(const Fraction& a, const Fraction& b) {
    return static_cast<__int128>(a.num) * b.den < static_cast<__int128>(b.num) * a.den;
}

bool greater_frac(const Fraction& a, const Fraction& b) {
    return static_cast<__int128>(a.num) * b.den > static_cast<__int128>(b.num) * a.den;
}

Fraction next_transition(const Fraction& r, int base_m, int margin, int& used_m) {
    auto compute_best = [&](int M, int& end_i, bool& overflow, int& best_i) {
        std::vector<std::uint64_t> p(static_cast<std::size_t>(M) + 1, 0);
        p[1] = 1;
        p[2] = 2;
        int d = 1;
        Fraction best{0, 1};
        bool best_set = false;
        best_i = -1;
        end_i = 2;
        overflow = false;

        for (int i = 3; i <= M; ++i) {
            while (d + 1 <= i - 1) {
                const std::uint64_t left = p[static_cast<std::size_t>(i - 1)];
                const std::uint64_t right = p[static_cast<std::size_t>(i - (d + 1))];
                if (static_cast<__int128>(left) * r.den <=
                    static_cast<__int128>(r.num) * right) {
                    ++d;
                } else {
                    break;
                }
            }

            if (i - (d + 1) >= 1) {
                Fraction cand{p[static_cast<std::size_t>(i - 1)],
                              p[static_cast<std::size_t>(i - (d + 1))]};
                if (greater_frac(cand, r) && (!best_set || less_frac(cand, best))) {
                    best = cand;
                    best_set = true;
                    best_i = i;
                }
            }

            const std::uint64_t a = p[static_cast<std::size_t>(i - 1)];
            const std::uint64_t b = p[static_cast<std::size_t>(i - d)];
            if (a > std::numeric_limits<std::uint64_t>::max() - b) {
                overflow = true;
                end_i = i - 1;
                break;
            }
            p[static_cast<std::size_t>(i)] = a + b;
            end_i = i;
        }
        if (!best_set) {
            best = Fraction{0, 1};
        }
        return best;
    };

    const int max_m = 1 << 20;
    Fraction prev_best{0, 1};
    bool have_prev = false;
    int M = base_m;
    while (true) {
        int end_i = 2;
        bool overflow = false;
        int best_i = -1;
        Fraction best = compute_best(M, end_i, overflow, best_i);

        if (overflow || M >= max_m) {
            used_m = end_i;
            return reduce(best);
        }
        if (have_prev && best.num == prev_best.num && best.den == prev_best.den &&
            best_i < end_i - margin) {
            used_m = end_i;
            return reduce(best);
        }
        prev_best = best;
        have_prev = true;
        M *= 2;
    }
}

std::vector<std::uint64_t> compute_p_brut(const Fraction& r, int nmax) {
    std::vector<int> g(nmax + 1, 0);
    g[0] = 1'000'000'000;
    g[1] = 1;
    for (int n = 2; n <= nmax; ++n) {
        for (int k = 1; k <= n; ++k) {
            if (static_cast<__int128>(k) * r.num <
                static_cast<__int128>(g[n - k]) * r.den) {
                g[n] = k;
                break;
            }
        }
        if (g[n] == 0) g[n] = n;
    }
    std::vector<std::uint64_t> p;
    for (int n = 1; n <= nmax; ++n) {
        if (g[n] == n) p.push_back(static_cast<std::uint64_t>(n));
    }
    return p;
}

std::vector<std::uint64_t> compute_p_recur(const Fraction& r, int nmax) {
    std::vector<std::uint64_t> p;
    p.reserve(64);
    p.push_back(1);
    p.push_back(2);
    int d = 1;
    std::vector<std::uint64_t> work;
    work.reserve(256);
    work.push_back(0);
    work.push_back(1);
    work.push_back(2);

    while (true) {
        int i = static_cast<int>(work.size());
        while (d + 1 <= i - 1) {
            std::uint64_t left = work[static_cast<std::size_t>(i - 1)];
            std::uint64_t right = work[static_cast<std::size_t>(i - (d + 1))];
            if (static_cast<__int128>(left) * r.den <=
                static_cast<__int128>(r.num) * right) {
                ++d;
            } else {
                break;
            }
        }
        std::uint64_t next = work[static_cast<std::size_t>(i - 1)] +
                             work[static_cast<std::size_t>(i - d)];
        if (next > static_cast<std::uint64_t>(nmax)) break;
        work.push_back(next);
        p.push_back(next);
    }
    return p;
}

void validate() {
    const int base_m = 512;
    const int margin = 32;
    int used_m = 0;
    Fraction r{1, 1};
    Fraction t2{0, 1};
    Fraction t22{0, 1};
    for (int i = 2; i <= 22; ++i) {
        r = next_transition(r, base_m, margin, used_m);
        if (i == 2) t2 = r;
        if (i == 22) t22 = r;
    }
    if (!(t2.num == 2 && t2.den == 1)) {
        std::cerr << "Validation failed for T(2).\n";
        std::exit(1);
    }
    if (!(t22.num == 145 && t22.den == 23)) {
        std::cerr << "Validation failed for T(22).\n";
        std::exit(1);
    }

    const Fraction checks[] = {
        Fraction{1, 1}, Fraction{2, 1}, Fraction{5, 2},
        Fraction{3, 1}, Fraction{7, 2}, Fraction{11, 3},
    };
    for (const auto& rr : checks) {
        auto p_brut = compute_p_brut(rr, 200);
        auto p_rec = compute_p_recur(rr, 200);
        if (p_brut != p_rec) {
            std::cerr << "Validation failed for P-positions at r="
                      << rr.num << "/" << rr.den << ".\n";
            std::exit(1);
        }
    }
}

}

int main(int argc, char** argv) {
    int target = 123456;
    int base_m = 6000;
    bool do_validate = true;
    bool target_set = false;
    bool base_set = false;

    for (int i = 1; i < argc; ++i) {
        std::string arg = argv[i];
        if (arg == "--no-validate") {
            do_validate = false;
            continue;
        }
        char* end = nullptr;
        long val = std::strtol(arg.c_str(), &end, 10);
        if (end && *end == '\0') {
            if (!target_set) {
                target = std::max(1L, val);
                target_set = true;
            } else if (!base_set) {
                base_m = std::max(128L, val);
                base_set = true;
            }
        }
    }

    if (do_validate) {
        validate();
    }

    Fraction ans{1, 1};
    if (target == 1) {
        ans = Fraction{1, 1};
    } else {
        Fraction r{1, 1};
        const int margin = 64;
        int used_m = base_m;
        for (int i = 2; i <= target; ++i) {
            ans = next_transition(r, base_m, margin, used_m);
            r = ans;
        }
    }

    long double value = static_cast<long double>(ans.num) / static_cast<long double>(ans.den);
    std::cout << std::fixed << std::setprecision(10) << value << '\n';
    return 0;
}

Python

from __future__ import annotations

import re
import shutil
import subprocess
from pathlib import Path

ANSWER_RE = re.compile(r"answer\s*:\s*(.+)$", re.IGNORECASE)
EQUAL_RE = re.compile(r"=\s*(.+)$")


def parse_output(stdout: str) -> str:
    lines = [line.strip() for line in stdout.splitlines() if line.strip()]
    if not lines:
        return ""
    answers = []
    equals = []
    for line in lines:
        m1 = ANSWER_RE.search(line)
        if m1:
            answers.append(m1.group(1).strip())
        m2 = EQUAL_RE.search(line)
        if m2:
            equals.append(m2.group(1).strip())
    if answers:
        return answers[-1]
    if equals:
        return equals[-1]
    return lines[-1]


def should_skip_cpp_checkpoints(src: Path) -> bool:
    try:
        text = src.read_text(encoding="utf-8", errors="ignore")
    except OSError:
        return False
    return "--skip-checkpoints" in text


def run_cpp(binary: Path, src: Path, root: Path) -> str:
    cmd = [str(binary)]
    if should_skip_cpp_checkpoints(src):
        cmd.append("--skip-checkpoints")

    try:
        return subprocess.check_output(cmd, text=True, cwd=root)
    except subprocess.CalledProcessError:
        return subprocess.check_output(cmd, text=True, cwd=src.parent)


def solve() -> str:
    problem_id = __file__.split("Euler")[-1].split(".")[0]
    root = Path(__file__).resolve().parent.parent
    src = root / "solutionsCpp" / f"Euler{problem_id}.cpp"
    binary = root / "solutionsCpp" / f".euler{problem_id}_py_bridge"

    if not binary.exists() or src.stat().st_mtime > binary.stat().st_mtime:
        compiler = shutil.which("clang++") or shutil.which("g++")
        if not compiler:
            raise RuntimeError("No C++ compiler found (clang++/g++).")
        subprocess.check_call([compiler, "-std=c++17", "-O2", str(src), "-o", str(binary)])

    output = run_cpp(binary=binary, src=src, root=root)
    parsed = parse_output(output)
    if not parsed:
        raise RuntimeError(f"Euler{problem_id} bridge produced empty output.")
    return parsed


if __name__ == "__main__":
    print(solve())

Java

import java.nio.file.*;
import java.util.*;
import java.util.regex.*;

public class Euler870 {
    private static final Pattern ANSWER_RE = Pattern.compile("answer\\s*:\\s*(.+)$", Pattern.CASE_INSENSITIVE);
    private static final Pattern EQUAL_RE = Pattern.compile("=\\s*(.+)$");

    private static String parseOutput(String stdout) {
        String[] lines = stdout.split("\\R");
        List<String> nonEmpty = new ArrayList<>();
        for (String line : lines) {
            String t = line.trim();
            if (!t.isEmpty()) {
                nonEmpty.add(t);
            }
        }
        if (nonEmpty.isEmpty()) {
            return "";
        }

        List<String> answers = new ArrayList<>();
        List<String> equals = new ArrayList<>();
        for (String line : nonEmpty) {
            Matcher m1 = ANSWER_RE.matcher(line);
            if (m1.find()) {
                answers.add(m1.group(1).trim());
            }
            Matcher m2 = EQUAL_RE.matcher(line);
            if (m2.find()) {
                equals.add(m2.group(1).trim());
            }
        }

        if (!answers.isEmpty()) {
            return answers.get(answers.size() - 1);
        }
        if (!equals.isEmpty()) {
            return equals.get(equals.size() - 1);
        }
        return nonEmpty.get(nonEmpty.size() - 1);
    }

    private static String pickCompiler() throws Exception {
        for (String compiler : List.of("clang++", "g++")) {
            Process probe = new ProcessBuilder("bash", "-lc", "command -v " + compiler)
                    .redirectErrorStream(true)
                    .start();
            String out = new String(probe.getInputStream().readAllBytes());
            int rc = probe.waitFor();
            if (rc == 0 && !out.trim().isEmpty()) {
                return compiler;
            }
        }
        throw new RuntimeException("No C++ compiler found (clang++/g++).");
    }

    private static Path cppSource(Path root) {
        return root.resolve("solutionsCpp").resolve("Euler870.cpp");
    }

    private static boolean shouldSkipCheckpoints(Path root) {
        Path src = cppSource(root);
        try {
            String text = Files.readString(src);
            return text.contains("--skip-checkpoints");
        } catch (Exception ex) {
            return false;
        }
    }

    private static Path ensureBridgeBinary() throws Exception {
        Path root = Paths.get(System.getProperty("user.dir"));
        Path src = cppSource(root);
        Path bin = root.resolve("solutionsCpp").resolve(".euler870_java_bridge");

        boolean rebuild = Files.notExists(bin)
                || Files.getLastModifiedTime(src).compareTo(Files.getLastModifiedTime(bin)) > 0;

        if (rebuild) {
            String compiler = pickCompiler();
            Process compile = new ProcessBuilder(
                    compiler,
                    "-std=c++17",
                    "-O2",
                    src.toString(),
                    "-o",
                    bin.toString())
                    .inheritIO()
                    .start();
            if (compile.waitFor() != 0) {
                throw new RuntimeException("Failed to compile Euler870 C++ bridge.");
            }
        }

        return bin;
    }

    private static String runBridge(Path bin, Path root, Path srcDir) throws Exception {
        List<String> cmd = new ArrayList<>();
        cmd.add(bin.toString());
        if (shouldSkipCheckpoints(root)) {
            cmd.add("--skip-checkpoints");
        }

        Process first = new ProcessBuilder(cmd)
                .directory(root.toFile())
                .redirectErrorStream(true)
                .start();
        String out = new String(first.getInputStream().readAllBytes());
        int rc = first.waitFor();
        if (rc == 0) {
            return out;
        }

        Process second = new ProcessBuilder(cmd)
                .directory(srcDir.toFile())
                .redirectErrorStream(true)
                .start();
        String out2 = new String(second.getInputStream().readAllBytes());
        int rc2 = second.waitFor();
        if (rc2 == 0) {
            return out2;
        }

        throw new RuntimeException("Euler870 C++ bridge failed.\n" + out + "\n" + out2);
    }

    private static String solveViaCppBridge() throws Exception {
        Path root = Paths.get(System.getProperty("user.dir"));
        Path src = cppSource(root);
        Path bin = ensureBridgeBinary();
        String out = runBridge(bin, root, src.getParent());
        String parsed = parseOutput(out);
        if (parsed.isEmpty()) {
            throw new RuntimeException("Euler870 C++ bridge produced empty output.");
        }
        return parsed;
    }

    public static void main(String[] args) throws Exception {
        System.out.println(solveViaCppBridge());
    }
}