Problem 337: Totient Stairstep Sequences

Project Euler Problem 337 Solution

EulerSolve provides an optimized solution for Project Euler Problem 337, Totient Stairstep Sequences, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary A totient stairstep sequence is a strictly increasing integer sequence $a_1,\dots,a_n$ such that $a_1=6$ and, for every adjacent pair, $$\varphi(a_i) \lt \varphi(a_{i+1}) \lt a_i \lt a_{i+1}.$$ The first inequality forces the totient values to rise, the middle inequality says the next totient must still stay below the previous raw integer, and the last inequality keeps the sequence increasing. We must count all such sequences whose final term is at most $N$, and report the answer modulo $10^8$. The final Project Euler numeric answer is intentionally omitted. Mathematical Approach Formal DP State Let $dp[x]$ be the number of valid sequences whose last term is exactly $x$. The one-term sequence $[6]$ is valid, so $$dp[6]=1.$$ For $x \ge 7$, every valid sequence ending at $x$ must come from a previous endpoint $i \lt x$ satisfying the step condition $$\varphi(i) \lt \varphi(x) \lt i.$$ Therefore $$dp[x]=\sum_{i=6}^{x-1} dp[i]\,\mathbf{1}_{\varphi(i) \lt \varphi(x) \lt i}.$$ The required total is then $$S(N)=\sum_{x=6}^{N} dp[x]\pmod{10^8}.$$ This recurrence already captures the full combinatorics of the problem, but evaluating it literally would compare each $x$ against all smaller $i$, leading to quadratic time. Why a Fixed Predecessor Creates an Interval Fix a predecessor $i$....

Detailed mathematical approach

Problem Summary

A totient stairstep sequence is an integer sequence $a_1,a_2,\dots,a_n$ of some length $n\ge1$ with $a_1=6$ and, for every $i$ with $1\le i\lt n$,

$$\varphi(a_i)\lt\varphi(a_{i+1})\lt a_i\lt a_{i+1},$$

where $\varphi(m)$ counts the integers in $1,\dots,m$ coprime to $m$. ^[2] Let $S(N)$ be the number of such sequences with $a_n\le N$, the one-term $\{6\}$ included. The statement gives $S(10)=4$ — namely $\{6\}$, $\{6,8\}$, $\{6,8,9\}$, $\{6,10\}$ — and the checkpoints $S(100)=482073668$ and $S(10\,000)\bmod 10^8=73808307$. We must find $S(20\,000\,000)\bmod 10^8$. ^[1]

The whole solution in one sentence: walk through the possible last terms in increasing order, counting the sequences that end at each one, which stays cheap because a term $y$ accepts exactly those continuations whose totient falls strictly between $\varphi(y)$ and $y$, so $y$ need only add its own count once to that one interval.

Mathematical Approach

1. Unpacking the step condition

The displayed chain is three separate inequalities. Writing $y=a_i$ for the current term and $x=a_{i+1}$ for the next one, it says

$$\varphi(y)\lt\varphi(x),\qquad \varphi(x)\lt y,\qquad y\lt x.$$

We write $y\to x$ when all three hold. For example $6\to8$, because $\varphi(6)=2$, $\varphi(8)=4$ and $2\lt4\lt6\lt8$. But $6\not\to9$: here $\varphi(9)=6$, and the middle inequality would need $6\lt6$. Yet $6\to8\to9$ is allowed, since $\varphi(8)=4\lt\varphi(9)=6\lt8\lt9$: the bound on $\varphi(9)$ is set by its own predecessor $8$, not by the earlier term $6$.

The condition involves adjacent terms only: nothing links $a_1$ to $a_3$, so a sequence is valid precisely when every consecutive pair is a step. The valid sequences are therefore the directed paths starting at $6$ in the graph on $\{6,\dots,N\}$ with edges $y\to x$. Every edge increases the number, so this graph has no cycles and the plain order $6,7,8,\dots$ lists its vertices with every edge pointing forwards.

2. Counting by the last term

Let $\mathrm{dp}(x)$ be the number of valid sequences whose last term is exactly $x$, and set $\mathrm{dp}(y)=0$ for $y\lt6$. The only sequence ending at $6$ is $\{6\}$, so $\mathrm{dp}(6)=1$. Deleting the last term of a longer sequence ending at $x$ leaves a valid sequence ending at some $y$ with $y\to x$, and by locality every such pair recombines into a valid sequence. This bijection gives, with

$$P(x)=\{\,y:\varphi(x)\lt y\lt x,\ \varphi(y)\lt\varphi(x)\,\}$$

and for $x\ge7$,

$$\mathrm{dp}(x)=\sum_{y\in P(x)}\mathrm{dp}(y),$$

while the quantity asked for is

$$S(N)=\sum_{x=6}^{N}\mathrm{dp}(x).$$

Because every element of $P(x)$ is smaller than $x$, sweeping $x$ upwards keeps every value on the right already known. Evaluated literally, though, the sum inspects every pair $y\lt x$: about $2\cdot10^{14}$ pairs at the target bound, which is hopeless.

3. Which numbers can occur at all

A number $x$ can be reached only if some integer $y$ satisfies $\varphi(x)\lt y\lt x$, and such a $y$ exists exactly when $\varphi(x)\le x-2$. Now $\varphi(x)=x-1$ holds precisely for prime $x$: a prime is coprime to all of $1,\dots,x-1$, ^[2] while a composite $x\ge4$ has a divisor $d$ with $1\lt d\lt x$, and this $d$ is not coprime to $x$, so at most $x-2$ of $1,\dots,x-1$ are.

Hence $\mathrm{dp}(x)=0$ for every prime $x\ge7$. Read for the earlier term, the same inequality shows that a prime has no continuation either, the open interval between $\varphi(y)$ and $y$ being empty. Primes are invisible from both sides, and every term of every sequence is composite. A brute-force check adds that for $6\le x\le100\,000$ the numbers with $\mathrm{dp}(x)=0$ are exactly the $9589$ primes there: every composite in that range is reachable.

4. One predecessor, one interval of totient values

Fix $y$ and ask which $x$ may follow it. The third inequality, $y\lt x$, needs no work: the sweep of Section 2 takes the numbers in increasing order, so every $y$ already entered satisfies it. The other two mention $x$ only through the single value $\varphi(x)$: they say $\varphi(y)\lt\varphi(x)$ and $\varphi(x)\lt y$. Since totients are integers, together they say

$$\varphi(y)+1\le\varphi(x)\le y-1.$$

So among the later numbers the continuations of $y$ are exactly those whose totient lands in one interval depending on $y$ alone, and rather than asking for every $x$ which $y$ are admissible we let every finished $y$ announce its interval once. Precisely: for $x\ge7$ and $1\le t\le N$ define the accumulator

$$A_x(t)=\sum_{6\le y\lt x}\mathrm{dp}(y)\,[\varphi(y)\lt t\lt y],$$

where $[C]$ is $1$ when the condition $C$ holds and $0$ otherwise. Claim: $\mathrm{dp}(x)=A_x(\varphi(x))$. Proof. Both sides sum $\mathrm{dp}(y)$ over a set of $y$, so it suffices that the sets agree. At $t=\varphi(x)$ the accumulator keeps exactly the $y\lt x$ with $\varphi(y)\lt\varphi(x)\lt y$, which are the conditions defining $P(x)$; restricting to $y\ge6$ changes nothing, as $\mathrm{dp}$ vanishes below $6$. $\square$

5. Range add and point query

By Section 4 we need two operations on an array indexed by totient values $1,\dots,N$: add a number to a whole interval, and read one entry. In a plain array the read is free but an interval addition touches every entry in it. The remedy is to store the difference array $D_t=A(t)-A(t-1)$, $A(0)=0$; adding $v$ on $[l,r]$ then changes only two entries,

$$D_l\leftarrow D_l+v,\qquad D_{r+1}\leftarrow D_{r+1}-v,$$

at the price of turning a read into a prefix sum,

$$A(t)=\sum_{u=1}^{t}D_u.$$

which is right because a pair of corrections at $l$ and $r+1$ contributes $v$ to that prefix sum exactly when $l\le t\le r$. Both operations become cheap in a Fenwick tree, or binary indexed tree, which stores sums of nested blocks of $D$. ^[6] Let $b(i)$ be the largest power of two dividing $i$, the bitwise AND of $i$ and $-i$; cell $i$ holds the sum of the $b(i)$ entries of $D$ ending at $i$. So cell $6$ holds $D_5+D_6$, cell $4$ holds $D_1+\dots+D_4$, and those two alone give $A(6)$. Since $i\mapsto i-b(i)$ clears the lowest set bit, a prefix sum reads at most $\log_2 N$ cells; and $i\mapsto i+b(i)$, which visits exactly the blocks containing index $i$, at least doubles $b(i)$, so an update costs the same. Two such updates with one prefix read are the standard range-add, point-query use of the structure. ^[7]

The algorithm is now short. Start the total at $1$ for $\{6\}$ and deposit weight $1$ on $[\varphi(6)+1,5]=[3,5]$. Then for $x=7,8,\dots,N$ in turn, read $\mathrm{dp}(x)$ as the point query at index $\varphi(x)$, add it to the total, and deposit $\mathrm{dp}(x)$ on $[\varphi(x)+1,x-1]$. The query precedes the deposit, so only $y\lt x$ have been entered, as $A_x$ requires; no number could contribute to itself anyway, its query index lying strictly below the interval it deposits on. Indices stay in range because $\varphi(x)\le x-1\le N-1$.

6. Worked example: $S(14)$ by hand

The totients are $\varphi(6)=2$, $\varphi(7)=6$, $\varphi(8)=4$, $\varphi(9)=6$, $\varphi(10)=4$, $\varphi(11)=10$, $\varphi(12)=4$, $\varphi(13)=12$ and $\varphi(14)=6$. ^[3] Running the algorithm:

$x=6$: $\mathrm{dp}(6)=1$; deposits $1$ on $[3,5]$.
$x=7$: query at $6$ finds nothing, $\mathrm{dp}(7)=0$; its interval $[7,6]$ is empty.
$x=8$: query at $4$ hits $[3,5]$ from $6$, so $\mathrm{dp}(8)=1$; deposits $1$ on $[5,7]$.
$x=9$: query at $6$ hits $[5,7]$ from $8$ but not $[3,5]$, so $\mathrm{dp}(9)=1$; deposits $1$ on $[7,8]$.
$x=10$: query at $4$ hits only $[3,5]$, so $\mathrm{dp}(10)=1$; deposits $1$ on $[5,9]$.
$x=11$: query at $10$ finds nothing, $\mathrm{dp}(11)=0$.
$x=12$: query at $4$ hits only $[3,5]$, so $\mathrm{dp}(12)=1$; deposits $1$ on $[5,11]$.
$x=13$: query at $12$ finds nothing, $\mathrm{dp}(13)=0$.
$x=14$: query at $6$ hits $[5,7]$, $[5,9]$ and $[5,11]$, from $8$, $10$ and $12$, and misses $[3,5]$ and $[7,8]$; so $\mathrm{dp}(14)=3$.

Adding up gives $S(10)=1+0+1+1+1=4$, the four sequences of the statement. The three counted by $\mathrm{dp}(14)$ are $\{6,8,14\}$, $\{6,10,14\}$ and $\{6,12,14\}$; with $\{6,12\}$ they take $S(10)=4$ to $S(14)=8$, as direct enumeration confirms. Note how $7$, $11$ and $13$ drop out as Section 3 predicts, and how the deposit of $9$ is never read by any $x\le14$.

7. All the totients at once

The algorithm needs $\varphi(m)$ for every $m\le N$, so they are built in one sweep. Euler's product formula states ^[2]

$$\varphi(m)=m\prod_{p\mid m}\Bigl(1-\frac1p\Bigr),$$

the product running over the distinct primes dividing $m$. A sieve realises it: initialise $\varphi[m]=m$, then for each prime $p$ in increasing order visit the multiples $m=p,2p,3p,\dots$ and replace $\varphi[m]$ by $\varphi[m]-\varphi[m]/p$. ^[4] The cell for $12$ starts at $12$, drops to $6$ at $p=2$ and to $4$ at $p=3$, giving $\varphi(12)=4$.

That integer division is exact, which is worth checking. Writing $q_1,\dots,q_k$ for the primes below $p$ that divide $m$, the cell holds, when $p$ is reached,

$$m\prod_{j=1}^{k}\Bigl(1-\frac1{q_j}\Bigr),$$

which equals

$$\frac{m}{q_1\cdots q_k}\,(q_1-1)\cdots(q_k-1).$$

Since $p$ divides $m$ and differs from each $q_j$, the factor $m/(q_1\cdots q_k)$ is still a multiple of $p$, so the cell is divisible by $p$ and the subtraction multiplies it by exactly $1-1/p$. After the last prime divisor the cell holds $\varphi(m)$.

The primes come from a sieve of Eratosthenes, which may begin crossing out at $p^2$ because smaller multiples of $p$ already carry a smaller prime factor. ^[5] Both sieves cost $O(N\log\log N)$ operations, the bound coming from the sum of $1/p$ over primes $p\le N$. ^[4] ^[5] The C++ version sieves the totients block by block, which adds a surcharge counted below.

8. Working modulo $10^8$

The counts grow ferociously. $S(100)=482073668$ still has nine digits, but the exact $S(10\,000)$ has $373$ decimal digits, and the smallest $N$ with $S(N)\gt2^{64}-1$ is $273$. Exact arithmetic in a fixed-width type is out of the question, which is why a residue is asked for.

Reducing is safe because the recurrence of Section 2 and the accumulator of Section 4 use nothing but addition, and congruence modulo $m$ is compatible with addition: from $a\equiv b$ and $c\equiv d$ follows $a+c\equiv b+d$. ^[8] Keeping every count, every difference-array cell and the total modulo $10^8$ therefore yields the correct residue of $S(N)$. One warning: only the counts are reduced. The comparisons $\varphi(y)\lt\varphi(x)\lt y$ and the interval endpoints use exact totients, and reducing those would change the problem.

The modulus also fixes the integer width: since $2\,(10^8-1)\lt2^{31}-1$, a stored cell plus an incoming delta never leaves the range of a signed or unsigned 32-bit integer, and one conditional subtraction of $10^8$ restores the residue without a division.

How the Code Works

The C++ program runs the totient sieve of Section 7 first and then the counting sweep of Section 5. parse_arguments reads the optional flags --n=, replacing the default bound $20\,000\,000$, and --skip-checks, suppressing the validation. sieve_primes is the sieve of Eratosthenes, crossing out from $p^2$ with 64-bit loop variables so that $p\cdot p$ cannot overflow. The totient routine initialises $\varphi[x]=x$ and applies cur -= cur / p to the multiples of every prime, in blocks of $2^{20}$ entries: it enters each block at the first multiple $\lceil l/p\rceil\,p$ of $p$ inside it and leaves the prime loop once $p$ exceeds the block end. The result equals a single pass, but the cells touched between consecutive primes span megabytes instead of the whole table.

FenwickRangeAddPointQueryMod is Section 5. Its private add_point walks upwards adding a delta, one conditional subtraction of the modulus per cell; point_query walks downwards accumulating the prefix sum; range_add clamps the interval to $[1,N]$, returns at once if it is empty or the delta zero, and otherwise calls add_point(left, delta) and add_point(right + 1, kMod - delta), the modular $-v$ correction.

count_sequences_mod is the sweep: total $1$, the tree seeded by range_add(phi[6] + 1, 5, 1) — the interval $[3,5]$ of the worked example — and then, for each $x$ from $7$ to $N$, one point_query(phi[x]), one addition to the total and one range_add(phi[x] + 1, x - 1, dp_x). count_sequences_exact_bruteforce is the quadratic recurrence of Section 2 written out literally, in exact 64-bit arithmetic without a modulus, and serves only to validate the fast routine on small bounds. validate_totients compares $\varphi(1),\dots,\varphi(10)$ with the stored list $1,1,2,2,4,2,6,4,6,4$, ^[3] and run_validations chains the checkpoints listed below before main prints the requested residue.

The Java program is a direct port: sievePrimes, computeTotients, the nested class Fenwick and countSequencesMod follow their C++ counterparts, except that the totients are built in one unblocked pass and the primes come back as a List<Integer>, so each of the $1\,270\,607$ primes up to the bound becomes a boxed object — tens of megabytes that C++ avoids with a std::vector<int>. Java has no unsigned types, but by Section 8 every sum stays below $2\cdot10^8$, so its signed int arrays are safe.

The Python file is a wrapper rather than a third implementation, the motive being speed: an interpreted mirror of the same loops already needs seconds at $N=10^6$. solve uses the repository layout to find solutionsCpp/Euler337.cpp one directory above its own and, whenever the cached binary .euler337_py_bridge is missing or older than that source, compiles it with clang++ or, failing that, g++ at -std=c++17 -O2. It then runs the binary with no arguments, so the C++ checks are enabled, and takes the result from parse_output: the last line matching ANSWER_RE, else EQUAL_RE, else the last non-empty line.

Complexity and Verification

The sieves cost $O(N\log\log N)$, and the main sweep performs $N-6$ point queries and $N-6$ calls to the range-add routine, each $O(\log N)$, so the running time is

$$O(N\log\log N)+O(N\log N)=O(N\log N).$$

The C++ blocks add a term of their own: in the block ending at $r$ the routine restarts the prime list and inspects every prime up to $r$, including those with no multiple in the block, so for a fixed block size $B$ the surcharge grows like $N^2/(B\log N)$. At the target bound it is $14\,338\,384$ inspections against $61\,128\,033$ totient cell updates — a quarter more work in that stage, still far below $N\log N$.

The sweep itself makes $19\,999\,994$ point queries and as many calls to range_add, of which $1\,270\,605$ return without touching the tree: $1\,270\,604$ for the primes above $6$, whose intervals are empty, and one for $x=5\,935\,626$, the only composite in the range whose count is divisible by $10^8$ — skipping it is harmless, since depositing a zero changes nothing. That leaves $18\,729\,389$ deposits, or $18\,729\,390$ counting the initial one for $6$. The inner loops of point_query and add_point visit $719\,389\,690$ cells in all, about $36$ per value of $x$: one query and two updates, each walking on average about half of the $\log_2(2\cdot10^7)\approx24$ levels.

Memory is $O(N)$. The totient table and the Fenwick array each hold $N+1$ cells of 32 bits, about $80$ MB apiece; the prime sieve needs a byte per number plus the prime list, but in C++ both are released before the Fenwick array is built, so the peak is the two large arrays. On the machine used to check this article the C++ run peaked near $180$ MB resident and took about $1.3$ seconds including its checks, the Java port about $1.6$ seconds inside a $200$ MB heap; such figures depend on the environment, the allocation sizes do not.

Only the C++ program checks itself, and only when --skip-checks is absent: it verifies $\varphi(1),\dots,\varphi(10)$, then $S(10)=4$ and $S(100)=482073668$ from the quadratic recurrence, then agreement of that recurrence with the Fenwick routine at $N=100$, then $S(10\,000)\bmod10^8=73808307$ — all three checkpoints printed in the statement. Python inherits these by running the binary with no arguments; Java fixes the bound in solve and carries no numerical checks. One honest limitation: the quadratic routine accumulates exactly in 64 bits, which by Section 8 overflows from $N=273$ on, so it is sound only on the small bounds where it is used.

The following was verified independently for this article, in arbitrary-precision arithmetic wherever exactness is claimed. The quadratic and the Fenwick recurrence give identical exact counts for every $x\le1500$, and their running totals, reduced modulo $10^8$, match the C++ and Java output for every bound $N\le1500$. Direct enumeration of all sequences reproduces the table of Section 6, with $S(10)=4$ and $S(14)=8$. A second implementation that replaces the difference-array Fenwick tree by a segment tree with lazy addition — a structurally different way to do range update and point query — agrees modulo $10^8$ at $N=100$, $N=10\,000$, $N=10^6$, where the residue is $23116575$, and at the target bound; an independent Python mirror gives that $N=10^6$ residue too. The blocked totient sieve matches a plain single-pass sieve for all $m\le3\,200\,000$ at the production block size $2^{20}$, a range crossing three block boundaries, and for smaller block sizes as well. Finally, for $6\le x\le100\,000$ the numbers with $\mathrm{dp}(x)=0$ are exactly the primes; Section 3 proves that half, while the reachability of every composite in this finite range is a computational observation only.

Footnotes and References

The references below supply the background facts used above; the interval reduction, the accumulator identity, the exactness of the totient sieve and the operation counts are derived here.

Project Euler 337 — Totient Stairstep Sequences. The official statement fixes a_1 = 6, imposes the chain of inequalities on adjacent terms, defines S(N) as the number of such sequences with last term at most N, gives S(10) = 4 with its four sequences and the checkpoints S(100) = 482073668 and S(10000) mod 10^8 = 73808307, and asks for the residue at N = 20000000.
Wikipedia — Euler's totient function. The article defines φ(n) as the count of integers up to n relatively prime to n, states that φ(p) = p − 1 for prime p, and gives Euler's product formula φ(n) = n ∏ (1 − 1/p) over the distinct primes dividing n. Sections 3 and 7 rest on these two facts.
OEIS A000010 — Euler totient function phi(n). The entry lists 1, 1, 2, 2, 4, 2, 6, 4, 6, 4, 10, 4, 12, 6, 8, … from n = 1 on: the totients of the worked example in Section 6 and of the C++ validation routine.
cp-algorithms — Euler's totient function. Its section on the totient for all numbers from 1 to n uses the idea of the sieve of Eratosthenes: start from phi[i] = i, and whenever an entry is still unchanged, so that i is prime, apply phi[j] -= phi[j] / i to every multiple j of i. The page gives the complexity O(n log log n) quoted in Section 7.
Wikipedia — Sieve of Eratosthenes. The article explains why crossing out may start at p², all smaller multiples of p being marked already, and derives O(n log log n) from the prime harmonic series approaching log log n; both points support Section 7. Its segmented variants, which sieve one portion of the range at a time with the same operational complexity, are the spirit of the C++ block structure, whose own surcharge is counted separately.
Wikipedia — Fenwick tree. The article describes a structure that stores an array and performs prefix sums and updates in O(log n) node accesses, built on the least significant set bit lsb(i) = i & −i, query indices moving by i − lsb(i) and update indices by i + lsb(i). It records that the structure was proposed by Boris Ryabko in 1989 and later named after Peter Fenwick's 1994 article on cumulative frequency tables.
cp-algorithms — Fenwick tree. Its section "Range Update and Point Query" is the variant used here: from a tree of zeros, adding x to [l, r] is done by the point updates add(l, x) and add(r + 1, −x), and a single entry is read as an ordinary prefix sum, because beyond r the second update cancels the first.
Wikipedia — Modular arithmetic. The article states that congruence modulo m is compatible with addition, subtraction and multiplication, and gives the rule that a₁ ≡ b₁ and a₂ ≡ b₂ imply a₁ + a₂ ≡ b₁ + b₂, which justifies the reduction in Section 8.

Problem 337 source code

C++

#include <algorithm>
#include <atomic>
#include <cstdint>
#include <iostream>
#include <limits>
#include <string>
#include <thread>
#include <vector>

namespace {

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

constexpr int kDefaultN = 20'000'000;
constexpr u32 kMod = 100'000'000U;

struct Options {
    int n = kDefaultN;
    bool run_checks = true;
    bool allow_multithreading = true;
    unsigned requested_threads = 0U;
};

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

    const std::string tail = arg.substr(p.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) return false;
        parsed = parsed * 10ULL + digit;
    }

    value = parsed;
    return true;
}

bool parse_unsigned_after_prefix(const std::string& arg,
                                 const char* prefix,
                                 unsigned& value) {
    u64 parsed = 0ULL;
    if (!parse_u64_after_prefix(arg, prefix, parsed)) return false;
    if (parsed > static_cast<u64>(std::numeric_limits<unsigned>::max())) return false;
    value = static_cast<unsigned>(parsed);
    return true;
}

bool parse_int_after_prefix(const std::string& arg, const char* prefix, int& value) {
    u64 parsed = 0ULL;
    if (!parse_u64_after_prefix(arg, prefix, parsed)) return false;
    if (parsed > static_cast<u64>(std::numeric_limits<int>::max())) return false;
    value = static_cast<int>(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-checks") {
            options.run_checks = false;
            continue;
        }
        if (arg == "--single-thread") {
            options.allow_multithreading = false;
            continue;
        }

        unsigned parsed_threads = 0U;
        if (parse_unsigned_after_prefix(arg, "--threads=", parsed_threads)) {
            options.requested_threads = parsed_threads;
            continue;
        }

        int parsed_n = 0;
        if (parse_int_after_prefix(arg, "--n=", parsed_n)) {
            options.n = parsed_n;
            continue;
        }

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

    if (options.n < 1) {
        std::cerr << "N must be positive.\n";
        return false;
    }
    return true;
}

unsigned choose_thread_count(const bool allow_multithreading,
                             const unsigned requested_threads,
                             const std::size_t workload_units) {
    if (!allow_multithreading || workload_units < 2ULL) return 1U;

    unsigned threads = requested_threads;
    if (threads == 0U) {
        threads = std::thread::hardware_concurrency();
        if (threads == 0U) threads = 1U;
    }

    return std::max(1U, std::min<unsigned>(threads, static_cast<unsigned>(workload_units)));
}

std::vector<int> sieve_primes(const int limit) {
    std::vector<int> primes;
    if (limit < 2) return primes;

    std::vector<std::uint8_t> is_composite(static_cast<std::size_t>(limit + 1), 0U);
    primes.reserve(static_cast<std::size_t>(limit / 10));

    for (int i = 2; i <= limit; ++i) {
        if (is_composite[static_cast<std::size_t>(i)] == 0U) {
            primes.push_back(i);
            if (static_cast<u64>(i) * static_cast<u64>(i) <= static_cast<u64>(limit)) {
                for (u64 j = static_cast<u64>(i) * static_cast<u64>(i);
                     j <= static_cast<u64>(limit);
                     j += static_cast<u64>(i)) {
                    is_composite[static_cast<std::size_t>(j)] = 1U;
                }
            }
        }
    }

    return primes;
}

std::vector<u32> compute_totients_parallel(const int n, const unsigned thread_count) {
    std::vector<u32> phi(static_cast<std::size_t>(n + 1), 0U);
    if (n >= 1) phi[1] = 1U;

    const std::vector<int> primes = sieve_primes(n);
    if (n < 2) return phi;

    constexpr int kBlockSize = 1 << 20;
    const int block_count = (n + kBlockSize - 1) / kBlockSize;
    std::atomic<int> next_block{0};

    const auto process_block = [&](const int block_id) {
        const int left = block_id * kBlockSize + 1;
        const int right = std::min(n, (block_id + 1) * kBlockSize);

        for (int x = left; x <= right; ++x) {
            phi[static_cast<std::size_t>(x)] = static_cast<u32>(x);
        }

        for (int p : primes) {
            if (p > right) break;
            const int start = ((left + p - 1) / p) * p;
            for (int x = start; x <= right; x += p) {
                u32& cur = phi[static_cast<std::size_t>(x)];
                cur -= cur / static_cast<u32>(p);
            }
        }
    };

    if (thread_count <= 1U || block_count <= 1) {
        for (int block = 0; block < block_count; ++block) {
            process_block(block);
        }
        return phi;
    }

    std::vector<std::thread> workers;
    workers.reserve(thread_count);
    for (unsigned t = 0; t < thread_count; ++t) {
        workers.emplace_back([&]() {
            while (true) {
                const int block_id = next_block.fetch_add(1, std::memory_order_relaxed);
                if (block_id >= block_count) break;
                process_block(block_id);
            }
        });
    }
    for (std::thread& th : workers) th.join();

    return phi;
}

class FenwickRangeAddPointQueryMod {
  public:
    explicit FenwickRangeAddPointQueryMod(const int n)
        : n_(n), bit_(static_cast<std::size_t>(n + 1), 0U) {}

    void range_add(int left, int right, const u32 delta) {
        if (left < 1) left = 1;
        if (right > n_) right = n_;
        if (left > right || delta == 0U) return;

        add_point(left, delta);
        if (right + 1 <= n_) {
            const u32 neg = (delta == 0U) ? 0U : (kMod - delta);
            add_point(right + 1, neg);
        }
    }

    u32 point_query(const int index) const {
        u32 out = 0U;
        int i = index;
        while (i > 0) {
            out += bit_[static_cast<std::size_t>(i)];
            if (out >= kMod) out -= kMod;
            i -= i & -i;
        }
        return out;
    }

  private:
    void add_point(int index, const u32 delta) {
        int i = index;
        while (i <= n_) {
            u32 next = bit_[static_cast<std::size_t>(i)] + delta;
            if (next >= kMod) next -= kMod;
            bit_[static_cast<std::size_t>(i)] = next;
            i += i & -i;
        }
    }

    int n_ = 0;
    std::vector<u32> bit_;
};

u32 count_sequences_mod(const int n, const std::vector<u32>& phi) {
    if (n < 6) return 0U;

    FenwickRangeAddPointQueryMod fenwick(n);
    u32 total = 1U;  // Sequence [6].
    fenwick.range_add(static_cast<int>(phi[6]) + 1, 5, 1U);

    for (int x = 7; x <= n; ++x) {
        const u32 dp_x = fenwick.point_query(static_cast<int>(phi[static_cast<std::size_t>(x)]));
        total += dp_x;
        if (total >= kMod) total -= kMod;
        fenwick.range_add(static_cast<int>(phi[static_cast<std::size_t>(x)]) + 1, x - 1, dp_x);
    }

    return total;
}

u64 count_sequences_exact_bruteforce(const int n, const std::vector<u32>& phi) {
    if (n < 6) return 0ULL;

    std::vector<u64> dp(static_cast<std::size_t>(n + 1), 0ULL);
    dp[6] = 1ULL;
    u64 total = 1ULL;

    for (int x = 7; x <= n; ++x) {
        const int phix = static_cast<int>(phi[static_cast<std::size_t>(x)]);
        u64 sum = 0ULL;
        for (int i = 6; i < x; ++i) {
            if (phi[static_cast<std::size_t>(i)] < static_cast<u32>(phix) && phix < i) {
                sum += dp[static_cast<std::size_t>(i)];
            }
        }
        dp[static_cast<std::size_t>(x)] = sum;
        total += sum;
    }

    return total;
}

bool validate_totients(const std::vector<u32>& phi) {
    if (phi.size() < 11U) return false;
    const std::vector<u32> expected = {0U, 1U, 1U, 2U, 2U, 4U, 2U, 6U, 4U, 6U, 4U};
    for (std::size_t i = 1; i <= 10U; ++i) {
        if (phi[i] != expected[i]) {
            std::cerr << "Totient validation failed at n=" << i << ": got " << phi[i]
                      << ", expected " << expected[i] << '\n';
            return false;
        }
    }
    return true;
}

bool run_validations(const bool allow_multithreading, const unsigned requested_threads) {
    {
        const std::vector<u32> phi_100 = compute_totients_parallel(100, 1U);
        if (!validate_totients(phi_100)) return false;

        const u64 s10 = count_sequences_exact_bruteforce(10, phi_100);
        if (s10 != 4ULL) {
            std::cerr << "Validation failed: S(10) = " << s10 << ", expected 4\n";
            return false;
        }

        const u64 s100 = count_sequences_exact_bruteforce(100, phi_100);
        if (s100 != 482'073'668ULL) {
            std::cerr << "Validation failed: S(100) = " << s100
                      << ", expected 482073668\n";
            return false;
        }

        const u32 s100_mod = count_sequences_mod(100, phi_100);
        if (s100_mod != static_cast<u32>(s100 % static_cast<u64>(kMod))) {
            std::cerr << "Validation failed: modular solver mismatch at N=100\n";
            return false;
        }
    }

    const unsigned threads_for_checks = choose_thread_count(
        allow_multithreading, requested_threads, std::max<std::size_t>(2ULL, 8ULL));

    {
        const std::vector<u32> phi_10k = compute_totients_parallel(10'000, threads_for_checks);
        const u32 s10k = count_sequences_mod(10'000, phi_10k);
        if (s10k != 73'808'307U) {
            std::cerr << "Validation failed: S(10000) mod 10^8 = " << s10k
                      << ", expected 73808307\n";
            return false;
        }
    }

    if (threads_for_checks > 1U) {
        const int test_n = 300'000;
        const std::vector<u32> phi_single = compute_totients_parallel(test_n, 1U);
        const std::vector<u32> phi_multi = compute_totients_parallel(test_n, threads_for_checks);
        if (phi_single != phi_multi) {
            std::cerr << "Validation failed: single-thread and multi-thread totients differ\n";
            return false;
        }

        const u32 s_single = count_sequences_mod(test_n, phi_single);
        const u32 s_multi = count_sequences_mod(test_n, phi_multi);
        if (s_single != s_multi) {
            std::cerr << "Validation failed: thread consistency mismatch in DP stage\n";
            return false;
        }
    }

    return true;
}

}  // namespace

int main(int argc, char** argv) {
    Options options;
    if (!parse_arguments(argc, argv, options)) return 1;

    if (options.run_checks && !run_validations(options.allow_multithreading, options.requested_threads)) {
        return 1;
    }

    const int n = options.n;
    const std::size_t totient_work_units =
        static_cast<std::size_t>((n + (1 << 20) - 1) / (1 << 20));
    const unsigned threads_for_totients = choose_thread_count(
        options.allow_multithreading, options.requested_threads, std::max<std::size_t>(1ULL, totient_work_units));

    const std::vector<u32> phi = compute_totients_parallel(n, threads_for_totients);
    const u32 answer = count_sequences_mod(n, phi);
    std::cout << answer << '\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 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 = subprocess.check_output([str(binary)], text=True)
    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.util.*;

public class Euler337 {

    static final int MOD = 100000000;

    static List<Integer> sievePrimes(int limit) {
        List<Integer> primes = new ArrayList<>();
        if (limit < 2)
            return primes;
        byte[] isComposite = new byte[limit + 1];

        for (int i = 2; i <= limit; i++) {
            if (isComposite[i] == 0) {
                primes.add(i);
                if ((long) i * i <= limit) {
                    for (long j = (long) i * i; j <= limit; j += i) {
                        isComposite[(int) j] = 1;
                    }
                }
            }
        }
        return primes;
    }

    static int[] computeTotients(int n) {
        int[] phi = new int[n + 1];
        for (int i = 0; i <= n; i++)
            phi[i] = i;

        List<Integer> primes = sievePrimes(n);
        for (int p : primes) {
            for (int x = p; x <= n; x += p) {
                phi[x] -= phi[x] / p;
            }
        }
        return phi;
    }

    static class Fenwick {
        int n;
        int[] bit;

        Fenwick(int n) {
            this.n = n;
            bit = new int[n + 1];
        }

        void addPoint(int index, int delta) {
            for (int i = index; i <= n; i += i & -i) {
                int next = bit[i] + delta;
                if (next >= MOD)
                    next -= MOD;
                bit[i] = next;
            }
        }

        void rangeAdd(int left, int right, int delta) {
            if (left < 1)
                left = 1;
            if (right > n)
                right = n;
            if (left > right || delta == 0)
                return;

            addPoint(left, delta);
            if (right + 1 <= n) {
                int neg = (delta == 0) ? 0 : (MOD - delta);
                addPoint(right + 1, neg);
            }
        }

        int pointQuery(int index) {
            int out = 0;
            for (int i = index; i > 0; i -= i & -i) {
                out += bit[i];
                if (out >= MOD)
                    out -= MOD;
            }
            return out;
        }
    }

    static int countSequencesMod(int n, int[] phi) {
        if (n < 6)
            return 0;
        Fenwick fenwick = new Fenwick(n);
        int total = 1;
        fenwick.rangeAdd(phi[6] + 1, 5, 1);

        for (int x = 7; x <= n; x++) {
            int dpX = fenwick.pointQuery(phi[x]);
            total += dpX;
            if (total >= MOD)
                total -= MOD;
            fenwick.rangeAdd(phi[x] + 1, x - 1, dpX);
        }
        return total;
    }

    public static String solve() {
        int n = 20000000;
        int[] phi = computeTotients(n);
        int ans = countSequencesMod(n, phi);
        return String.valueOf(ans);
    }

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