Problem 308: An Amazing Prime-generating Automaton

Project Euler Problem 308 Solution

EulerSolve provides an optimized solution for Project Euler Problem 308, An Amazing Prime-generating Automaton, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary The automaton in Problem 308 is a graphical encoding of Conway's prime-generating FRACTRAN program. If we simulate it literally, the number of steps explodes long before we reach the $10001$-st prime output. The key is to understand one complete candidate cycle analytically and count steps without walking through every transition. Mathematical Approach 1) State interpretation: prime exponents, not raw integers The FRACTRAN machine can be represented by exponent counts of the primes $$2,3,5,7,11,13,17,19,23,29.$$ So a machine state is equivalent to a vector of exponents, or to the integer $$N=2^a3^b5^c7^d11^e13^f17^g19^h23^i29^j.$$ The automaton “prints” a value exactly when the state is a pure power of two: $$N=2^m.$$ Then the printed number is the exponent $m$. The brute-force checkpoint in the C++ code confirms that the first outputs are $$2,3,5,7,11,13,17,19,23,29,\dots$$ with step indices $$19,69,281,710,2375,3893,8102,11361,19268,36981,\dots$$ 2) Canonical start-state for candidate $n$ The program does not test primes by trial division in the usual sense....

Detailed mathematical approach

Problem Summary

The automaton in Problem 308 is a graphical encoding of Conway's prime-generating FRACTRAN program. If we simulate it literally, the number of steps explodes long before we reach the $10001$-st prime output. The key is to understand one complete candidate cycle analytically and count steps without walking through every transition.

Mathematical Approach

1) State interpretation: prime exponents, not raw integers

The FRACTRAN machine can be represented by exponent counts of the primes

$$2,3,5,7,11,13,17,19,23,29.$$

So a machine state is equivalent to a vector of exponents, or to the integer

$$N=2^a3^b5^c7^d11^e13^f17^g19^h23^i29^j.$$

The automaton “prints” a value exactly when the state is a pure power of two:

$$N=2^m.$$

Then the printed number is the exponent $m$. The brute-force checkpoint in the C++ code confirms that the first outputs are

$$2,3,5,7,11,13,17,19,23,29,\dots$$

with step indices

$$19,69,281,710,2375,3893,8102,11361,19268,36981,\dots$$

2) Canonical start-state for candidate $n$

The program does not test primes by trial division in the usual sense. Instead, the FRACTRAN dynamics repeatedly reaches a canonical state that encodes “we are now testing candidate $n$.” In the optimized derivation, the first step index of that canonical state is denoted by

$$T_n.$$

The source code uses

$$T_2=2,$$

meaning that after two initial transitions from the seed, the machine has entered the standard start configuration for candidate $2$.

3) What one candidate cycle does

Once candidate $n$ starts, the automaton behaves like a deterministic divisor-testing routine.

Setup: it builds the control structure for the current candidate.

Middle phase: it tests divisors $d$ in descending order.

Exit: if a divisor is found, the machine moves on to candidate $n+1$; if no divisor exists, it follows a special prime branch and emits $2^n$.

This is why the whole process can be summarized by two exact arithmetic quantities:

$$\Delta(n)=T_{n+1}-T_n,$$

the total step count from the start of candidate $n$ to the start of candidate $n+1$, and

$$H(p),$$

the offset from $T_p$ to the moment where prime candidate $p$ is actually emitted.

4) Why the largest proper divisor appears

Let

$$s=\operatorname{spf}(n)$$

be the smallest prime factor of $n$. Then the first divisor encountered while scanning downward is the largest proper divisor

$$d_*=\frac{n}{s}.$$

All integers

$$d=n-1,n-2,\dots,d_*+1$$

are guaranteed to be non-divisors, and the scan stops exactly at $d_*$. That is why the composite-case formula naturally splits into “all non-divisor tests above $d_*$” plus “the final successful divisor branch at $d_*$.”

5) Exact cost formulas

The code has already flattened the entire state graph of one candidate cycle into exact transition counts.

For a general candidate $n$, the total jump to the next candidate is

$$\Delta(n)=(2n-1)+\sum_{d=d_*+1}^{n-1}\bigl(6n+2+2\lfloor n/d\rfloor\bigr)+\bigl(5n+d_*+2\lfloor n/d_*\rfloor+2\bigr).$$

The three parts are:

Setup: $2n-1$.

Rejected divisor tests: one term $6n+2+2\lfloor n/d\rfloor$ for every non-divisor above $d_*$.

Terminal composite branch: the final successful branch at $d_*$.

If $p$ is prime, there is no divisor in $2,\dots,p-1$, so the machine instead reaches the prime-output branch:

$$H(p)=(2p-1)+\sum_{d=2}^{p-1}\bigl(6p+2+2\lfloor p/d\rfloor\bigr)+(6p+2).$$

The last term $6p+2$ is the prime branch itself.

6) Accumulating global time

Once $\Delta(n)$ is known, the candidate-start times satisfy the prefix recurrence

$$T_{n+1}=T_n+\Delta(n),\qquad T_2=2.$$

If $p_k$ is the $k$-th prime, then the answer to the Project Euler problem is simply

$$T_{p_k}+H(p_k).$$

For example, candidate $2$ starts at step $T_2=2$, and

$$H(2)=17,$$

so the first prime hit occurs at

$$2+17=19.$$

Likewise, $\Delta(2)=20$, so $T_3=22$; then $H(3)=47$, giving the second prime hit at

$$22+47=69.$$

These match the brute-force checkpoints exactly.

7) Fast evaluation of the floor sums

The expensive-looking part is

$$\sum \lfloor n/d\rfloor.$$

The implementation evaluates such sums in quotient blocks: if several consecutive divisors have the same quotient $q=\lfloor n/d\rfloor$, they are added at once. This is the standard harmonic-grouping trick and makes each floor-sum sublinear instead of linear.

How the Code Works

1) Prime bound. The solver first estimates an upper bound for the $k$-th prime using $n(\log n+\log\log n)$-type asymptotics, then enlarges it if necessary.

2) Linear sieve. A sieve produces both the prime list and the smallest-prime-factor table $\operatorname{spf}(n)$.

3) Floor-sum helpers. floor_sum_upto() and floor_sum_range() implement quotient-block summation.

4) Closed formulas. candidate_delta() computes $\Delta(n)$, and prime_hit_offset() computes $H(p)$.

5) Prefix walk. The main loop advances $T_n$ candidate by candidate until it reaches the target prime $p_k$, then returns $T_{p_k}+H(p_k)$.

6) Strong validation. The C++ code also contains a literal automaton simulator and checks that the first 10 prime hits and several small target indices agree with the closed formulas.

Complexity Analysis

The sieve costs $O(B)$ time and memory up to the chosen prime bound $B$. Each candidate then needs only a constant amount of algebra plus a few harmonic floor-sums, each taking about $O(\sqrt n)$ block steps. This is enormously faster than literal FRACTRAN simulation, which would need to traverse about $1.5\times10^{15}$ transitions for the final answer.

Problem 308 source code

C++

#include <algorithm>
#include <cmath>
#include <cstdint>
#include <exception>
#include <iostream>
#include <limits>
#include <stdexcept>
#include <string>
#include <utility>
#include <vector>

namespace {

using u64 = std::uint64_t;
using u32 = std::uint32_t;
using u128 = __uint128_t;

struct Options {
    u64 target_index = 10001ULL;
    bool run_checkpoints = true;
};

bool parse_u64_after_prefix(const std::string& arg,
                            const std::string& prefix,
                            u64& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }

    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }

    u64 parsed = 0ULL;
    for (const char c : tail) {
        if (c < '0' || c > '9') {
            return false;
        }
        const u64 digit = static_cast<u64>(c - '0');
        if (parsed > (std::numeric_limits<u64>::max() - digit) / 10ULL) {
            throw std::overflow_error("Argument overflow");
        }
        parsed = parsed * 10ULL + digit;
    }

    value = parsed;
    return true;
}

bool parse_arguments(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;
            continue;
        }
        if (parse_u64_after_prefix(arg, "--target=", options.target_index)) {
            continue;
        }

        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }

    return true;
}

std::string to_string_u128(u128 value) {
    if (value == 0U) {
        return "0";
    }

    std::string digits;
    while (value > 0U) {
        const int digit = static_cast<int>(value % 10U);
        digits.push_back(static_cast<char>('0' + digit));
        value /= 10U;
    }

    std::reverse(digits.begin(), digits.end());
    return digits;
}

u64 estimate_upper_bound_for_nth_prime(const u64 n) {
    if (n < 6ULL) {
        return 15ULL;
    }

    const long double nd = static_cast<long double>(n);
    const long double estimate =
        nd * (std::log(nd) + std::log(std::log(nd))) + 32.0L;

    if (!std::isfinite(static_cast<double>(estimate)) || estimate <= 0.0L) {
        throw std::runtime_error("Prime bound estimate failed");
    }

    const u64 bound = static_cast<u64>(std::ceil(estimate));
    return std::max<u64>(bound, 15ULL);
}

struct SieveData {
    std::vector<u32> primes;
    std::vector<u32> spf;
};

SieveData build_linear_sieve(const u64 limit) {
    if (limit > static_cast<u64>(std::numeric_limits<u32>::max())) {
        throw std::runtime_error("Limit too large for 32-bit SPF table");
    }

    SieveData out;
    out.spf.assign(static_cast<std::size_t>(limit) + 1ULL, 0U);
    out.primes.reserve(static_cast<std::size_t>(limit / 10ULL + 16ULL));

    for (u32 i = 2U; i <= limit; ++i) {
        if (out.spf[static_cast<std::size_t>(i)] == 0U) {
            out.spf[static_cast<std::size_t>(i)] = i;
            out.primes.push_back(i);
        }

        for (const u32 p : out.primes) {
            const u64 composite = static_cast<u64>(p) * static_cast<u64>(i);
            if (composite > limit) {
                break;
            }

            out.spf[static_cast<std::size_t>(composite)] = p;
            if (p == out.spf[static_cast<std::size_t>(i)]) {
                break;
            }
        }
    }

    return out;
}

SieveData build_sieve_with_prime_count(const u64 target_index) {
    u64 limit = estimate_upper_bound_for_nth_prime(target_index);

    while (true) {
        SieveData sieve = build_linear_sieve(limit);
        if (sieve.primes.size() >= target_index) {
            return sieve;
        }
        if (limit > std::numeric_limits<u64>::max() / 2ULL) {
            throw std::runtime_error("Prime sieve limit overflow");
        }
        limit *= 2ULL;
    }
}

u128 floor_sum_upto(const u64 n, u64 r) {
    if (r == 0ULL) {
        return 0U;
    }
    if (r > n) {
        r = n;
    }

    u128 sum = 0U;
    u64 left = 1ULL;
    while (left <= r) {
        const u64 q = n / left;
        u64 right = n / q;
        if (right > r) {
            right = r;
        }

        sum += static_cast<u128>(q) * static_cast<u128>(right - left + 1ULL);
        left = right + 1ULL;
    }

    return sum;
}

u128 floor_sum_range(const u64 n, const u64 l, const u64 r) {
    if (l > r || r == 0ULL) {
        return 0U;
    }
    return floor_sum_upto(n, r) - floor_sum_upto(n, l - 1ULL);
}

u128 candidate_delta_to_next_start(const u64 n, const u64 spf_n) {
    const u64 max_proper_divisor = n / spf_n;

    const u64 non_divisor_count =
        (n > max_proper_divisor + 1ULL) ? (n - max_proper_divisor - 1ULL) : 0ULL;
    const u128 non_divisor_floor_sum =
        floor_sum_range(n, max_proper_divisor + 1ULL, n - 1ULL);

    // Candidate processing decomposition:
    // setup in S11, all non-divisor tests, and terminal divisor branch.
    const u128 setup = static_cast<u128>(2ULL * n - 1ULL);
    const u128 non_divisor_total =
        static_cast<u128>(6ULL * n + 2ULL) * static_cast<u128>(non_divisor_count) +
        2U * non_divisor_floor_sum;
    const u128 terminal =
        static_cast<u128>(5ULL * n + max_proper_divisor + 2ULL * (n / max_proper_divisor) +
                          2ULL);

    return setup + non_divisor_total + terminal;
}

u128 prime_hit_offset_from_start(const u64 prime_candidate) {
    const u128 non_divisor_floor_sum = floor_sum_range(prime_candidate, 2ULL, prime_candidate - 1ULL);

    const u128 setup = static_cast<u128>(2ULL * prime_candidate - 1ULL);
    const u128 non_divisor_total =
        static_cast<u128>(6ULL * prime_candidate + 2ULL) *
            static_cast<u128>(prime_candidate - 2ULL) +
        2U * non_divisor_floor_sum;
    const u128 prime_branch_to_hit = static_cast<u128>(6ULL * prime_candidate + 2ULL);

    return setup + non_divisor_total + prime_branch_to_hit;
}

u128 solve_iterations_to_prime_index(const u64 target_index) {
    if (target_index == 0ULL) {
        throw std::invalid_argument("--target must be positive");
    }

    const SieveData sieve = build_sieve_with_prime_count(target_index);
    const u64 prime_target =
        static_cast<u64>(sieve.primes[static_cast<std::size_t>(target_index - 1ULL)]);

    // T_n is the first step index where the machine is in canonical start-state
    // for candidate n. T_2 = 2 from seed 2 after applying rules 12 and 14 once each.
    u128 start_step = 2U;
    u64 found_primes = 0ULL;

    for (u64 n = 2ULL; n <= prime_target; ++n) {
        const u64 spf_n = static_cast<u64>(sieve.spf[static_cast<std::size_t>(n)]);
        const bool is_prime = (spf_n == n);

        if (is_prime) {
            const u128 hit_step = start_step + prime_hit_offset_from_start(n);
            ++found_primes;
            if (found_primes == target_index) {
                return hit_step;
            }
        }

        start_step += candidate_delta_to_next_start(n, spf_n);
    }

    throw std::runtime_error("Prime target was not reached");
}

bool is_prime_small(const u64 x) {
    if (x < 2ULL) {
        return false;
    }
    if ((x % 2ULL) == 0ULL) {
        return x == 2ULL;
    }
    for (u64 d = 3ULL; d * d <= x; d += 2ULL) {
        if ((x % d) == 0ULL) {
            return false;
        }
    }
    return true;
}

struct FractranState {
    u64 a = 1ULL;  // exponent of 2
    u64 b = 0ULL;  // exponent of 3
    u64 c = 0ULL;  // exponent of 5
    u64 d = 0ULL;  // exponent of 7
    u64 e = 0ULL;  // exponent of 11
    u64 f = 0ULL;  // exponent of 13
    u64 g = 0ULL;  // exponent of 17
    u64 h = 0ULL;  // exponent of 19
    u64 i = 0ULL;  // exponent of 23
    u64 j = 0ULL;  // exponent of 29
};

bool is_pure_power_of_two(const FractranState& s) {
    return s.b == 0ULL && s.c == 0ULL && s.d == 0ULL && s.e == 0ULL && s.f == 0ULL &&
           s.g == 0ULL && s.h == 0ULL && s.i == 0ULL && s.j == 0ULL;
}

void apply_one_fractran_step(FractranState& s) {
    if (s.d > 0ULL && s.f > 0ULL) {
        --s.d;
        --s.f;
        ++s.g;
        return;
    }
    if (s.c > 0ULL && s.g > 0ULL) {
        ++s.a;
        ++s.b;
        --s.c;
        ++s.f;
        --s.g;
        return;
    }
    if (s.b > 0ULL && s.g > 0ULL) {
        --s.b;
        --s.g;
        ++s.h;
        return;
    }
    if (s.a > 0ULL && s.h > 0ULL) {
        --s.a;
        --s.h;
        ++s.i;
        return;
    }
    if (s.b > 0ULL && s.e > 0ULL) {
        --s.b;
        --s.e;
        ++s.j;
        return;
    }
    if (s.j > 0ULL) {
        ++s.d;
        ++s.e;
        --s.j;
        return;
    }
    if (s.i > 0ULL) {
        ++s.c;
        ++s.h;
        --s.i;
        return;
    }
    if (s.h > 0ULL) {
        ++s.d;
        ++s.e;
        --s.h;
        return;
    }
    if (s.g > 0ULL) {
        --s.g;
        return;
    }
    if (s.f > 0ULL) {
        ++s.e;
        --s.f;
        return;
    }
    if (s.e > 0ULL) {
        --s.e;
        ++s.f;
        return;
    }
    if (s.a > 0ULL) {
        --s.a;
        ++s.b;
        ++s.c;
        return;
    }
    if (s.d > 0ULL) {
        --s.d;
        return;
    }

    ++s.c;
    ++s.e;
}

std::vector<std::pair<u64, u64>> brute_prime_hits(const std::size_t count) {
    std::vector<std::pair<u64, u64>> hits;
    hits.reserve(count);

    FractranState state;
    u64 steps = 0ULL;

    while (hits.size() < count) {
        if (is_pure_power_of_two(state) && is_prime_small(state.a)) {
            hits.push_back({state.a, steps});
        }

        apply_one_fractran_step(state);
        ++steps;
        if (steps > 2000000000ULL) {
            throw std::runtime_error("Brute-force checkpoint exceeded step cap");
        }
    }

    return hits;
}

bool check_equal_u64(const u64 actual, const u64 expected, const std::string& label) {
    if (actual == expected) {
        return true;
    }
    std::cerr << "Checkpoint failed for " << label << ": expected " << expected
              << ", got " << actual << '\n';
    return false;
}

bool check_equal_u128(const u128 actual,
                      const u128 expected,
                      const std::string& label) {
    if (actual == expected) {
        return true;
    }
    std::cerr << "Checkpoint failed for " << label << ": expected "
              << to_string_u128(expected) << ", got " << to_string_u128(actual)
              << '\n';
    return false;
}

bool run_checkpoints() {
    const std::vector<std::pair<u64, u64>> brute10 = brute_prime_hits(10U);
    const std::vector<std::pair<u64, u64>> expected10 = {
        {2ULL, 19ULL},   {3ULL, 69ULL},    {5ULL, 281ULL},   {7ULL, 710ULL},
        {11ULL, 2375ULL}, {13ULL, 3893ULL}, {17ULL, 8102ULL}, {19ULL, 11361ULL},
        {23ULL, 19268ULL}, {29ULL, 36981ULL}};

    if (brute10.size() != expected10.size()) {
        std::cerr << "Checkpoint failed: brute sequence size mismatch\n";
        return false;
    }
    for (std::size_t i = 0; i < expected10.size(); ++i) {
        if (!check_equal_u64(brute10[i].first, expected10[i].first,
                             "prime exponent #" + std::to_string(i + 1U))) {
            return false;
        }
        if (!check_equal_u64(brute10[i].second, expected10[i].second,
                             "iteration count #" + std::to_string(i + 1U))) {
            return false;
        }
    }

    const std::vector<std::pair<u64, u64>> brute30 = brute_prime_hits(30U);
    const std::vector<u64> targets = {1ULL, 2ULL, 3ULL, 10ULL, 20ULL, 30ULL};
    for (const u64 target : targets) {
        const u128 formula = solve_iterations_to_prime_index(target);
        const u128 expected =
            static_cast<u128>(brute30[static_cast<std::size_t>(target - 1ULL)].second);
        if (!check_equal_u128(formula, expected,
                              "formula-vs-brute for target=" + std::to_string(target))) {
            return false;
        }
    }

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    try {
        Options options;
        if (!parse_arguments(argc, argv, options)) {
            return 1;
        }
        if (options.target_index == 0ULL) {
            std::cerr << "--target must be positive.\n";
            return 1;
        }

        if (options.run_checkpoints && !run_checkpoints()) {
            return 1;
        }

        const u128 answer = solve_iterations_to_prime_index(options.target_index);
        std::cout << to_string_u128(answer) << '\n';
    } catch (const std::exception& ex) {
        std::cerr << "Error: " << ex.what() << '\n';
        return 1;
    }

    return 0;
}

Python

import math

def solve():
    target_index = 10001

    def estimate_upper(n):
        if n < 6:
            return 15
        nd = float(n)
        return int(math.ceil(nd * (math.log(nd) + math.log(math.log(nd))) + 32))

    def build_sieve(limit):
        spf = list(range(limit + 1))
        primes = []
        spf[0] = 0
        if limit >= 1:
            spf[1] = 1
        for i in range(2, limit + 1):
            if spf[i] == i:
                primes.append(i)
            for p in primes:
                if p * i > limit or p > spf[i]:
                    break
                spf[p * i] = p
        return primes, spf

    limit = estimate_upper(target_index)
    primes, spf = build_sieve(limit)
    while len(primes) < target_index:
        limit *= 2
        primes, spf = build_sieve(limit)

    prime_target = primes[target_index - 1]

    def floor_sum_upto(n, r):
        if r == 0:
            return 0
        r = min(r, n)
        s = 0
        left = 1
        while left <= r:
            q = n // left
            right = min(n // q, r)
            s += q * (right - left + 1)
            left = right + 1
        return s

    def floor_sum_range(n, l, r):
        if l > r or r == 0:
            return 0
        return floor_sum_upto(n, r) - floor_sum_upto(n, l - 1)

    def candidate_delta(n, spf_n):
        max_pd = n // spf_n
        ndc = (n - max_pd - 1) if n > max_pd + 1 else 0
        ndfs = floor_sum_range(n, max_pd + 1, n - 1)
        setup = 2 * n - 1
        ndt = (6 * n + 2) * ndc + 2 * ndfs
        terminal = 5 * n + max_pd + 2 * (n // max_pd) + 2
        return setup + ndt + terminal

    def prime_hit_offset(pc):
        ndfs = floor_sum_range(pc, 2, pc - 1)
        setup = 2 * pc - 1
        ndt = (6 * pc + 2) * (pc - 2) + 2 * ndfs
        pb = 6 * pc + 2
        return setup + ndt + pb

    start_step = 2
    found = 0
    for n in range(2, prime_target + 1):
        spf_n = spf[n]
        is_p = (spf_n == n)
        if is_p:
            hit_step = start_step + prime_hit_offset(n)
            found += 1
            if found == target_index:
                return str(hit_step)
        start_step += candidate_delta(n, spf_n)

    return str(start_step)

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

Java

import java.util.*;

public class Euler308 {
    static long estimateUpperBoundForNthPrime(long n) {
        if (n < 6)
            return 15;
        double nd = (double) n;
        double estimate = nd * (Math.log(nd) + Math.log(Math.log(nd))) + 32.0;
        long bound = (long) Math.ceil(estimate);
        return Math.max(bound, 15);
    }

    static class SieveData {
        int[] primes;
        int primeCount;
        int[] spf;
    }

    static SieveData buildLinearSieve(int limit) {
        SieveData out = new SieveData();
        out.spf = new int[limit + 1];
        out.primes = new int[limit / 10 + 16];
        out.primeCount = 0;

        for (int i = 2; i <= limit; ++i) {
            if (out.spf[i] == 0) {
                out.spf[i] = i;
                if (out.primeCount == out.primes.length) {
                    out.primes = Arrays.copyOf(out.primes, out.primes.length * 2);
                }
                out.primes[out.primeCount++] = i;
            }

            for (int pIndex = 0; pIndex < out.primeCount; ++pIndex) {
                int p = out.primes[pIndex];
                long composite = (long) p * i;
                if (composite > limit)
                    break;
                out.spf[(int) composite] = p;
                if (p == out.spf[i])
                    break;
            }
        }
        return out;
    }

    static SieveData buildSieveWithPrimeCount(long targetIndex) {
        long limit = estimateUpperBoundForNthPrime(targetIndex);
        while (true) {
            if (limit > Integer.MAX_VALUE)
                throw new RuntimeException("Sieve too large");
            SieveData sieve = buildLinearSieve((int) limit);
            if (sieve.primeCount >= targetIndex) {
                return sieve;
            }
            limit *= 2;
        }
    }

    static long floorSumUpto(long n, long r) {
        if (r == 0)
            return 0;
        if (r > n)
            r = n;

        long sum = 0;
        long left = 1;
        while (left <= r) {
            long q = n / left;
            long right = n / q;
            if (right > r)
                right = r;
            sum += q * (right - left + 1);
            left = right + 1;
        }
        return sum;
    }

    static long floorSumRange(long n, long l, long r) {
        if (l > r || r == 0)
            return 0;
        return floorSumUpto(n, r) - floorSumUpto(n, l - 1);
    }

    static long candidateDeltaToNextStart(long n, long spfN) {
        long maxProperDivisor = n / spfN;
        long nonDivisorCount = (n > maxProperDivisor + 1) ? (n - maxProperDivisor - 1) : 0;
        long nonDivisorFloorSum = floorSumRange(n, maxProperDivisor + 1, n - 1);

        long setup = 2 * n - 1;
        long nonDivisorTotal = (6 * n + 2) * nonDivisorCount + 2 * nonDivisorFloorSum;
        long terminal = 5 * n + maxProperDivisor + 2 * (n / maxProperDivisor) + 2;

        return setup + nonDivisorTotal + terminal;
    }

    static long primeHitOffsetFromStart(long primeCandidate) {
        long nonDivisorFloorSum = floorSumRange(primeCandidate, 2, primeCandidate - 1);

        long setup = 2 * primeCandidate - 1;
        long nonDivisorTotal = (6 * primeCandidate + 2) * (primeCandidate - 2) + 2 * nonDivisorFloorSum;
        long primeBranchToHit = 6 * primeCandidate + 2;

        return setup + nonDivisorTotal + primeBranchToHit;
    }

    static long solveIterationsToPrimeIndex(long targetIndex) {
        if (targetIndex == 0)
            return 0;

        SieveData sieve = buildSieveWithPrimeCount(targetIndex);
        long primeTarget = sieve.primes[(int) (targetIndex - 1)];

        long startStep = 2;
        long foundPrimes = 0;

        for (long n = 2; n <= primeTarget; ++n) {
            long spfN = sieve.spf[(int) n];
            boolean isPrime = (spfN == n);

            if (isPrime) {
                long hitStep = startStep + primeHitOffsetFromStart(n);
                foundPrimes++;
                if (foundPrimes == targetIndex) {
                    return hitStep;
                }
            }

            startStep += candidateDeltaToNextStart(n, spfN);
        }

        return 0;
    }

    public static String solve() {
        return String.valueOf(solveIterationsToPrimeIndex(10001));
    }

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