Project Euler Problem 451

Solution

Answer: 153651073760956

We want all integers $m$ such that

$$m^2 \equiv 1 \pmod n,$$

because an element is its own modular inverse iff $m^{-1}\equiv m\pmod n$.

The problem defines $I(n)$ as the largest such $m<n-1$.

---

1. Mathematical analysis

We must solve

$$x^2 \equiv 1 \pmod n.$$

Equivalently,

$$(x-1)(x+1)\equiv 0 \pmod n.$$

The key structure comes from the Chinese Remainder Theorem.

---

Step 1: Solutions modulo prime powers

For an odd prime power $p^k$:

$$x^2\equiv 1 \pmod{p^k}$$

has exactly two solutions:

$$x\equiv \pm 1 \pmod{p^k}.$$

For powers of $2$:

• modulo $2$: one solution
• modulo $4$: two solutions
• modulo $2^k$ for $k\ge 3$: four solutions

$$x\equiv \pm1,\quad \pm(1+2^{k-1}) \pmod{2^k}.$$

---

Step 2: CRT decomposition

If

$$n=\prod p_i^{a_i},$$

then every solution modulo $n$ is obtained by independently choosing one solution modulo each prime power and combining them with CRT.

Thus the nontrivial solutions arise by assigning some prime-power factors the residue $+1$ and others the residue $-1$.

---

Step 3: Characterizing the largest solution

Suppose

$$n=ab,\qquad \gcd(a,b)=1.$$

Take the CRT system

$$x\equiv -1 \pmod a,\qquad x\equiv 1 \pmod b.$$

CRT gives a unique solution modulo $n$.

Every nontrivial involution comes from such a split.

If $x$ is a solution, then so is $n-x$.

Therefore the largest solution below $n-1$ is obtained from the smallest nontrivial solution.

So define

$$s(n)=\min{x>1 : x^2\equiv1\pmod n},$$

then

$$I(n)=n-s(n).$$

---

Step 4: Efficient computation

A direct search is impossible up to $2\times10^7$.

Instead:

1. Sieve smallest prime factors.
2. Factor each $n$.
3. Generate all CRT involutions from the prime-power decomposition.
4. Keep the smallest nontrivial solution $s(n)$.
5. Add $n-s(n)$.

The number of involutions is tiny because the number of distinct prime factors is tiny.

The total runtime is easily feasible in optimized Python/PyPy.

---

2. Python implementation

from math import prod

LIMIT = 20_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 Euclid modular inverse
# ------------------------------------------------------------

def inv_mod(a, 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

# ------------------------------------------------------------
# CRT merge:
# x ≡ r1 (mod m1)
# x ≡ r2 (mod m2)
# gcd(m1,m2)=1
# ------------------------------------------------------------

def crt(r1, m1, r2, m2):
    t = ((r2 - r1) * inv_mod(m1, m2)) % m2
    return r1 + m1 * t

# ------------------------------------------------------------
# Generate prime powers
# ------------------------------------------------------------

def factorize(n):
    out = []
    while n > 1:
        p = spf[n]
        e = 0
        pp = 1
        while n % p == 0:
            n //= p
            e += 1
            pp *= p
        out.append((p, e, pp))
    return out

# ------------------------------------------------------------
# Solutions modulo p^k
# ------------------------------------------------------------

def local_roots(p, e, pe):
    if p == 2:
        if e == 1:
            return [1]
        elif e == 2:
            return [1, 3]
        else:
            h = 1 << (e - 1)
            return [1, pe - 1, 1 + h, h - 1]
    else:
        return [1, pe - 1]

# ------------------------------------------------------------
# Main
# ------------------------------------------------------------

answer = 0

for n in range(3, LIMIT + 1):

    fac = factorize(n)

    roots = [(0, 1)]

    # build all CRT combinations
    for p, e, pe in fac:
        local = local_roots(p, e, pe)

        new_roots = []

        for r, mod in roots:
            for t in local:
                nr = crt(r, mod, t, pe)
                new_roots.append((nr, mod * pe))

        roots = new_roots

    best = 1

    for r, _ in roots:
        r %= n
        if r != n - 1:
            best = max(best, r)

    answer += best

print(answer)

---

3. Code walkthrough

SPF sieve

We precompute the smallest prime factor of every number up to $2\times10^7$.

This allows factorization in nearly constant amortized time.

---

factorize(n)

Produces

$$n=\prod p_i^{e_i}$$

together with the actual prime powers $p_i^{e_i}$.

---

local_roots

Returns all solutions to

$$x^2\equiv1\pmod{p^e}.$$

For odd primes: only $\pm1$.

For powers of two ($e\ge3$):

$$1,\ -1,\ 1+2^{e-1},\ -1+2^{e-1}.$$

---

CRT combination

We recursively combine local solutions into global solutions modulo $n$.

Because the number of distinct prime factors is tiny, the total number of roots is very small.

---

Extracting $I(n)$

Among all roots, $n-1$ is always present (the trivial solution $-1$).

We choose the largest root strictly smaller than $n-1$.

---

4. Verification on small examples

$n=15$

Roots are

$$1,4,11,14.$$

Largest below $14$ is

$$I(15)=11.$$

Correct.

---

$n=100$

Roots are

$$1, 49, 51, 99.$$

Largest below $99$ is

$$I(100)=51.$$

Correct.

---

$n=7$

Only roots:

$$1,6.$$

Since $6=n-1$ is excluded:

$$I(7)=1.$$

Correct.

---

5. Final verification

The sum is on the order of

$$\sum_{n\le 2\times10^7} O(n),$$

so a result around $10^{14}$–$10^{15}$ is reasonable.

The computation is entirely deterministic and based on the full classification of involutions modulo $n$.

---

Answer: 153651073760956