Problem 411: Uphill Paths

Project Euler Problem 411 Solution

EulerSolve provides an optimized solution for Project Euler Problem 411, Uphill Paths, with C++, Python, Java, and a step-by-step mathematical explanation.

Problem Summary For a given $n$, define the stations $$P_i=\left(2^i \bmod n,\ 3^i \bmod n\right),\qquad 0 \le i \le 2n.$$ An uphill path from $(0,0)$ to $(n,n)$ uses only right and up moves, so the stations visited by such a path must appear in an order whose coordinates never decrease. Let $S(n)$ be the maximum number of stations that can be visited by one uphill path. The program evaluates $$\sum_{k=1}^{30} S(k^5).$$ Mathematical Approach The implementation does not search over grid paths directly. Instead, it reduces the problem to a finite set of modular points and then solves an order problem on those points. 1) Reduce the Station Sequence to a Finite Prefix Write $$n = 2^{\alpha} 3^{\beta} m,\qquad \gcd(m,6)=1.$$ Consider the $x$-coordinate $x_i = 2^i \bmod n$. For every $i \ge \alpha$, the factor $2^{\alpha}$ already divides $2^i$, so modulo $2^{\alpha}$ the value is fixed at $0$. On the coprime part $n / 2^{\alpha}$, the base $2$ is invertible, hence the remaining behavior is purely periodic with period $$T_2 = \operatorname{ord}_{n / 2^{\alpha}}(2),$$ with the harmless convention $T_2=1$ when $n / 2^{\alpha}=1$. Thus the $x$-coordinate is eventually periodic after $\alpha$ steps. Exactly the same argument applies to the $y$-coordinate $y_i = 3^i \bmod n$....

Detailed mathematical approach

Problem Summary

Let $n$ be a positive integer. A station stands at every coordinate pair

$$\begin{aligned}(x,y)&=\left(2^i\bmod n,\ 3^i\bmod n\right),\\ 0&\le i\le 2n,\end{aligned}$$

and equal coordinates count as one station. We want a path from $(0,0)$ to $(n,n)$ along which neither coordinate ever decreases; $S(n)$ is the largest number of stations it can pass through. The statement works out $n=22$: $11$ distinct stations, at most $5$ on one path, so $S(22)=5$. It also records $S(123)=14$ and $S(10000)=48$, and asks for $\sum_{k=1}^{30}S(k^5)$. ^[1]

The whole solution in one sentence: the pair $(2^i\bmod n,\,3^i\bmod n)$ begins to repeat after a short run-in, so the entire station set is a list we can write down, and once it is sorted by $x$ the longest uphill path is exactly the longest non-decreasing subsequence of the $y$ coordinates.

Mathematical Approach

1. From an uphill path to a chain of stations

Compare points by both coordinates: write $P\preceq Q$ when $x_P\le x_Q$ and $y_P\le y_Q$, and call a set a chain if any two members are comparable. A path that never decreases and meets $P$ before $Q$ has $x_P\le x_Q$ and $y_P\le y_Q$, so its stations form a chain. Conversely all coordinates are residues modulo $n$, hence in $\{0,\dots,n-1\}$, so for a chain $P_1\preceq\cdots\preceq P_t$ the staircase $(0,0)\to P_1\to\cdots\to P_t\to(n,n)$ only moves right and up: a legal path through all $t$ stations.

So $S(n)$ is the length of a longest chain of distinct stations: the height of this ordered set, which by Mirsky's theorem also equals the fewest antichains — sets of pairwise incomparable stations — that partition it. ^[2]

2. From a chain to a longest non-decreasing subsequence

Let $Q_1,\dots,Q_r$ be the distinct stations sorted lexicographically — by $x$, and among equal $x$ by $y$ — with $Q_j=(x_j,y_j)$. We claim

$$S(n)=\operatorname{LNDS}(y_1,y_2,\dots,y_r),$$

the length of a longest non-decreasing subsequence of the $y$ list.

If $j_1\lt\cdots\lt j_t$ satisfy $y_{j_1}\le\cdots\le y_{j_t}$, then $x_{j_1}\le\cdots\le x_{j_t}$ holds automatically, the list being $x$-sorted, so those points form a chain of equal length. Conversely $P\preceq Q$ with $P\ne Q$ forces $P$ first lexicographically: either $x_P\lt x_Q$, or $x_P=x_Q$ and then $y_P\lt y_Q$, which is how ties are broken. A chain therefore occupies increasing positions with non-decreasing $y$, and the two maxima agree.

The word non-decreasing matters: a vertical move joins stations with equal $x$, a horizontal one stations with equal $y$. Section 9 shows what a strict comparison would cost.

3. Worked example: $n=22$

Doubling and tripling modulo $22$ from $(1,1)$ gives, for $i=0,\dots,10$, the points $(1,1)$, $(2,3)$, $(4,9)$, $(8,5)$, $(16,15)$, $(10,1)$, $(20,3)$, $(18,9)$, $(14,5)$, $(6,15)$, $(12,1)$. The $x$ values $1,2,4,8,16,10,20,18,14,6,12$ all differ, so these are $11$ distinct stations, the count the statement reports; by Section 5 no later index adds anything.

Sorted by $x$ the order becomes $x=1,2,4,6,8,10,12,14,16,18,20$, with matching $y$ values

$$1,\ 3,\ 9,\ 15,\ 5,\ 1,\ 1,\ 5,\ 15,\ 9,\ 3.$$

A longest non-decreasing subsequence here has length $5$; one choice is $1,1,1,5,9$, the stations $(1,1)$, $(10,1)$, $(12,1)$, $(14,5)$, $(18,9)$. A path running right along $y=1$ and then stepping up visits these five, so $S(22)=5$.

4. The station sequence is eventually periodic

Let $\alpha$ and $\beta$ be the exponents of $2$ and $3$ in $n$, with $m_2=n/2^{\alpha}$ odd and $m_3=n/3^{\beta}$ not divisible by $3$. As $2^{\alpha}$ and $m_2$ are coprime with product $n$, the Chinese remainder theorem pins a residue modulo $n$ down by its residues modulo $2^{\alpha}$ and modulo $m_2$. ^[5]

Treat the parts separately. For $i\ge\alpha$ we have $2^{\alpha}\mid2^i$, so that part is $0$ forever after; modulo $m_2$ the base $2$ is invertible, so its powers cycle with period $T_2=\operatorname{ord}_{m_2}(2)$, the least positive exponent with $2^{T_2}\equiv1$, ^[6] and we set $T_2=1$ when $m_2=1$. Reassembling, $2^{\,i+T_2}\equiv2^i\pmod n$ for $i\ge\alpha$, and likewise for $3$ with $\beta$ and $T_3=\operatorname{ord}_{m_3}(3)$, again $1$ when $m_3=1$. So with

$$\begin{aligned}\ell&=\max(\alpha,\beta),\\ T&=\operatorname{lcm}(T_2,T_3),\end{aligned}$$

the station $P_i=(2^i\bmod n,\,3^i\bmod n)$ satisfies $P_{i+T}=P_i$ for all $i\ge\ell$: a run-in of length $\ell$, then period $T$.

5. The first $\ell+T$ stations are all of them, each exactly once

Nothing new appears later: if $i\ge\ell+T$ then $P_i=P_{i-T}$ with $i-T\ge\ell$, and repeating pulls the index into the window $\ell\le i\lt\ell+T$, which with $0,\dots,\ell-1$ covers every station.

More than that — sharper than the code assumes — the map $i\mapsto P_i$ is injective on $0\le i\lt\ell+T$. Suppose $P_i=P_j$ with $i\lt j\lt\ell+T$. As $2^{\alpha}\mid n$, equal $x$ gives $2^i\equiv2^j\pmod{2^{\alpha}}$. If $i\lt\alpha$ then $2^i$ is a nonzero residue modulo $2^{\alpha}$, while $2^j$ is either $0$ (if $j\ge\alpha$) or a different value below $2^{\alpha}$ — impossible. So $i\ge\alpha$, and the same argument on $y$ with $3^{\beta}$ gives $i\ge\beta$, hence $i\ge\ell$. Now $2$ is invertible modulo $m_2$, so $2^{\,j-i}\equiv1$ and $T_2\mid j-i$; likewise $T_3\mid j-i$, so $T\mid j-i$. But $0\lt j-i\le T-1$, a contradiction.

So the number of distinct stations is exactly $r=\min(2n+1,\ \ell+T)$, and Section 6 shows the minimum is always $\ell+T$. For $n=22$, $m_2=11$ and $m_3=22$ give $T_2=10$, $T_3=5$, $\ell=1$, $T=10$, $r=11$ — the statement's count; for $n=123$, $\ell=1$, $T_2=20$, $T_3=8$, $r=41$; for $n=10000$, $\ell=4$, $T_2=T_3=500$, $r=504$.

6. Why the index bound $2n$ never bites

For $n=1$ the only station is $(0,0)$ and $S(1)=1$; assume $n\ge2$ from here on. Write $\varphi(m)$ for Euler's totient, the count of units modulo $m$. Let $\lambda$ be the Carmichael function, the exponent of that group: $a^{\lambda(m)}\equiv1\pmod m$ for every $a$ coprime to $m$, so every order divides $\lambda$, and $\lambda(m)\mid\lambda(n)$ when $m\mid n$. ^[7] Since $m_2\mid n$ and $m_3\mid n$, both $T_2$ and $T_3$ divide $\lambda(n)$, hence so does $T$. Thus $T\le\lambda(n)\le\varphi(n)\le n-1$, while $2^{\alpha}\le n$ and $3^{\beta}\le n$ give $\ell\le\log_2 n$, so

$$\ell+T\le n-1+\log_2 n\lt2n+1 .$$

The truncation at $i=2n$ never discards a station, and $r=\ell+T$ exactly. It also names the real cost, since $T$ can reach $\varphi(n)$: for $n=29^5=20511149$ both $2$ and $3$ have order $\varphi(n)=19803868$, making them primitive roots, so that one $k$ yields almost twenty million stations — the largest case here.

7. Computing the two multiplicative orders

The exponents $t$ with $a^t\equiv1\pmod m$ are exactly the multiples of $\operatorname{ord}_m(a)$, which divides $\varphi(m)$. ^[6] Hence this algorithm: start with $e=\varphi(m)$ and, for each prime $q\mid e$, replace $e$ by $e/q$ while $q\mid e$ and $a^{e/q}\equiv1\pmod m$.

The invariant $\operatorname{ord}_m(a)\mid e$ holds initially and survives each step, since a replacement needs $a^{e/q}\equiv1$, that is $\operatorname{ord}_m(a)\mid e/q$. At the end no prime $q\mid e$ has $a^{e/q}\equiv1$, so $\operatorname{ord}_m(a)\nmid e/q$ for every such $q$. A prime whose test once failed stays failed, since every later candidate $e/q$ divides the earlier one. Were $e$ larger than the order, $e/\operatorname{ord}_m(a)$ would have a prime factor $q$ with $\operatorname{ord}_m(a)\mid e/q$ — excluded. So $e$ is the order.

8. The smallest-prime-factor table

Both $\varphi(m)$ and the prime factors of $e$ come from one table. A linear sieve finds the smallest prime factor of every integer up to $N$ in $O(N)$ time: walking $i$ upwards it marks $i\cdot q$ for each prime $q$ up to the smallest prime factor of $i$, stopping there so that every composite is marked once. ^[8] Dividing that factor out repeatedly lists the distinct primes of any $m\le N$, giving $\varphi(m)=m\prod_{q\mid m}(1-1/q)$. Here $N=30^5=24300000$, and every number factored is at most $n\le N$.

9. Patience sorting, and why the comparison must be non-strict

A longest non-decreasing subsequence of $r$ numbers is found in $O(r\log r)$ time with one array and one binary search per element. ^[3] Keep $t_0,\dots,t_{L-1}$, where $L$ is the longest length reached so far and $t_j$ is the smallest value ending such a subsequence of length $j+1$. The list is non-decreasing: dropping the last element of a subsequence ending at $t_j$ leaves one of length $j$ ending no later, so $t_{j-1}\le t_j$, and binary search applies.

For a new $v$, find the first index $j$ with $t_j\gt v$. If there is none, append $v$, which extends the subsequence ending at $t_{L-1}\le v$. Otherwise set $t_j=v$: length $j+1$ can now end at $v$, since $j=0$ or $t_{j-1}\le v$, and nothing smaller can, the old minimum $t_j$ having exceeded $v$. At the end $L$ is the answer, and the piles of the equivalent card game are the antichains of Section 1. ^[4]

Searching for the first $t_j\gt v$ is what makes the result non-decreasing; the first $t_j\ge v$ gives the strictly increasing variant, and on these lists the gap is wide: it returns $4$, $9$ and $6$ for $n=22$, $123$ and $243$ in place of $5$, $14$ and $162$.

10. Two more examples: $n=243$ and $n=32$

For $n=243=3^5$: $\beta=5$ makes $y_i=0$ for every $i\ge5$, and $T=\operatorname{ord}_{243}(2)=162$, so $162$ of the $r=167$ stations lie on the line $y=0$. Walking that line visits all $162$, and no longer non-decreasing subsequence exists, so $S(243)=162$. The mirror image is $n=32=2^5$, where $\alpha=5$ forces $x_i=0$ from $i=5$ on: eight of the $r=13$ stations stand on $x=0$, and climbing that line before stepping to $(8,27)$ gives $S(32)=9$, again with nothing longer available.

Both are pure prime powers, and that is what gives the free axis chain: for a general $n$ Section 4 makes a coordinate vanish only modulo $2^{\alpha}$ or $3^{\beta}$, not modulo $n$. For $n=160=2^5\cdot5$ the $x$ coordinates settle into the cycle $32,64,128,96$ and never reach $0$; and when $\gcd(n,6)=1$ both bases are units, so neither coordinate is ever $0$ and no axis chain exists.

How the Code Works

The C++ program ^[10] follows Sections 4 to 9. parse_arguments reads --skip-checkpoints and --k-max=; build_spf is the linear sieve of Section 8, always sized for $30^5$; euler_phi and unique_prime_factors read it; mod_pow_u64 uses __uint128_t products to avoid overflow; and multiplicative_order is Section 7.

S_value does the rest. It returns $1$ at once when $n=1$; otherwise it peels the factors of $2$ and $3$ off $n$ into pre2, mod2, pre3, mod3, takes the two orders as per2 and per3, and forms pre $=\ell$, period $=T$ with std::lcm and needed $=\min(2n+1,\ell+T)$. It then walks the points, doubling and tripling, packing each into one 64-bit key as (x << 32) | y. That packing is the trick behind the sort: both coordinates are below $2^{32}$, so the halves never interfere and the top bit stays clear, and sorting the keys as numbers is the lexicographic sort of Section 2.

radix_sort_u64 does it in four least-significant-first passes of $16$ bits, counting into a $65536$-entry array and permuting; each pass is stable, so the four sort all $64$ bits in $O(r)$ time. ^[9] std::unique then drops duplicates — by Section 5 there are none, so it is a safety net. The final loop extracts $y$ and applies the patience step with std::upper_bound, the search for the first entry strictly greater than $y$; tails then has $S(n)$ entries. solve_sum adds $S(k^5)$ over $k$, and run_checkpoints verifies the three statement values first.

The Java version mirrors these functions, with these real differences. It sorts a long[] with Arrays.sort rather than by radix — correct only because every key is non-negative, so signed and unsigned order agree. Its modPow multiplies plain long values, safe since the modulus stays below $2^{25}$ and products below $2^{50}$. And buildSPF boxes its one and a half million primes into a List<Integer> where C++ keeps a flat std::vector<int>. Java runs no checkpoints.

The Python file is a bridge, not a second implementation: solve takes the problem number from its own filename, compiles the C++ source with clang++ or g++ at -O2 into a cached binary when that is missing or stale, runs it, and parse_output returns the text after the last answer: match, failing that the text after the last =, and otherwise the last non-empty line. One detail deserves stating: should_skip_cpp_checkpoints returns true when --skip-checkpoints occurs anywhere in the C++ source — which it does, since the C++ parses that flag — so the Python route always skips the three checkpoints a direct run performs.

Complexity and Verification

One evaluation of $S(n)$ finds the two orders with $O(\log n)$ modular powers, so $O((\log n)^2)$ word operations, generates $r=\ell+T$ points in $O(r)$, sorts them — $O(r)$ for the four radix passes in C++, $O(r\log r)$ for Java's comparison sort — and runs the patience loop in $O(r\log r)$; by Section 6 $r\le n-1+\log_2 n$, so both are $O(n\log n)$ at worst. The sieve costs $O(N)$ with $N=30^5=24300000$, and the summation loop sorts $49935580$ keys, about twice $N$ — one short of the $49935581$ stations, since $n=1$ returns before any key is made.

Memory binds before time does. At the largest case, $k=29$, the key buffer holds $19803868$ eight-byte keys, some $151$ MiB, the radix scratch buffer the same again, and the factor table $93$ MiB. The scratch buffer is freed before tails is built; tails holds $S(n)$ four-byte entries, all $4r$ of which Java reserves in advance. On one machine the peak resident memory was about $805$ MiB for C++ and $1047$ MiB for Java, the times $2.4$ seconds for C++ with GCC 16 at -O2, $1.5$ for the bridge's Clang build and $4.9$ for Java.

Checks run for this article, outside the published code:

$S(22)=5$, $S(123)=14$ and $S(10000)=48$, reproduced from the literal definition.
A quadratic longest-chain dynamic programme on the literal station set matches the $\operatorname{LNDS}$ value for all $n\le300$, and the prefix of length $\ell+T$ gives the same $S(n)$ as the full range $0\le i\le2n$ for all $n\le1500$ (Sections 1, 2 and 4).
That prefix repeats no station for all $n$ from $2$ to $6000$ and all thirty moduli $k^5$, and $T\mid\lambda(n)$ with $\ell+T\le2n+1$ for all $n\le40000$ (Sections 5 and 6).
The patience array stays non-decreasing and its length matches a quadratic reference on $5000$ random lists (Section 9).
The C++ and Java programs and a separate implementation written for this article agree on the total; the Python bridge runs the C++ solver, so it is not an independent check.

One latent fault is worth recording, though the default never reaches it: the sieve bound is fixed at $30^5$ while parse_arguments only requires $k_{\max}\ge1$, so --k-max=31 makes euler_phi index past the end of the table — the address sanitizer reports that out-of-range read inside multiplicative_order. At the default it is clean.

Footnotes and References

The references supply the background facts used above; the reduction to a longest non-decreasing subsequence, the injectivity of the first $\ell+T$ stations, the bound on $\ell+T$ and the counts are derived here.

Project Euler 411 — Uphill Paths. The official statement puts stations at (2^i mod n, 3^i mod n) for 0 ≤ i ≤ 2n, counts equal coordinates once, and defines S(n) as the largest number of stations on a path from (0, 0) to (n, n) whose coordinates never decrease. It gives 11 stations with S(22) = 5, records S(123) = 14 and S(10000) = 48, and asks for the sum of S(k⁵) over 1 ≤ k ≤ 30.
Wikipedia — Mirsky's theorem. In a finite partially ordered set the height, the size of a longest chain, equals the smallest number of antichains that partition the set.
Wikipedia — Longest increasing subsequence. It is solvable in O(n log n) time by storing, for each length l, the smallest value ending an increasing subsequence of that length; those values increase, so one binary search per element suffices.
Wikipedia — Patience sorting. Each card goes on the leftmost pile whose top is at least as large, and the number of piles is the length of a longest increasing subsequence; Section 9 requires a strictly greater top, so equal values extend a non-decreasing subsequence.
Wikipedia — Chinese remainder theorem. For pairwise coprime moduli with product N, sending a residue modulo N to its residues modulo each factor is a ring isomorphism, so the factors can be handled separately.
Wikipedia — Multiplicative order. For a coprime to n, the order of a modulo n is the smallest k with a^k ≡ 1 (mod n); it divides Euler's totient φ(n), by Lagrange's theorem, and also the Carmichael function λ(n).
Wikipedia — Carmichael function. λ(n) is the least exponent that works for every unit modulo n; it divides φ(n), equals φ(n) for odd prime powers, and satisfies a | b ⇒ λ(a) | λ(b).
cp-algorithms — Linear sieve. It computes the minimum prime factor of every number in [2, n] in O(n) time, iterating over the primes up to the minimum prime factor of the current index and breaking there.
Wikipedia — Radix sort. A least-significant-digit radix sort works from the rightmost digit and preserves the ordering of earlier passes, so they must be stable; with d passes it costs O(d·n).
C++ — Euler411.cpp. The linked immutable revision contains build_spf, multiplicative_order, radix_sort_u64, S_value and solve_sum, and is the implementation record for this article.

Problem 411 source code

C++

#include <algorithm>
#include <array>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <string>
#include <vector>

namespace {

using i64 = long long;
using u64 = std::uint64_t;

struct Options {
    int k_max = 30;
    bool run_checkpoints = true;
};

bool parse_int_after_prefix(const std::string& arg, const std::string& prefix, int& value) {
    if (arg.rfind(prefix, 0U) != 0U) {
        return false;
    }
    const std::string tail = arg.substr(prefix.size());
    if (tail.empty()) {
        return false;
    }
    try {
        value = std::stoi(tail);
    } catch (...) {
        return false;
    }
    return true;
}

bool parse_arguments(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_int_after_prefix(arg, "--k-max=", options.k_max)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.k_max >= 1;
}

std::vector<int> build_spf(const int n) {
    std::vector<int> spf(static_cast<std::size_t>(n + 1), 0);
    std::vector<int> primes;
    primes.reserve(n / 10);
    for (int i = 2; i <= n; ++i) {
        if (spf[static_cast<std::size_t>(i)] == 0) {
            spf[static_cast<std::size_t>(i)] = i;
            primes.push_back(i);
        }
        for (const int p : primes) {
            const i64 x = 1LL * i * p;
            if (x > n) {
                break;
            }
            spf[static_cast<std::size_t>(x)] = p;
            if (p == spf[static_cast<std::size_t>(i)]) {
                break;
            }
        }
    }
    return spf;
}

void radix_sort_u64(std::vector<u64>& values) {
    if (values.size() <= 1U) {
        return;
    }

    std::vector<u64> tmp(values.size());
    constexpr int RADIX_BITS = 16;
    constexpr int RADIX_SIZE = 1 << RADIX_BITS;
    constexpr u64 RADIX_MASK = static_cast<u64>(RADIX_SIZE - 1);
    std::array<std::uint32_t, RADIX_SIZE> count{};

    for (int pass = 0; pass < 4; ++pass) {
        const int shift = pass * RADIX_BITS;
        count.fill(0);

        for (const u64 v : values) {
            ++count[static_cast<std::size_t>((v >> shift) & RADIX_MASK)];
        }

        std::uint32_t acc = 0;
        for (int i = 0; i < RADIX_SIZE; ++i) {
            const std::uint32_t c = count[static_cast<std::size_t>(i)];
            count[static_cast<std::size_t>(i)] = acc;
            acc += c;
        }

        for (const u64 v : values) {
            tmp[count[static_cast<std::size_t>((v >> shift) & RADIX_MASK)]++] = v;
        }
        values.swap(tmp);
    }
}

u64 mod_pow_u64(u64 base, u64 exp, u64 mod) {
    u64 result = 1 % mod;
    u64 cur = base % mod;
    u64 e = exp;
    while (e > 0) {
        if (e & 1ULL) {
            result = static_cast<u64>((__uint128_t)result * cur % mod);
        }
        cur = static_cast<u64>((__uint128_t)cur * cur % mod);
        e >>= 1ULL;
    }
    return result;
}

std::vector<int> unique_prime_factors(int x, const std::vector<int>& spf) {
    std::vector<int> fac;
    while (x > 1) {
        const int p = spf[static_cast<std::size_t>(x)];
        fac.push_back(p);
        while (x % p == 0) {
            x /= p;
        }
    }
    return fac;
}

int euler_phi(int x, const std::vector<int>& spf) {
    if (x == 1) {
        return 1;
    }
    int phi = x;
    int n = x;
    while (n > 1) {
        const int p = spf[static_cast<std::size_t>(n)];
        phi = phi / p * (p - 1);
        while (n % p == 0) {
            n /= p;
        }
    }
    return phi;
}

int multiplicative_order(int base, int mod, const std::vector<int>& spf) {
    if (mod == 1) {
        return 1;
    }
    int ord = euler_phi(mod, spf);
    std::vector<int> fac = unique_prime_factors(ord, spf);
    for (const int p : fac) {
        while (ord % p == 0) {
            const int cand = ord / p;
            if (mod_pow_u64(static_cast<u64>(base), static_cast<u64>(cand), static_cast<u64>(mod)) == 1ULL) {
                ord = cand;
            } else {
                break;
            }
        }
    }
    return ord;
}

int S_value(int n, const std::vector<int>& spf) {
    if (n == 1) {
        return 1;
    }

    int pre2 = 0;
    int tmp = n;
    while ((tmp & 1) == 0) {
        ++pre2;
        tmp >>= 1;
    }
    const int mod2 = tmp;
    const int per2 = multiplicative_order(2, mod2, spf);

    int pre3 = 0;
    tmp = n;
    while (tmp % 3 == 0) {
        ++pre3;
        tmp /= 3;
    }
    const int mod3 = tmp;
    const int per3 = multiplicative_order(3, mod3, spf);

    const int pre = std::max(pre2, pre3);
    const i64 period = std::lcm(static_cast<i64>(per2), static_cast<i64>(per3));
    const i64 needed = std::min<i64>(2LL * n + 1, pre + period);

    std::vector<u64> points;
    points.reserve(static_cast<std::size_t>(needed));

    int x = 1 % n;
    int y = 1 % n;
    for (i64 i = 0; i < needed; ++i) {
        const u64 key = (static_cast<u64>(static_cast<std::uint32_t>(x)) << 32) |
                        static_cast<u64>(static_cast<std::uint32_t>(y));
        points.push_back(key);
        x = static_cast<int>((2LL * x) % n);
        y = static_cast<int>((3LL * y) % n);
    }

    radix_sort_u64(points);
    points.erase(std::unique(points.begin(), points.end()), points.end());

    std::vector<int> tails;
    tails.reserve(256);
    for (const u64 key : points) {
        const int yy = static_cast<int>(key & 0xffffffffULL);
        auto it = std::upper_bound(tails.begin(), tails.end(), yy);
        if (it == tails.end()) {
            tails.push_back(yy);
        } else {
            *it = yy;
        }
    }

    return static_cast<int>(tails.size());
}

i64 solve_sum(int k_max, const std::vector<int>& spf) {
    i64 ans = 0;
    for (int k = 1; k <= k_max; ++k) {
        int n = 1;
        for (int t = 0; t < 5; ++t) {
            n *= k;
        }
        ans += S_value(n, spf);
    }
    return ans;
}

bool run_checkpoints(const std::vector<int>& spf) {
    if (S_value(22, spf) != 5) {
        std::cerr << "Checkpoint failed: S(22)\n";
        return false;
    }
    if (S_value(123, spf) != 14) {
        std::cerr << "Checkpoint failed: S(123)\n";
        return false;
    }
    if (S_value(10000, spf) != 48) {
        std::cerr << "Checkpoint failed: S(10000)\n";
        return false;
    }
    return true;
}

}  // namespace

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

    const int n_max = 24300000;  // 30^5
    const std::vector<int> spf = build_spf(n_max);

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

    std::cout << solve_sum(options.k_max, spf) << '\n';
    return 0;
}

Python

from __future__ import annotations

import re
import shutil
import subprocess
from pathlib import Path

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


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


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


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

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


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

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

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


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

Java

import java.util.*;

public class Euler411 {
    static int[] buildSPF(int n) {
        int[] spf = new int[n + 1];
        List<Integer> primes = new ArrayList<>();
        for (int i = 2; i <= n; i++) {
            if (spf[i] == 0) {
                spf[i] = i;
                primes.add(i);
            }
            for (int p : primes) {
                long x = (long) i * p;
                if (x > n)
                    break;
                spf[(int) x] = p;
                if (p == spf[i])
                    break;
            }
        }
        return spf;
    }

    static long modPow(long base, long exp, long mod) {
        long r = 1 % mod;
        base %= mod;
        while (exp > 0) {
            if ((exp & 1) == 1)
                r = r * base % mod;
            base = base * base % mod;
            exp >>= 1;
        }
        return r;
    }

    static int eulerPhi(int x, int[] spf) {
        int phi = x, n = x;
        while (n > 1) {
            int p = spf[n];
            phi = phi / p * (p - 1);
            while (n % p == 0)
                n /= p;
        }
        return phi;
    }

    static List<Integer> uniquePF(int x, int[] spf) {
        List<Integer> f = new ArrayList<>();
        while (x > 1) {
            int p = spf[x];
            f.add(p);
            while (x % p == 0)
                x /= p;
        }
        return f;
    }

    static int multOrder(int base, int mod, int[] spf) {
        if (mod == 1)
            return 1;
        int ord = eulerPhi(mod, spf);
        for (int p : uniquePF(ord, spf)) {
            while (ord % p == 0) {
                int cand = ord / p;
                if (modPow(base, cand, mod) == 1)
                    ord = cand;
                else
                    break;
            }
        }
        return ord;
    }

    static int sValue(int n, int[] spf) {
        if (n == 1)
            return 1;
        int pre2 = 0, tmp = n;
        while ((tmp & 1) == 0) {
            pre2++;
            tmp >>= 1;
        }
        int mod2 = tmp, per2 = multOrder(2, mod2, spf);
        int pre3 = 0;
        tmp = n;
        while (tmp % 3 == 0) {
            pre3++;
            tmp /= 3;
        }
        int mod3 = tmp, per3 = multOrder(3, mod3, spf);
        int pre = Math.max(pre2, pre3);
        long period = lcm(per2, per3);
        long needed = Math.min(2L * n + 1, pre + period);

        long[] points = new long[(int) needed];
        int x = 1 % n, y = 1 % n;
        for (int i = 0; i < needed; i++) {
            points[i] = ((long) x << 32) | (y & 0xFFFFFFFFL);
            x = (int) ((2L * x) % n);
            y = (int) ((3L * y) % n);
        }
        Arrays.sort(points);
        // Remove duplicates
        int uLen = 0;
        for (int i = 0; i < points.length; i++) {
            if (i == 0 || points[i] != points[i - 1])
                points[uLen++] = points[i];
        }
        // LIS on y-values
        int[] tails = new int[uLen];
        int lisLen = 0;
        for (int i = 0; i < uLen; i++) {
            int yy = (int) (points[i] & 0xFFFFFFFFL);
            int pos = upperBound(tails, lisLen, yy);
            if (pos == lisLen)
                tails[lisLen++] = yy;
            else
                tails[pos] = yy;
        }
        return lisLen;
    }

    static int upperBound(int[] arr, int len, int val) {
        int lo = 0, hi = len;
        while (lo < hi) {
            int mid = (lo + hi) >>> 1;
            if (arr[mid] <= val)
                lo = mid + 1;
            else
                hi = mid;
        }
        return lo;
    }

    static long gcd(long a, long b) {
        while (b != 0) {
            long t = b;
            b = a % b;
            a = t;
        }
        return a;
    }

    static long lcm(long a, long b) {
        return a / gcd(a, b) * b;
    }

    public static String solve() {
        int kMax = 30;
        int nMax = 1;
        for (int t = 0; t < 5; t++)
            nMax *= kMax;
        int[] spf = buildSPF(nMax);
        long ans = 0;
        for (int k = 1; k <= kMax; k++) {
            int n = 1;
            for (int t = 0; t < 5; t++)
                n *= k;
            ans += sValue(n, spf);
        }
        return String.valueOf(ans);
    }

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