Problem 407: Idempotents

Project Euler Problem 407 Solution

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

Problem Summary For each integer $n \le N$, define $M(n)$ as the largest residue $a \lt n$ satisfying $$a^2 \equiv a \pmod{n}.$$ The goal is to compute $$\sum_{n=1}^{N} M(n)$$ for the large bound $N=10^7$. The key observation is that these residues are exactly the idempotents modulo $n$, and their structure is controlled by the prime-power factorization of $n$. Mathematical Approach Step 1: Rewrite the Congruence The defining condition is equivalent to $$a^2 \equiv a \pmod{n}\iff a(a-1)\equiv 0 \pmod{n}.$$ This form is more informative because the two factors differ by $1$, so $$\gcd(a,a-1)=1.$$ Therefore every prime power dividing $n$ must divide one of these two coprime factors completely. Step 2: Prime-Power Blocks Force Local Choices Write the factorization of $n$ as $$n=\prod_{i=1}^{k} q_i,\qquad q_i=p_i^{e_i},$$ where the numbers $q_1,\dots,q_k$ are pairwise coprime prime powers. If $q_i\mid a(a-1)$, then the prime underlying $q_i$ cannot divide both $a$ and $a-1$, so the whole block $q_i$ must divide exactly one of them....

Detailed mathematical approach

Problem Summary

Call a residue $a$ with $0\le a\lt n$ idempotent modulo $n$ when

$$a^2\equiv a\pmod n,$$

and let $M(n)$ be the largest idempotent residue modulo $n$. The statement works out $n=6$: the values of $a^2\bmod 6$ for $a=0,1,\dots,5$ are $0,1,4,3,4,1$, and they agree with $a$ itself exactly for $a=0,1,3,4$, so $M(6)=4$. We must compute

$$\sum_{n=1}^{N}M(n),\qquad N=10^{7}.$$

The case $n=1$ is degenerate: the only residue below $1$ is $a=0$, and $0^2\equiv 0\pmod 1$, so $M(1)=0$. ^[1]

The whole solution in one sentence: a residue is idempotent modulo $n$ exactly when it is $0$ or $1$ modulo every prime power occurring in $n$, so each splitting of $n$ into coprime factors $n=uv$ yields one idempotent, and for every $n$ we build all of them and keep the largest.

Mathematical Approach

Throughout, $\omega(n)$ denotes the number of distinct primes dividing $n$. We show that the idempotents modulo $n$ correspond to the coprime splittings of $n$, give a closed formula for the idempotent of one splitting, and use a symmetry to bound $M(n)$. The computation is then a smallest-prime-factor sieve, which yields the blocks by repeated division, plus a loop over the $2^{\omega(n)}$ splittings of each $n$.

1. From a congruence to divisibility

Subtracting $a$ from both sides, $a^2\equiv a\pmod n$ says precisely that

$$n\mid a(a-1).$$

The two factors are consecutive integers, so they are coprime: a common divisor of $a$ and $a-1$ divides their difference $1$. No prime of $n$ can divide both, which is the source of all the structure below.

2. Prime-power blocks decide locally

By the uniqueness of prime factorisation write

$$n=q_1q_2\cdots q_k,\qquad q_i=p_i^{e_i},\ k=\omega(n),$$

with distinct primes $p_i$. We call $q_1,\dots,q_k$ the blocks of $n$; they are pairwise coprime.

Claim. A residue $a$ satisfies $n\mid a(a-1)$ if and only if for every $i$ either $a\equiv 0$ or $a\equiv 1$ modulo $q_i$.

Proof. Suppose $n\mid a(a-1)$ and fix $i$. Then $q_i\mid a(a-1)$. The prime $p_i$ cannot divide both $a$ and $a-1$; say it misses $a$ (the other case is symmetric). Then $\gcd(q_i,a)=1$, and $q_i\mid a(a-1)$ forces $q_i\mid a-1$, that is $a\equiv 1\pmod{q_i}$. Conversely, if each $q_i$ divides $a$ or $a-1$, then each $q_i$ divides the product $a(a-1)$; as the $q_i$ are pairwise coprime, their product $n$ divides it too. End of proof.

The divisibility has become a free and independent choice of $0$ or $1$ on each of the $k$ blocks.

3. Counting the idempotents

The Chinese remainder theorem says that for pairwise coprime moduli $q_1,\dots,q_k$ with product $n$, every choice of residues modulo the $q_i$ is realised by exactly one $a$ with $0\le a\lt n$. ^[3] Applied to the $2^k$ patterns of Section 2 it gives exactly $2^{\omega(n)}$ idempotents modulo $n$, all distinct. For $n=6$ the blocks are $2$ and $3$, and the four patterns give $0,1,3,4$, precisely the residues of the statement's table. For $n=12=4\cdot 3$ they are $0,1,4,9$, so $M(12)=9$.

It is convenient to record a pattern by the set $U$ of blocks on which $a\equiv 0$. Putting

$$u=\prod_{i\in U}q_i,\qquad v=\frac nu,$$

we get $uv=n$ and $\gcd(u,v)=1$, since no prime of $n$ appears in both products; conversely each such factorisation comes from exactly one $U$, as every block lies wholly in $u$ or wholly in $v$. These $u$ are the unitary divisors of $n$, and there are $2^{\omega(n)}$ of them. ^[4]

4. The idempotent of one splitting

Let $n\ge 2$ and fix a splitting $n=uv$ with $\gcd(u,v)=1$. We must solve

$$a\equiv 0\pmod u,\qquad a\equiv 1\pmod v.$$

The first congruence says $a=ut$, and the second then reads $ut\equiv 1\pmod v$. Since $\gcd(u,v)=1$, the residue $u\bmod v$ is invertible, so $t$ is determined modulo $v$. Taking the representative $0\le t\lt v$ gives

$$a=u\,t,\qquad t\equiv u^{-1}\pmod v,$$

and $0\le a\le u(v-1)\lt n$ automatically, so no reduction modulo $n$ is needed. The inverse comes from the extended Euclidean algorithm, which produces $x,y$ with $ux+vy=1$ in $O(\log v)$ divisions; then $t=x\bmod v$. ^[5]

Two splittings are degenerate: $u=1$, $v=n$ gives $t=1$ and $a=1$, while $u=n$, $v=1$ gives $t=0$ and $a=0$. These two trivial idempotents exist for every $n\ge 2$; any other is called nontrivial. For $n=1$ the two coincide in the single residue $0$.

5. A worked example: $n=60$

Here $60=2^2\cdot 3\cdot 5$, so the blocks are $4$, $3$ and $5$ and there are $2^3=8$ splittings. Listing them as $u$, $v$, $t$, $a=ut$:

$u=1$, $v=60$: $t=1$, $a=1$.
$u=4$, $v=15$: $t=4$, $a=16$.
$u=3$, $v=20$: $t=7$, $a=21$.
$u=5$, $v=12$: $t=5$, $a=25$.
$u=12$, $v=5$: $t=3$, $a=36$.
$u=20$, $v=3$: $t=2$, $a=40$.
$u=15$, $v=4$: $t=3$, $a=45$.
$u=60$, $v=1$: $t=0$, $a=0$.

The idempotents modulo $60$ are thus $0,1,16,21,25,36,40,45$ and $M(60)=45$; the largest one checks directly, $45^2=2025=33\cdot 60+45$. Reading one line, say $u=12$, $v=5$: the inverse of $12\bmod 5=2$ modulo $5$ is $3$, and $a=12\cdot 3=36$ is divisible by $12$ and one more than a multiple of $5$.

6. The complementary idempotent

If $a$ is idempotent then so is $1-a$, because

$$(1-a)^2=1-2a+a^2\equiv 1-2a+a=1-a.$$

This map is an involution, so it permutes the idempotents. ^[2] For nontrivial $a$ we have $2\le a\le n-1$, so the representative of $1-a$ is $n+1-a$, again nontrivial. In the language of Section 4 the involution merely exchanges $u$ and $v$: as $n$ vanishes modulo both, $a\equiv 0\pmod u$ and $a\equiv 1\pmod v$ become $n+1-a\equiv 1\pmod u$ and $n+1-a\equiv 0\pmod v$. The $n=60$ list shows the pairs $16+45$, $21+40$, $25+36$, each $61=n+1$; for $n=6$ it is $3+4=7$.

Write $m(n)$ for the smallest nontrivial idempotent. When $\omega(n)\ge 2$, the pairing turns the largest nontrivial idempotent into the smallest, so

$$M(n)=n+1-m(n).$$

For each complementary pair, therefore, one inverse suffices: compute $a$, take its partner $n+1-a$, and compare both with the current maximum.

7. Prime powers contribute nothing but $1$

If $\omega(n)=1$ there is a single block $q_1=n$, so the only patterns are $a\equiv 0$ and $a\equiv 1$ modulo $n$: the idempotents are $0$ and $1$, and $M(n)=1$ for every prime power $n\ge 2$. Up to $10^{7}$ there are $665134$ such $n$, namely the $664579$ primes and $555$ higher powers, so about $6.7\%$ of all $n$ contribute only $1$. Every other $n\ge 2$ has nontrivial idempotents and hence $M(n)\ge 2$.

8. How large $M(n)$ can be

Let $a$ be nontrivial, coming from $n=uv$ with $u,v\ge 2$. From $u\mid a$ and $a\ge 1$ we get $a\ge u$; from $v\mid a-1$ and $a-1\ge 1$ we get $a\ge v+1$. Multiplying,

$$a^2\ge u(v+1)\gt uv=n,$$

so every nontrivial idempotent exceeds $\sqrt n$. Hence $m(n)\ge\lfloor\sqrt n\rfloor+1$ and, for $\omega(n)\ge 2$,

$$M(n)\le n-\lfloor\sqrt n\rfloor .$$

Equality holds exactly for the oblong numbers $n=s(s+1)$, $s\ge 2$. ^[8] Indeed, writing $s=\lfloor\sqrt n\rfloor$, equality means $m(n)=s+1$, hence $n\mid (s+1)s$; every divisor of $s(s+1)$ except $s(s+1)$ itself is at most $s(s+1)/2$, which is smaller than $s^2$ once $s\ge 2$, while $n\ge s^2$ — so $n=s(s+1)$. Conversely, for $n=s(s+1)$ the factors $u=s+1$ and $v=s$ are coprime and at least $2$, and $u\equiv 1\pmod v$ gives $t=1$ and $a=s+1$, the least value allowed. Then

$$M\bigl(s(s+1)\bigr)=s(s+1)-s=s^2,$$

which produces the squares $M(6)=4$, $M(12)=9$, $M(20)=16$, $M(30)=25$, $M(42)=36$.

9. First values, and why a search is needed

For $n=1,\dots,12$ the values are $0,1,1,1,1,4,1,1,1,6,1,9$, so the sum up to $10$ is $17$; the sequence is catalogued. ^[7]

The number of distinct prime factors alone does not determine $M(n)$: the nearby numbers $60$, $66=2\cdot 3\cdot 11$ and $84=2^2\cdot 3\cdot 7$ all have $\omega=3$, yet $m(60)=16$, $m(66)=12$ and $m(84)=21$, giving $M=45$, $55$ and $64$. The minimum depends on how the individual inverses fall, and no simpler rule for it is known, so for each $n$ that is not a prime power we evaluate all $2^{\omega(n)}-2$ nontrivial splittings — far fewer than the $n$ residues a direct search would test.

How the Code Works

All three programs follow Sections 2 to 5 literally. build_spf (buildSpf in Java, inline in Python) runs a linear sieve: for $i=2,\dots,n$ it records a new prime when spf[i] is still zero, then walks the primes found so far, sets spf[i*p] = p and breaks as soon as p == spf[i]. ^[6] The break is what makes the sieve linear: each composite is written once, by its smallest prime factor. C++ and Java form i*p in 64 bits before comparing with $n$, so the index cannot overflow.

The main loop of sum_M (solve in Java and Python) lets $n=1$ contribute $0$ and then peels off blocks: it reads $p=$ spf[x], divides $x$ by $p$ while it divides and stores the full power pe, so prime_powers[0..k-1] holds the blocks $q_i$ in increasing order of their primes. If $k$ is $1$ it adds $1$ at once, which is Section 7.

Otherwise it builds subset_prod, the table of the $2^k$ subset products $u$: the lowest set bit of mask is isolated as mask & -mask, its index found with __builtin_ctz in C++, Integer.numberOfTrailingZeros in Java and lsb.bit_length() - 1 in Python, and subset_prod[mask] is subset_prod[mask ^ lsb] times that block — one multiplication per nonempty subset, the empty product being initialised to $1$.

The candidate loop runs over mask from $1$ to full - 2, the nonempty proper subsets, skipping exactly the two trivial splittings. For each it takes $u=$ subset_prod[mask], $v=n/u$, calls inv_mod(u % v, v) and forms $a$; best starts at the trivial value $1$ and ends as $M(n)$.

inv_mod (invMod in Java) is the iterative extended Euclidean algorithm carrying only the Bézout coefficient of its first argument: it replaces $a$ by $a-tb$ with $t=\lfloor a/b\rfloor$, swaps, and treats the coefficient pair the same way. C++ and Java add the modulus once if the final coefficient is negative; Python reduces it with u % m instead. ^[5] The reduction u % v is harmless, the inverse depending only on $u$ modulo $v$, and it cannot give $0$, since $\gcd(u,v)=1$ and $v\ge 2$.

Integer widths differ. C++ forms $u\cdot t$ as 1LL * u * inv, although $ut\lt n\le 10^7$ makes 32 bits enough, and keeps the total in a uint64_t. Java accumulates into a long, which is safe because $M(n)\lt n$ gives the crude bound $\sum M(n)\lt N^2/2=5\cdot 10^{13}$; it maps each $n$ to its $M(n)$ with IntStream.rangeClosed(1, limit), mapToLong and sum, allocating its two small arrays inside the mapping function. Python needs no such analysis, its integers being arbitrary precision.

Only the C++ program takes arguments: parse_arguments understands --limit= and --skip-checkpoints, and run_checkpoints recovers the statement's value as sum_M(6) - sum_M(5), comparing it with $4$, and checks $M(1)+M(2)=1$. The Python and Java ports fix the bound at $10^7$ and carry no self-test.

Complexity and Verification

The sieve is $O(N)$ in time, and the auxiliary storage stays $O(N)$ machine words, dominated by the table of $N+1$ entries and the list of the $664579$ primes: with 32-bit entries, $40$ MB, and the measured peak resident size of the C++ run is about $45$ MB. Python's lists of references to integer objects bring it to about $122$ MB.

The work for one $n$ is $O(2^{\omega(n)}\log n)$: extracting the blocks takes $O(\log n)$ divisions, and each subset costs one multiplication and, apart from the two trivial ones, one modular inverse. Now $2^{\omega(n)}$ is also the number of squarefree divisors of $n$, a squarefree divisor being a product of a subset of the distinct primes, so counting the pairs $d\mid n$ in the other order,

$$\sum_{n\le N}2^{\omega(n)}=\sum_{d\le N}\mu^2(d)\lfloor N/d\rfloor,$$

where $\mu^2(d)$ is $1$ for squarefree $d$ and $0$ otherwise. The right-hand side is at most $N\sum_{d\le N}1/d=O(N\log N)$, and the positive density $6/\pi^2$ of the squarefree numbers ^[9] makes it $\Theta(N\log N)$. At $N=10^{7}$ both sides equal $105854997$, about $10.6$ splittings per $n$. For the $9334865$ values of $n$ with $\omega(n)\ge 2$ the code fills $104524728$ table entries, of which $95189863$ are multiplications and the rest the initial $1$, and it performs $85854998$ modular inverses. With the $O(\log n)$ bound per inverse the computation is $O(N\log^2 N)$ operations. ^[5] The run takes $460941721$ Euclidean iterations, $5.37$ per inverse; the classical average-case estimate assumes uniformly random inputs and predicts $6.46$ here, so it is only a rough guide for these structured pairs.

For $n\le 10^{7}$ the exponent $\omega(n)$ never exceeds $8$, and $8$ occurs for the single number $9699690$, the product of the first eight primes, the next primorial being $223092870$. Over $2\le n\le 10^{7}$ the counts by $\omega=1,2,\dots,8$ are $665134$, $2536838$, $3642766$, $2389433$, $691209$, $72902$, $1716$ and $1$, so no $n$ needs more than $2^8-2=254$ candidates. Compiled with -O2, the C++ program finishes in about $2.0$ seconds and the Python one in about $44$ seconds on the same machine; both figures are machine-dependent.

Two refinements follow from the mathematics above and are not used by the code. Restricting the nonempty proper subsets to those that keep one fixed block on the $u$ side examines one member of each complementary pair and reads off both $a$ and $n+1-a$, which halves the inverses to $42927499$. Taking instead the $k$ primitive idempotents $e_i$, each $1$ modulo $q_i$ and $0$ modulo the other blocks and so given by Section 4 with $u=n/q_i$, needs one inverse per block, $29465183$ in all, after which every idempotent is a subset sum $\sum_{i\notin U}e_i\bmod n$, because on each block only its own $e_i$ is $1$. Both variants reproduce $M(n)$ for every $n\le 60000$.

Beyond the C++ checkpoint, the following was verified for this article. A brute force over every $a\lt n$ reproduces $M(n)$ for all $n\le 3000$, where the idempotent count also equals $2^{\omega(n)}$, and the first twenty values match the catalogued sequence. ^[7] The map $a\mapsto n+1-a$ permutes the nontrivial idempotents for $n\le 20000$; every nontrivial idempotent exceeds $\sqrt n$ for $n\le 10^{5}$; the bound of Section 8 holds for $n\le 2\cdot 10^{5}$, with equality at exactly the $445$ oblong numbers there. The partial sums are $17$, $2549$, $314034$, $34981569$, $3717852515$ and $386778427463$ for $N=10,10^2,\dots,10^6$, and the three programs' outputs at the full bound are identical.

Footnotes and References

The bijection with coprime splittings, the complement symmetry, the bound involving $\sqrt n$ with its extremal cases, and the operation counts are derived in this article; the references supply the background facts.

Project Euler 407 — Idempotents. The official statement lists a² mod 6 for 0 ≤ a ≤ 5 as 0, 1, 4, 3, 4, 1, defines M(n) as the largest a < n with a² ≡ a mod n, states M(6) = 4, and asks for the sum of M(n) over 1 ≤ n ≤ 10⁷.
Wikipedia — Idempotent (ring theory). An element with a² = a is idempotent, and if a is idempotent so is 1 − a. It works out the integers modulo 6, whose idempotents are 0, 1, 3 and 4 with 3 + 4 ≡ 1 mod 6, and notes that a product of m factors having only 0 and 1 as idempotents has 2^m idempotents.
Wikipedia — Chinese remainder theorem. For pairwise coprime moduli with product N, the system x ≡ a_i mod n_i has exactly one solution with 0 ≤ x < N, and reduction is a ring isomorphism onto the product of the integers modulo each n_i.
Wikipedia — Unitary divisor. A divisor d of n is unitary when d and n/d share no factor greater than 1; since each is then a product of full prime powers over a subset of the primes, their number is 2 raised to the number of distinct prime factors. OEIS A034444 tabulates it; its first terms 1, 2, 2, 2, 2, 4, 2, 2, 2, 4, 2, 4 count the idempotents modulo 1, …, 12.
Wikipedia — Extended Euclidean algorithm. Besides the greatest common divisor the algorithm yields the Bézout coefficients; for coprime inputs the coefficient of the first, reduced modulo the second, is the modular inverse. For the number of divisions, Wikipedia — Euclidean algorithm gives Lamé's worst case, at most five times the decimal digit count of the smaller input, and an average of (12/π²) ln 2 · ln n + 0.06 for inputs drawn uniformly up to n — a distribution the coprime factor pairs here do not follow.
cp-algorithms — Linear sieve. It computes the minimum prime factor of every integer in the range, and the list of primes, in O(n) time: the inner loop stops once the current prime equals the minimum prime factor of i, so each composite is assigned once.
OEIS A182665. Defined as the greatest x < n with n dividing x(x − 1), which is M(n); its terms for n = 1, …, 20 are 0, 1, 1, 1, 1, 4, 1, 1, 1, 6, 1, 9, 1, 8, 10, 1, 1, 10, 1, 16. It notes that the value is 1 exactly for prime powers, and records M(n) = n − A344005(n) with A344005(n) the least positive m with n dividing m(m + 1) — the symmetry of Section 6.
OEIS A002378. The oblong (or pronic) numbers are the products of consecutive integers, s(s + 1), beginning 0, 2, 6, 12, 20, 30, 42, 56, 72, 90; Section 8 shows that those with s ≥ 2 are exactly the n with ω(n) ≥ 2 attaining M(n) = n − ⌊√n⌋.
Wikipedia — Square-free integer. The count of squarefree integers up to x is 6x/π² + O(√x), so their natural density is 6/π² ≈ 0.6079, which is the fact used in the complexity section.

Problem 407 source code

C++

#include <cstdint>
#include <iostream>
#include <string>
#include <vector>

namespace {

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

struct Options {
    int limit = 10000000;
    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, "--limit=", options.limit)) {
            continue;
        }
        std::cerr << "Unknown argument: " << arg << '\n';
        return false;
    }
    return options.limit >= 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;
}

int inv_mod(int a, int mod) {
    int b = mod;
    int u = 1;
    int v = 0;
    while (b != 0) {
        const int t = a / b;
        a -= t * b;
        std::swap(a, b);
        u -= t * v;
        std::swap(u, v);
    }
    if (u < 0) {
        u += mod;
    }
    return u;
}

u64 sum_M(const int limit) {
    const std::vector<int> spf = build_spf(limit);

    std::vector<int> prime_powers(10, 0);
    std::vector<int> subset_prod(1 << 10, 1);

    u64 ans = 0ULL;
    for (int n = 1; n <= limit; ++n) {
        if (n == 1) {
            continue;
        }

        int x = n;
        int k = 0;
        while (x > 1) {
            const int p = spf[static_cast<std::size_t>(x)];
            int pe = 1;
            while (x % p == 0) {
                x /= p;
                pe *= p;
            }
            prime_powers[static_cast<std::size_t>(k++)] = pe;
        }

        if (k == 1) {
            ans += 1ULL;
            continue;
        }

        const int full = (1 << k);
        subset_prod[0] = 1;
        for (int mask = 1; mask < full; ++mask) {
            const int lsb = mask & -mask;
            const int bit = __builtin_ctz(static_cast<unsigned>(lsb));
            subset_prod[mask] = subset_prod[mask ^ lsb] * prime_powers[static_cast<std::size_t>(bit)];
        }

        int best = 1;
        for (int mask = 1; mask < full - 1; ++mask) {
            const int u = subset_prod[mask];
            const int v = n / u;
            const int inv = inv_mod(u % v, v);
            const int a = static_cast<int>((1LL * u * inv) % n);
            if (a > best) {
                best = a;
            }
        }

        ans += static_cast<u64>(best);
    }

    return ans;
}

bool run_checkpoints() {
    // From statement: M(6)=4.
    const u64 s6 = sum_M(6) - sum_M(5);
    if (s6 != 4ULL) {
        std::cerr << "Checkpoint failed: M(6)\n";
        return false;
    }
    // Small consistency: M(1)=0 and M(prime)=1.
    const u64 s2 = sum_M(2);
    if (s2 != 1ULL) {
        std::cerr << "Checkpoint failed: sum up to 2\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_checkpoints && !run_checkpoints()) {
        return 2;
    }

    std::cout << sum_M(options.limit) << '\n';
    return 0;
}

Python

def solve():
    LIMIT = 10_000_000

    # Build SPF via linear sieve
    spf = [0] * (LIMIT + 1)
    primes = []
    for i in range(2, LIMIT + 1):
        if spf[i] == 0:
            spf[i] = i
            primes.append(i)
        for p in primes:
            if i * p > LIMIT: break
            spf[i * p] = p
            if p == spf[i]: break

    def inv_mod(a, m):
        b, u, v = m, 1, 0
        while b:
            t = a // b
            a, b = b, a - t * b
            u, v = v, u - t * v
        return u % m

    ans = 0
    for n in range(2, LIMIT + 1):
        x = n
        k = 0
        pp = []
        while x > 1:
            p = spf[x]
            pe = 1
            while x % p == 0:
                x //= p
                pe *= p
            pp.append(pe)
            k += 1

        if k == 1:
            ans += 1
            continue

        full = 1 << k
        subset_prod = [0] * full
        subset_prod[0] = 1
        for mask in range(1, full):
            lsb = mask & -mask
            bit = lsb.bit_length() - 1
            subset_prod[mask] = subset_prod[mask ^ lsb] * pp[bit]

        best = 1
        for mask in range(1, full - 1):
            u = subset_prod[mask]
            v = n // u
            inv = inv_mod(u % v, v)
            a = (u * inv) % n
            if a > best:
                best = a

        ans += best

    return str(ans)

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

Java

import java.util.ArrayList;
import java.util.List;
import java.util.stream.IntStream;

public class Euler407 {
    static int[] buildSpf(int n) {
        int[] spf = new int[n + 1];
        List<Integer> primes = new ArrayList<>(n / 10);

        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 int invMod(int a, int mod) {
        int b = mod;
        int u = 1;
        int v = 0;
        while (b != 0) {
            int t = a / b;
            int tempA = a - t * b;
            a = b;
            b = tempA;
            int tempU = u - t * v;
            u = v;
            v = tempU;
        }
        if (u < 0)
            u += mod;
        return u;
    }

    static String solve() {
        int limit = 10000000;
        int[] spf = buildSpf(limit);

        long total = IntStream.rangeClosed(1, limit)
                .parallel()
                .mapToLong(n -> {
                    if (n == 1)
                        return 0;

                    int[] primePowers = new int[10];
                    int x = n;
                    int k = 0;

                    while (x > 1) {
                        int p = spf[x];
                        int pe = 1;
                        while (x % p == 0) {
                            x /= p;
                            pe *= p;
                        }
                        primePowers[k++] = pe;
                    }

                    if (k == 1)
                        return 1;

                    int full = 1 << k;
                    int[] subsetProd = new int[full];
                    subsetProd[0] = 1;

                    for (int mask = 1; mask < full; mask++) {
                        int lsb = mask & -mask;
                        int bit = Integer.numberOfTrailingZeros(lsb);
                        subsetProd[mask] = subsetProd[mask ^ lsb] * primePowers[bit];
                    }

                    int best = 1;
                    for (int mask = 1; mask < full - 1; mask++) {
                        int u = subsetProd[mask];
                        int v = n / u;
                        int inv = invMod(u % v, v);
                        int a = (int) (((long) u * inv) % n);
                        if (a > best)
                            best = a;
                    }

                    return best;
                })
                .sum();

        return Long.toString(total);
    }

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