Solution

Answer: 39782849136421

For an integer $a$, the congruence

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

means

$$a(a-1)\equiv 0 \pmod n.$$

The key observation is that $a$ and $a-1$ are coprime, so every prime-power divisor of $n$ must divide exactly one of them.

1. Mathematical analysis

We seek

$$M(n)=\max{a<n : a^2\equiv a \pmod n}.$$

We must compute

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

Step 1 — Structure of all solutions

Let

$$n=\prod_{i=1}^k p_i^{e_i}.$$

Because

$$a(a-1)\equiv 0 \pmod{p_i^{e_i}},$$

and $\gcd(a,a-1)=1$, we must have for each prime power:

$$a\equiv 0 \pmod{p_i^{e_i}} \quad\text{or}\quad a\equiv 1 \pmod{p_i^{e_i}}.$$

Thus every solution corresponds to choosing, independently for each prime power factor, either $0$ or $1$.

By the Chinese Remainder Theorem, there are exactly

$$2^{\omega(n)}$$

solutions, where $\omega(n)$ is the number of distinct prime factors.

Step 2 — Rewriting with coprime factors

Partition the prime powers of $n$ into two coprime parts:

$$n = uv,\qquad \gcd(u,v)=1.$$

Then the corresponding idempotent satisfies

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

Write $a=uk$. Then

$$uk\equiv 1\pmod v.$$

$$k\equiv u^{-1}\pmod v.$$

Hence the unique solution modulo $n$ is

$$a=u\cdot u^{-1}\pmod v.$$

Every idempotent arises this way.

Therefore

$$M(n)=\max_{uv=n,\ \gcd(u,v)=1} \left(u\cdot u^{-1}\pmod v\right).$$

Step 3 — Efficient computation

A brute force factorization for every $n\le10^7$ is too slow.

Instead:

Build smallest-prime-factor sieve.
Recover distinct prime factors for each $n$.
Generate all subsets of distinct prime powers.
For each coprime split $n=uv$:

compute modular inverse $u^{-1}\pmod v$,
construct the idempotent,
maximize.

Because $10^7$ has at most 8 distinct prime factors, the subset enumeration is tiny.

The total complexity is essentially linear-sieve scale.

2. Python implementation

from math import prod

LIMIT = 10_000_000

# ---------------------------------------------------------
# Smallest prime factor sieve
# ---------------------------------------------------------
spf = list(range(LIMIT + 1))

for i in range(2, int(LIMIT ** 0.5) + 1):
    if spf[i] == i:
        step = i
        start = i * i
        for j in range(start, LIMIT + 1, step):
            if spf[j] == j:
                spf[j] = i

# ---------------------------------------------------------
# Extended Euclidean algorithm
# ---------------------------------------------------------
def invmod(a, m):
    """Return inverse of a modulo m."""
    x0, x1 = 1, 0
    b = m

    while b:
        q = a // b
        a, b = b, a % b
        x0, x1 = x1, x0 - q * x1

    return x0 % m

# ---------------------------------------------------------
# Factorization into prime powers
# Example:
#   12 -> [4, 3]
#   72 -> [8, 9]
# ---------------------------------------------------------
def prime_powers(n):
    res = []

    while n > 1:
        p = spf[n]
        pp = 1

        while n % p == 0:
            n //= p
            pp *= p

        res.append(pp)

    return res

# ---------------------------------------------------------
# Compute M(n)
# ---------------------------------------------------------
def M(n):
    if n == 1:
        return 0

    pp = prime_powers(n)
    k = len(pp)

    best = 1

    # Enumerate all subset products u
    for mask in range(1 << k):
        u = 1

        for i in range(k):
            if mask & (1 << i):
                u *= pp[i]

        v = n // u

        # trivial solutions 0 and 1
        if v == 1 or u == 1:
            continue

        # Solve:
        #   a ≡ 0 (mod u)
        #   a ≡ 1 (mod v)
        inv = invmod(u % v, v)

        a = u * inv

        if a > best:
            best = a

    return best

# ---------------------------------------------------------
# Main summation
# ---------------------------------------------------------
ans = 0

for n in range(1, LIMIT + 1):
    ans += M(n)

print(ans)

3. Code walkthrough

Sieve

spf = list(range(LIMIT + 1))

Initially every number is assumed prime.

The sieve fills spf[n] with the smallest prime dividing $n$.

This allows very fast factorizations.

`prime_powers(n)`

Example:

$$72 = 2^3\cdot 3^2$$

returns

[8, 9]

because the CRT works naturally with prime powers.

Enumerating CRT choices

Each prime power can belong either to:

the factor dividing $a$,
or the factor dividing $a-1$.

So we enumerate all subset products $u$.

Then:

$$v = n/u.$$

Since $u$ and $v$ are coprime, we solve

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

Writing $a=uk$:

$$uk \equiv 1 \pmod v$$

$$k = u^{-1}\pmod v.$$

Thus

a = u * invmod(u % v, v)

which is the CRT solution.

4. Small-example verification

For $n=6$:

Prime powers are:

$$6=2\cdot3.$$

Possible CRT choices:

$a \bmod 2$	$a \bmod 3$	solution
0	0	0
1	1	1
0	1	4
1	0	3

Largest is:

$$M(6)=4.$$

matching the statement.

5. Final verification

$M(p)=1$ for primes $p$, since only $0,1$ are idempotent modulo a prime.
Composite numbers with multiple distinct prime factors often have large nontrivial idempotents close to $n$.
The total sum should therefore be on the order of $10^7 \times 10^7 = 10^{14}$, which matches the magnitude of the computed result.

The computation yields:

$$\sum_{n=1}^{10^7} M(n) = 39782849136421.$$

Answer: 39782849136421

Project Euler Problem 407