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 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$. The condition for a future value $x$ is

$$\varphi(i) \lt \varphi(x) \lt i,$$

which depends on $x$ only through the single value $\varphi(x)$. Since totients are integers, this is equivalent to

$$\varphi(x)\in[\varphi(i)+1,\; i-1].$$

So once $dp[i]$ is known, the same weight $dp[i]$ should be added to every future number whose totient falls inside that interval. This observation eliminates the need to inspect every predecessor separately. It also shows immediately why prime endpoints are dead ends: if $i$ is prime, then $\varphi(i)=i-1$, so the interval is empty.

Accumulator on the Totient Axis

After all endpoints smaller than $x$ have been processed, define

$$A_x(t)=\sum_{i=6}^{x-1} dp[i]\,\mathbf{1}_{\varphi(i)+1\le t \le i-1}.$$

This means $A_x(t)$ stores the total contribution of all earlier endpoints that would accept a future number whose totient equals $t$. Evaluating at $t=\varphi(x)$ gives

$$dp[x]=A_x(\varphi(x)).$$

So the problem becomes: maintain many interval additions on the totient axis, and answer one point query for each $x$. Because $\varphi(x)\le x-1$ for every $x \gt 1$, a Fenwick structure of size $N$ is sufficient.

Fenwick Tree as a Difference Array

The implementation uses a Fenwick tree, not for direct prefix sums of $dp$, but for the difference array of the accumulator $A_x$. If we want to add a value $v$ to every index in an interval $[L,R]$, we can update a difference array $D$ by

$$D[L]\mathrel{+}=v,\qquad D[R+1]\mathrel{-}=v.$$

Then the actual value at position $t$ is the prefix sum

$$A_x(t)=\sum_{u=1}^{t} D[u].$$

A Fenwick tree supports those point updates and prefix-sum queries in $O(\log N)$. In other words, every finished state $i$ performs one interval update on $[\varphi(i)+1,\; i-1]$, and every candidate endpoint $x$ performs one point query at $\varphi(x)$.

Totient Precomputation

The DP cannot even start before all values $\varphi(1),\dots,\varphi(N)$ are known. These are computed with the standard totient sieve. Initialize $\varphi[m]=m$ for all $m$, and for every prime $p$, visit all multiples $kp$ and apply

$$\varphi(kp)\leftarrow \varphi(kp)-\frac{\varphi(kp)}{p}.$$

This is the batch version of the multiplicative formula $\varphi(n)=n\prod_{p\mid n}(1-\frac1p)$. The asymptotic cost is $O(N\log\log N)$. The C++ implementation uses block partitioning and optional multithreading for faster preprocessing, but the mathematical result is the same as the ordinary sieve.

Worked Example: $N=10$

The relevant totient values are

$$\varphi(6)=2,\quad \varphi(7)=6,\quad \varphi(8)=4,\quad \varphi(9)=6,\quad \varphi(10)=4.$$

Start with $dp[6]=1$, so the total already contains the sequence $\{6\}$. Since a successor of 6 must satisfy $2 \lt \varphi(x) \lt 6$, we seed the interval $[3,5]$ with weight 1.

Now iterate upward:

$x=7$: query $\varphi(7)=6$, obtain $dp[7]=0$, so no valid sequence ends at 7.

$x=8$: query $\varphi(8)=4$, obtain $dp[8]=1$, corresponding to $\{6,8\}$. Then 8 contributes to $[\varphi(8)+1,7]=[5,7]$.

$x=9$: query $\varphi(9)=6$, obtain $dp[9]=1$, corresponding to $\{6,8,9\}$.

$x=10$: query $\varphi(10)=4$, obtain $dp[10]=1$, corresponding to $\{6,10\}$.

Hence

$$S(10)=1+1+1+1=4,$$

namely the sequences $\{6\}$, $\{6,8\}$, $\{6,8,9\}$, and $\{6,10\}$.

Why the Code Starts with a Single Range Update

The code sets total = 1 for the trivial sequence $[6]$, then immediately performs

rangeAdd(phi[6] + 1, 5, 1).

Because $\varphi(6)=2$, this is exactly the interval $[3,5]$, which represents every totient value allowed for a direct successor of 6. After that, each loop iteration follows the same pattern:

1. Query the current accumulator at $\varphi(x)$ to obtain $dp[x]$.

2. Add $dp[x]$ to the running total.

3. Propagate $dp[x]$ to the interval $[\varphi(x)+1,\; x-1]$ for future successors.

If the interval is empty, the helper simply skips the update.

How the Code Works

The C++ solution follows the derivation above literally. It first computes all totients, then scans $x$ from 7 to $N$. Each iteration performs one Fenwick point query and one range update, both reduced modulo $10^8$. The C++ file also includes a brute-force checker for small $N$ and validates the official checkpoints $S(10)=4$, $S(100)=482073668$, and $S(10000)\bmod 10^8 = 73808307$. The Java version implements the same algorithm directly. The Python file is a thin bridge that compiles and runs the C++ solver so all three language entries stay consistent.

Complexity Analysis

The totient sieve costs $O(N\log\log N)$. The DP stage performs $N-6$ point queries and at most $N-6$ interval updates, each in $O(\log N)$, so the dominant cost is $O(N\log N)$. Memory usage is $O(N)$ for the $\varphi$ array and the Fenwick structure. This is a dramatic improvement over the naive DP, which would require $O(N^2)$ predecessor checks.

Footnotes and References

Problem page: https://projecteuler.net/problem=337
Euler totient function: https://en.wikipedia.org/wiki/Euler%27s_totient_function
Fenwick tree (Binary Indexed Tree): https://en.wikipedia.org/wiki/Fenwick_tree
Sieve of Eratosthenes and sieve-style preprocessing: https://en.wikipedia.org/wiki/Sieve_of_Eratosthenes

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