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

A Fractran program is a list of positive fractions. The state of the machine is a positive integer, set at the start to a seed, and one iteration multiplies it by the first fraction of the list that leaves it an integer. We are given John Horton Conway's fourteen-fraction prime generator

$$\begin{aligned}&\tfrac{17}{91},\ \tfrac{78}{85},\ \tfrac{19}{51},\ \tfrac{23}{38},\ \tfrac{29}{33},\\&\tfrac{77}{29},\ \tfrac{95}{23},\ \tfrac{77}{19},\ \tfrac{1}{17},\ \tfrac{11}{13},\\&\tfrac{13}{11},\ \tfrac{15}{2},\ \tfrac{1}{7},\ \tfrac{55}{1}\end{aligned}$$

together with the seed $2$, whose iterations produce 15, 825, 725, 1925, 2275, 425, …, 68, 4, 30, …, 136, 8, 60, …, 544, 32, 240, …. The powers of two that occur are $2^2,2^3,2^5,\dots$: all have prime exponents, and every prime appears, in increasing order. If the program is used to answer Problem 7, finding the $10001$-st prime $p$, how many iterations are needed until its state first equals $2^{p}$? ^[1]

The whole solution in one sentence: the fourteen fractions are a register machine that divides each candidate $n$ by the trial divisors $n-1,n-2,\dots$ through repeated subtraction, so we write the cost of one trial division in closed form, add those costs over the candidates below the target prime, and add the steps until that prime's power of two appears.

Mathematical Approach

1. Fractran as a register machine

A positive integer factors uniquely into primes, so a state is the list of its prime exponents. Only the ten primes $2,3,5,7,11,13,17,19,23,29$ occur here; write $a,b,c,d,e,f,g,h,i,j$ for their exponents, so the state is

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

All fourteen fractions are in lowest terms, and for a reduced $u/v$ the product is an integer exactly when $v$ divides the state, that is, when each register named in $v$ holds at least the required amount; the multiplication then subtracts the exponents of $v$ and adds those of $u$. The state $425=5^2\cdot17$, say, can use $78/85$ and becomes $390=2\cdot3\cdot5\cdot13$: one unit leaves the registers of $5$ and $17$ and one enters those of $2$, $3$ and $13$. A fraction is thus a guarded instruction, and the order of the list is a priority order: the first applicable instruction runs. ^[2] Written as "consumed $\to$ produced" and numbered $R_1,\dots,R_{14}$, the program is

$R_1=17/91$: $7\cdot13\to17$
$R_2=78/85$: $5\cdot17\to2\cdot3\cdot13$
$R_3=19/51$: $3\cdot17\to19$
$R_4=23/38$: $2\cdot19\to23$
$R_5=29/33$: $3\cdot11\to29$
$R_6=77/29$: $29\to7\cdot11$
$R_7=95/23$: $23\to5\cdot19$
$R_8=77/19$: $19\to7\cdot11$
$R_9=1/17$: $17\to1$
$R_{10}=11/13$: $13\to11$
$R_{11}=13/11$: $11\to13$
$R_{12}=15/2$: $2\to3\cdot5$
$R_{13}=1/7$: $7\to1$
$R_{14}=55/1$: $1\to5\cdot11$

No instruction can test a register for zero: that test must be built by letting a guard fail so a later instruction takes over, which is how every loop below ends. ^[2]

2. What the registers mean

Sections 3 to 7 justify this reading step by step. While a candidate $n$ is examined, register $c$ (the exponent of $5$) holds the part of $n$ still to be subtracted and register $a$ (the exponent of $2$) the part already subtracted; register $d$ (the exponent of $7$) holds the trial divisor, which the subtraction eats into and a later phase rebuilds; register $b$ (the exponent of $3$) counts the units transferred in the current round. The exponents of $11,13,17,19,23,29$ are flags saying which phase runs: never more than one of the six is occupied, and never by more than one token. The conserved quantity is $a+c+i=n$, with $i$ the exponent of $23$ parking a unit in transit; in the subtraction loop $i=0$, so there $a+c=n$.

3. The state that starts a candidate

Call

$$C_n=3^{n-1}\cdot5^n\cdot11$$

the start state of candidate $n$, and let $T_n$ be the iteration at which the machine reaches it. From the seed $2$ the first applicable instruction is $R_{12}$ — $R_{14}$ fits every state but has the lowest priority — giving $15=3\cdot5$; on $15$ every guard fails except the last, and $R_{14}$ gives $15\cdot55=825=3\cdot5^2\cdot11=C_2$. So

$$T_2=2,$$

and $825$ is the second term of the sequence displayed in the statement, the third state once the seed is counted.

4. Setting up the first trial divisor

On $C_n$ the flag is $11$. The pair $R_5$ ($3\cdot11\to29$) and $R_6$ ($29\to7\cdot11$) moves one unit from register $b$ to register $d$ and returns the flag to $11$, two iterations each time. It repeats while $b\ge1$, that is $n-1$ times. Then $b=0$, every earlier guard fails, and $R_{11}$ ($11\to13$) switches the flag in one iteration. After

$$2(n-1)+1=2n-1$$

iterations the state is $5^n\cdot7^{n-1}\cdot13$: the candidate in register $c$, the first trial divisor $D=n-1$ in $d$, and the flag $13$.

5. The subtraction loop

Call a state $2^{a}5^{c}7^{D}13$ with $a+c=n$ a pass start. Here $R_1$ ($7\cdot13\to17$) fires, taking one unit off the divisor and setting the flag to $17$; then $R_2$ ($5\cdot17\to2\cdot3\cdot13$) fires if $c\ge1$, moving one unit from $c$ to $a$, copying it into $b$ and restoring the flag. Neither changes $a+c$, as Section 2 promised, and $b$ counts the units subtracted in this pass.

Suppose $c\ge D$. The alternation runs $D$ times each, $2D$ iterations, and stops because $d=0$ fails the guard of $R_1$, leaving $2^{a+D}3^{D}5^{c-D}13$. Now $R_{10}$ ($13\to11$) fires, the $R_5,R_6$ pair moves all $D$ units of $b$ back into $d$ in $2D$ iterations, and $R_{11}$ restores the flag. One pass therefore costs

$$2D+1+2D+1=4D+2$$

iterations and turns $2^{a}5^{c}7^{D}13$ into $2^{a+D}5^{c-D}7^{D}13$: it subtracts $D$ from $c$ and restores the divisor. From $c=n$, after $q=\lfloor n/D\rfloor$ passes the register $c$ holds $r=n\bmod D$.

6. A failed trial division

Suppose $r\ge1$, so $D$ does not divide $n$. The last, partial pass runs the alternation $r$ times each, $2r$ iterations, leaving $c=0$, $b=r$, $d=D-r$. Since $r\le D-1$ we still have $d\ge1$, so $R_1$ fires once more, leaving $d=D-r-1$ and the flag $17$; now $c=0$ blocks $R_2$ and $R_3$ ($3\cdot17\to19$) takes over, spending one unit of $b$ and setting the flag to $19$. That is $2r+2$ iterations.

What follows restores the candidate. The pair $R_4$ ($2\cdot19\to23$) and $R_7$ ($23\to5\cdot19$) moves one unit from $a$ to $c$ and returns the flag; since $a=n$ it repeats until $a=0$, so $2n$ iterations leave $c=n$. Then $R_8$ ($19\to7\cdot11$) adds one to $d$ and sets the flag to $11$; $R_5,R_6$ move the $r-1$ units still in $b$ into $d$ in $2(r-1)$ iterations, and $R_{11}$ restores the flag $13$. The divisor register now holds

$$(D-r-1)+1+(r-1)=D-1,$$

so the machine simply continues with the next smaller trial divisor. Adding the three parts, with $q=\lfloor n/D\rfloor$ and $qD=n-r$,

$$q(4D+2)+(2r+2)+(2n+2r)=6n+2+2q.$$

The remainder cancels: a failed trial division costs exactly $6n+2+2\lfloor n/D\rfloor$ iterations.

7. A successful trial division, and the output

Suppose $D\mid n$. The $q$-th pass empties $c$ with the divisor restored, so the state is $2^{n}7^{D}13$. Then $R_1$ gives $2^{n}7^{D-1}17$; here $b=c=0$ blocks $R_2$ and $R_3$, and $R_4$ to $R_8$ have no flag to use, so $R_9$ ($17\to1$) clears the flag, leaving $2^{n}7^{D-1}$ two iterations later.

If $D=1$ that state is the pure power $2^{n}$ — the machine has printed $n$; if $D\ge2$ the factor $7^{D-1}$ remains and nothing is printed. Either way $R_{12}$ now runs $n$ times (emptying $a$ into $b$ and $c$), $R_{13}$ runs $D-1$ times (clearing $d$) and $R_{14}$ once, producing $3^{n}5^{n+1}11=C_{n+1}$ after $n+D$ further iterations. With $q=n/D$, a successful trial division therefore costs

$$q(4D+2)+2+n+D=5n+D+2(n/D)+2$$

iterations from the pass start to $C_{n+1}$; when $D=1$ the output $2^{n}$ appears $6n+2$ iterations in, $n+1$ before $C_{n+1}$.

8. Why exactly the primes are printed

The trial divisors run downwards from $n-1$, and the scan stops at the first one that divides $n$. That is the largest proper divisor

$$d^{*}=n/\operatorname{spf}(n),$$

where $\operatorname{spf}(n)$ is the smallest prime factor: a divisor $d\lt n$ of $n$ has $n/d\ge\operatorname{spf}(n)$, hence $d\le d^{*}$. For composite $n$, $\operatorname{spf}(n)\le\sqrt n$, so $2\le d^{*}\lt n$ and Section 7 leaves a factor $7^{d^{*}-1}$ with $d^{*}-1\ge1$: nothing is printed. For prime $n$ the scan passes $n-1,\dots,2$ in vain and ends at $D=1$, where $2^{n}$ appears. Apart from the seed, the exponents of the powers of two in the sequence are therefore exactly the primes, in increasing order.

No other power of two can slip through, and the phases alone show it. Sections 4 to 7 account for every state the machine passes while candidate $n$ is examined, and all of them carry one of the flag tokens $11,13,17,19,23,29$, hence an odd prime factor, except in the cleanup that $R_9$ opens. That cleanup starts at $2^{n}7^{d^{*}-1}$, a pure power of two exactly when $d^{*}=1$, and its first $R_{12}$ introduces factors $3$ and $5$ that every later cleanup state keeps until $R_{14}$ builds $C_{n+1}$. So the only power of two candidate $n$ can produce is $2^{n}$, and only if $n$ is prime. The scan always ends, since the divisors decrease and $D=1$ divides $n$; from $C_2$ the machine therefore works through the candidates in order, printing exactly the primes as exponents.

9. The two closed-form costs

Write

$$G(n,l)=\sum_{d=l}^{n-1}\bigl(6n+2+2\lfloor n/d\rfloor\bigr)$$

for the total cost of the failed trial divisions from $n-1$ down to $l$, and

$$E(n)=5n+d^{*}+2\lfloor n/d^{*}\rfloor+2$$

for the successful one. Sections 4, 6 and 7 then give the step from one candidate to the next,

$$\Delta(n)=(2n-1)+G(n,d^{*}+1)+E(n),$$

and, for a prime $p$, the offset from its start state to its output,

$$H(p)=(2p-1)+G(p,2)+(6p+2).$$

For a prime $d^{*}=1$, so $E(p)-(6p+2)=p+1$ as Section 7 said. With $T_2=2$ and

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

the quantity asked for is $T_{p}+H(p)$ for $p=p_{10001}=104743$.

10. A worked example

Take $n=5$, with $T_5=136$. The setup costs $2\cdot5-1=9$, so the first trial division, at $D=4$, begins at iteration $145$. A failed one costs $6\cdot5+2+2\lfloor5/D\rfloor$: that is $34$ for $D=4$ and $D=3$ and $36$ for $D=2$, reaching iterations $179$, $213$ and $249$. At $D=1$ the output comes $6\cdot5+2=32$ iterations later, so

$$H(5)=9+34+34+36+32=145$$

and $2^5=32$ is printed at iteration $136+145=281$ — the $32$ in the fragment "$544$, $32$, $240$" of the statement. Continuing to the next start state replaces the last $32$ by $E(5)=25+1+10+2=38$, so $\Delta(5)=151$ and $T_6=287$.

Take $n=6$ next, with $\operatorname{spf}(6)=2$ and $d^{*}=3$. The setup costs $11$, reaching $D=5$ at $298$; the failures at $D=5$ and $D=4$ cost $6\cdot6+2+2\cdot1=40$ each, reaching $338$ and $378$; the success at $D=3$ costs $E(6)=30+3+4+2=39$. So $\Delta(6)=130$ and $T_7=417$, with no power of two, because the exit left the factor $7^{2}$.

11. Summing the quotients in blocks

The constant part of $G(n,l)$ is one multiplication; the work is in its quotient sum. For $n=12$ the quotients $\lfloor12/d\rfloor$ at $d=1,\dots,12$ are $12,6,4,3,2,2,1,1,1,1,1,1$, so the last two groups can be added as $2\cdot2$ and $1\cdot6$ rather than one denominator at a time. The map $d\mapsto\lfloor n/d\rfloor$ takes $O(\sqrt n)$ distinct values, since $d\le\sqrt n$ gives at most $\sqrt n$ of them and $d\gt\sqrt n$ forces $\lfloor n/d\rfloor\lt\sqrt n$; the divisors sharing one quotient $q$ form the block

$$\lfloor n/(q+1)\rfloor+1\le d\le\lfloor n/q\rfloor.$$

Walking the blocks evaluates $\sum_{d\le r}\lfloor n/d\rfloor$ in $O(\sqrt n)$ time, the standard device behind the $O(\sqrt x)$ evaluation of the divisor summatory function. ^[5] For $n=104743$ the full range uses $646$ blocks, just under $2\sqrt n\approx647$.

How the Code Works

The C++ program ^[9] parses --target=K for the $K$-th prime, rejects a zero target, and runs run_checkpoints before the requested calculation unless --skip-checkpoints was given. estimate_upper_bound_for_nth_prime returns $\lceil n(\ln n+\ln\ln n)\rceil+32$ (and $15$ below $n=6$), the classical shape of a bound on the $n$-th prime, ^[7] while build_sieve_with_prime_count doubles that limit until enough primes are found, so nothing depends on the estimate really being a bound. build_linear_sieve fills spf[x] with the smallest prime factor of every $x$ and collects the primes in order, writing each composite once; ^[6] that table gives both the target prime and the $\operatorname{spf}(n)$ of Section 8.

floor_sum_upto(n, r) implements Section 11: after capping $r$ at $n$ it starts at left $=1$, takes $q=\lfloor n/\mathtt{left}\rfloor$, sets right to $\min(\lfloor n/q\rfloor,r)$, adds $q\,(\mathtt{right}-\mathtt{left}+1)$ and continues at right $+1$; floor_sum_range is the difference of two such prefix sums, and $0$ when $l\gt r$. Then candidate_delta_to_next_start(n, spf_n) is $\Delta(n)$, built from setup $=2n-1$, non_divisor_total $=(6n+2)$ times the number $n-d^{*}-1$ of failed divisors plus twice their quotient sum, and terminal $=E(n)$; prime_hit_offset_from_start is $H(p)$, identical but for the last term $6p+2$.

solve_iterations_to_prime_index runs the recurrence of Section 9. It keeps start_step $=T_n$, initialised to $2$, and loops $n=2,\dots,p$: when spf[n] equals $n$ the candidate is prime, so start_step $+\;H(n)$ is its output iteration and a counter of primes found decides whether this is the one wanted; then start_step advances by $\Delta(n)$. Nothing is kept per candidate, so the loop needs no memory beyond the sieve.

The rest of the file is the self-check. apply_one_fractran_step is the fourteen guards of Section 1, in list order, acting on a FractranState that holds the ten exponents $a,\dots,j$ as named fields, so the simulation is the definition rather than a model of it. is_prime_small admits prime exponents only, keeping the seed $2=2^1$ out but also hiding an unexpected composite exponent, so the claim that every power of two has a prime exponent needs the unfiltered simulation below. brute_prime_hits simulates until it has the requested number of outputs, capped at $2\cdot10^9$ iterations, and run_checkpoints compares its first ten against $(2,19)$, $(3,69)$, $(5,281)$, $(7,710)$, $(11,2375)$, $(13,3893)$, $(17,8102)$, $(19,11361)$, $(23,19268)$, $(29,36981)$ and its first thirty against the closed form at the targets $1,2,3,10,20,30$.

The three languages differ only in arithmetic and packaging. C++ accumulates in __uint128_t, a GCC and Clang extension ^[8] that this target does not need: the bound $\Delta(n)\lt7n^2$ proved below keeps the whole count below $2+7p^3$, under $2^{63}$ for $p\lt110000$, so Java's long suffices. Python counts in exact integers of unlimited width, though its estimate still uses floating-point logarithms. The ports rename the pieces — candidate_delta and prime_hit_offset inside Python's solve, candidateDeltaToNextStart and primeHitOffsetFromStart in Java — and neither carries the simulator or the checkpoints: both compute the closed form only, at the fixed target $10001$.

Complexity and Verification

The sieve costs $O(B)$ time and memory for its limit $B$. Here the estimate gives $B=114352$, a range holding $10820$ primes, so no doubling occurs; in C++ and Java its table is $114353$ four-byte entries, $457412$ bytes or about $447$ KiB, and dominates the memory use, while Python holds object references and pays more per entry. The main loop visits every $n\le p=104743$ once and spends $O(\sqrt n)$ on the quotient blocks, so the computation runs in

$$O(B)+O\bigl(p^{3/2}\bigr)$$

time, a few times $10^{7}$ elementary operations, each quotient block costing a constant number of them. All three programs finish in seconds, the C++ time being dominated by its literal checkpoints rather than by the count itself.

Simulating the program literally is hopeless by comparison. Since $d^{*}\le n/2$, at least $n/2-1$ trial divisions fail, each costing more than $6n$, so $\Delta(n)\gt3n^2-6n$. The other way round, $d^{*}\ge1$ leaves at most $n-2$ failures and $E(n)\le7n+3$, so

$$\begin{aligned}\Delta(n)&\le6n^2-n-2\\&\quad+2\sum_{d=2}^{n-1}\lfloor n/d\rfloor;\end{aligned}$$

and $\sum_{d=2}^{n-1}1/d\le\ln n$ turns this into $\Delta(n)\le6n^2+2n\ln n$, hence $\Delta(n)\lt7n^2$. So $\Delta(n)$ has order $n^2$ and the literal count order $p^3$; for the target prime it exceeds $10^{15}$, which the closed form reaches in about $10^{5}$ turns of its loop. The growth shows in the checkpoints: the $20$-th prime is printed at iteration $508284$ and the $30$-th at $2021938$, but the $100$-th at $213945763$ and the $1000$-th at $665728089948$.

The following was checked by running code.

An independent simulator of the fourteen fractions reproduces the statement's terms $2,15,825,725,1925,2275,425$ and its fragments $68,4,30$, $136,8,60$ and $544,32,240$ at the iterations $18$–$20$, $68$–$70$ and $280$–$282$. Its first ten outputs fall at $19,69,281,710,2375$, $3893,8102$, $11361,19268,36981$, matching the C++ checkpoints and the published sequences. ^[3] ^[4]
Apart from the seed $2=2^1$, the simulator's pure powers of two in the first $400000$ iterations have the exponents $2,3,5,7,11,13,17$, $19,23,29,31,37$, $41,43,47,53,59,61$ — exactly the primes up to $61$, in order.
The start states $C_n$ occur where $T_2=2$ and $T_{n+1}=T_n+\Delta(n)$ predict for every $n\le66$ — the first values are $2,22,73,136,287,417,718$ — and $T_p+H(p)$ is the output iteration of every prime $p\le61$. The closed form matches the simulator at the targets $1,2,3,10,20,30$ too, the last being the prime $113$ printed at $2021938$.
Within the first $200000$ iterations the simulator completed $52$ setup phases, $865$ failed trial divisions and $51$ successful ones, every cost matching Sections 4, 6 and 7 exactly; each failure continued with $D-1$, each success fell at $D=d^{*}$, a full pass cost $4D+2$. Of those iterations the $65109$ whose next instruction is $R_1$ or $R_2$ satisfied $a+c=n$, every one satisfied $a+c+i=n$, and no state held more than one flag token, the $1841$ cleanup states holding none.
The C++, Python and Java programs agree at the prime indices $1,2,3,10,20,30,100,1000$ and at the target $10001$, as does a fourth, independently written evaluation of $\Delta$ and $H$; the C++ checkpoints pass.

Footnotes and References

The references below supply the background facts used above. The register reading of the fourteen fractions, the four phase costs, the proof that only primes are printed and the recurrence for $T_n$ are derived explicitly in this article.

Project Euler 308 — An Amazing Prime-generating Automaton. The official statement defines the Fractran machine, lists Conway's fourteen fractions, shows the sequence the seed 2 produces with its powers of two 2², 2³, 2⁵, and asks for the number of iterations needed to produce 2 raised to the 10001st prime.
Wikipedia — FRACTRAN. The article attributes the language to John Conway and explains that each variable is held as the exponent of a prime, so a fraction tests the variables in its denominator and, if the test passes, decrements them and increments those of its numerator; it notes that a fraction cannot test for zero. It lists the same fourteen fractions as PRIMEGAME and calls it a quotient-and-remainder algorithm inside two loops, dividing by successively smaller numbers to find the largest divisor.
OEIS A007542 — successive integers produced by Conway's PRIMEGAME. The entry lists 2, 15, 825, 725, 1925, 2275, 425, 390, …, states that the terms 2 raised to a prime appear in increasing order, and records the positions 0, 19, 69, 281, 710, 2375, … of the powers of two, used as checkpoints above. It cites Conway's 1987 FRACTRAN paper and R. K. Guy, "Conway's prime producing machine", Mathematics Magazine 56 (1983).
OEIS A007547 — steps to compute the n-th prime in PRIMEGAME. The entry gives 19, 69, 281, 710, 2375, 3893, 8102, 11361, 19268, 36981, … as the iteration at which the n-th prime is printed, the quantity this problem asks for at n = 10001.
Wikipedia — Divisor summatory function. The article writes D(x) as the sum of ⌊x/k⌋ for k ≤ x and evaluates it as 2·Σ_{k≤u} ⌊x/k⌋ − u² with u = ⌊√x⌋, in O(√x) time — a different square-root-time route to the same floor sum, where Section 11 walks the runs of denominators sharing one quotient.
cp-algorithms — Linear sieve. The page presents a sieve that computes the minimum prime factor of every integer up to n in O(n) time, storing it in an array beside the list of primes found so far — the spf table used here.
Wikipedia — Prime number theorem. Its section on the n-th prime gives Dusart's bounds pₙ > n(log n + log log n − 1) for n ≥ 2 and pₙ < n(log n + log log n − 0.9484) for n ≥ 39017 — the shape of the code's estimate, though that threshold belongs to the sharper inequality, not to the looser expression the code uses.
GCC manual — 128-bit integers. The page documents __int128 and unsigned __int128 as extensions available on targets with an integer mode wide enough for 128 bits, the type of the C++ accumulators.
C++ — Euler308.cpp. The linked immutable revision contains floor_sum_upto, candidate_delta_to_next_start, prime_hit_offset_from_start, apply_one_fractran_step and run_checkpoints. It is the implementation record for this article.

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());
    }
}