Project Euler

Modular Inverses

For a positive integer n, define I(n) as the largest positive integer m < n - 1 such that m^2 equiv 1 (mod n). If no such m exists, set I(n) = 1. Known values: I(7) = 1, I(15) = 11, I(100) = 51. Co...

Source sync Apr 19, 2026

Problem #0451

Level Level 11

Solved By 1,718

Languages C++, Python

Answer 153651073760956

Length 333 words

modular_arithmeticnumber_theorylinear_algebra

Consider the number .
There are eight positive numbers less than which are coprime to : .
The modular inverses of these numbers modulo are:
because

Let be the largest positive number smaller than such that the modular inverse of modulo equals itself.
So .
Also and .

Find for .

Problem 451: Modular Inverses

Mathematical Foundation

Definition. A number $m$ with $1 \le m < n$ and $\gcd(m, n) = 1$ is a self-inverse modulo $n$ if $m^2 \equiv 1 \pmod{n}$ .

Theorem 1 (CRT factorization of self-inverses). Let $n = p_1^{a_1} p_2^{a_2} \cdots p_k^{a_k}$ . Then

m^2 \equiv 1 \pmod{n} \iff m^2 \equiv 1 \pmod{p_i^{a_i}} \text{ for all } i = 1, \ldots, k.

Proof. This is an immediate consequence of the Chinese Remainder Theorem: $\mathbb{Z}/n\mathbb{Z} \cong \prod_{i=1}^{k} \mathbb{Z}/p_i^{a_i}\mathbb{Z}$ , and the equation $m^2 - 1 \equiv 0$ holds in the product ring if and only if it holds in each factor. $\square$

Lemma 1 (Solutions modulo prime powers). The solutions of $x^2 \equiv 1 \pmod{p^a}$ are:

For $p$ odd: exactly 2 solutions, $x \equiv \pm 1 \pmod{p^a}$ .
For $p = 2$ , $a = 1$ : 1 solution, $x = 1$ .
For $p = 2$ , $a = 2$ : 2 solutions, $x \in \{1, 3\}$ .
For $p = 2$ , $a \ge 3$ : 4 solutions, $x \in \{1,\; 2^{a-1}-1,\; 2^{a-1}+1,\; 2^a-1\}$ .

Proof. For odd $p$ , the group $(\mathbb{Z}/p^a\mathbb{Z})^*$ is cyclic of order $\varphi(p^a) = p^{a-1}(p-1)$ . In a cyclic group, the equation $x^2 = e$ has at most 2 solutions. Since $x = \pm 1$ are always solutions, there are exactly 2.

For $p = 2$ and $a \ge 3$ , we have $(\mathbb{Z}/2^a\mathbb{Z})^* \cong \mathbb{Z}/2 \times \mathbb{Z}/2^{a-2}$ . The equation $x^2 = 1$ in a product of cyclic groups $\mathbb{Z}/2 \times \mathbb{Z}/2^{a-2}$ has $2 \times 2 = 4$ solutions (the identity and the unique element of order 2 in each factor). Direct verification:

1^2 \equiv 1, \quad (2^{a-1}-1)^2 = 2^{2a-2} - 2^a + 1 \equiv 1, \quad (2^{a-1}+1)^2 = 2^{2a-2} + 2^a + 1 \equiv 1, \quad (2^a - 1)^2 = 2^{2a} - 2^{a+1} + 1 \equiv 1 \pmod{2^a}.

The cases $a = 1$ and $a = 2$ are verified by exhaustion. $\square$

Corollary. The total number of self-inverses modulo $n = \prod p_i^{a_i}$ is $\prod_{i=1}^{k} s(p_i^{a_i})$ , where $s(p^a)$ denotes the count from Lemma 1.

Theorem 2 (Characterization of $I(n)$ ). Given the complete set of self-inverses modulo $n$ (constructed via CRT from Theorem 1 and Lemma 1), $I(n)$ equals the second-largest element of this set. The largest is always $n - 1$ ; if the only solutions are $\{1, n-1\}$ , then $I(n) = 1$ .

Proof. Note that $m = n - 1$ satisfies $(n-1)^2 = n^2 - 2n + 1 \equiv 1 \pmod{n}$ , so $n-1$ is always a self-inverse. By definition, $I(n)$ is the largest self-inverse strictly less than $n - 1$ . If no non-trivial self-inverse exists, the only solution below $n-1$ is $m = 1$ , giving $I(n) = 1$ . $\square$

Editorial

Approach: Sieve-based CRT. For each n, we build up the set of self-inverses by processing one prime-power factor at a time. We store only the set of known self-inverses so far (from partial factorization). We sieve smallest prime factor spf[i] for 2 <= i <= N. Finally, return total.

Pseudocode

Sieve smallest prime factor spf[i] for 2 <= i <= N
total <- 0
For n = 3 to N:
Return total

Complexity Analysis

Time: $O(N \log \log N)$ for the SPF sieve. For each $n$ , factoring via SPF takes $O(\log n)$ . The CRT reconstruction processes $\omega(n)$ prime factors with up to $2^{\omega(n)}$ solutions; since $\omega(n) \le O(\log n / \log \log n)$ and on average $\omega(n) \approx \log \log n$ , the amortized cost per $n$ is small. Total: $O(N \log \log N)$ .
Space: $O(N)$ for the SPF array.

Answer

\boxed{153651073760956}

C++ project_euler/problem_451/solution.cpp

#include <bits/stdc++.h>
using namespace std;

/*
 * Project Euler Problem 451: Modular Inverses
 *
 * Find sum of I(n) for 3 <= n <= 2*10^7, where I(n) is the largest m < n-1
 * such that m^2 = 1 (mod n).
 *
 * Approach: Sieve-based. For each n, we store a vector of all self-inverse
 * residues. We build these up by processing prime powers: for each prime power
 * q = p^a, for each multiple n of q, we combine existing solutions with
 * solutions mod q via CRT.
 *
 * Answer: 153651073760956
 */

const int MAXN = 20000001;

// spf[i] = smallest prime factor of i
int spf[MAXN];

// For each n, store all self-inverse solutions found so far
// We use a vector of vectors; memory-intensive but manageable for 2*10^7
// Actually, we'll use a more memory-efficient approach:
// Process prime factors one at a time via SPF, combining CRT incrementally.

long long extgcd(long long a, long long b, long long &x, long long &y) {
    if (a == 0) { x = 0; y = 1; return b; }
    long long x1, y1;
    long long g = extgcd(b % a, a, x1, y1);
    x = y1 - (b / a) * x1;
    y = x1;
    return g;
}

// CRT: x = r1 mod m1, x = r2 mod m2, gcd(m1,m2)=1
// Returns x mod m1*m2
long long crt2(long long r1, long long m1, long long r2, long long m2) {
    long long x, y;
    extgcd(m1, m2, x, y);
    long long mod = m1 * m2;
    long long res = (r1 + m1 * ((r2 - r1) % m2 + m2) % m2 * x) % mod;
    // More carefully:
    long long diff = ((r2 - r1) % m2 + m2) % m2;
    res = r1 + m1 * (diff * ((x % m2 + m2) % m2) % m2);
    res = ((res % mod) + mod) % mod;
    return res;
}

int main() {
    ios_base::sync_with_stdio(false);
    cin.tie(NULL);

    const int N = 20000000;

    // Sieve smallest prime factor
    for (int i = 0; i <= N; i++) spf[i] = i;
    for (int i = 2; (long long)i * i <= N; i++) {
        if (spf[i] == i) { // prime
            for (int j = i * i; j <= N; j += i) {
                if (spf[j] == j) spf[j] = i;
            }
        }
    }

    // For each n, factor using spf and compute all self-inverses via CRT
    // Then I(n) = second largest (largest < n-1)

    long long total = 0;

    for (int n = 3; n <= N; n++) {
        // Factor n
        int temp = n;
        // Collect prime powers
        // Use small static array since n <= 2*10^7 has at most ~8 prime factors
        int pp[10]; // prime powers p^a
        int npp = 0;

        while (temp > 1) {
            int p = spf[temp];
            int pa = 1;
            while (temp % p == 0) {
                pa *= p;
                temp /= p;
            }
            pp[npp++] = pa;
            // Store p as well for the 2-special case
            // We need to know if p==2 and what a is
            // Let's store both
            // Redo: store (p, pa) pairs
            pp[npp - 1] = pa; // we'll handle p==2 by checking pa
            // Actually let's use a struct approach
            // Simpler: just redo with separate arrays
            (void)p; // we'll fix below
        }

        // Redo factoring properly
        temp = n;
        npp = 0;
        int primes[10], pas[10];
        while (temp > 1) {
            int p = spf[temp];
            int pa = 1;
            while (temp % p == 0) {
                pa *= p;
                temp /= p;
            }
            primes[npp] = p;
            pas[npp] = pa;
            npp++;
        }

        // Build all self-inverse solutions via CRT
        // Start with solutions = {1} mod 1, current_mod = 1
        // Max solutions: 2^(npp) * possible factor of 2 for p=2, a>=3
        // Worst case about 2^9 = 512 solutions, but typically much fewer
        static long long sols[1024], new_sols[1024];
        int nsols = 1;
        sols[0] = 1;
        long long cur_mod = 1;

        for (int i = 0; i < npp; i++) {
            int p = primes[i];
            long long pa = pas[i];

            // Local solutions mod pa
            long long local_s[4];
            int nloc;
            if (p == 2) {
                if (pa == 1) {
                    local_s[0] = 1;
                    nloc = 1;
                } else if (pa == 2) {
                    // Actually for mod 4: 1^2=1, 3^2=9=1 mod 4. So {1, 3}.
                    local_s[0] = 1; local_s[1] = 3;
                    nloc = 2;
                } else {
                    // pa >= 8
                    long long half = pa / 2;
                    local_s[0] = 1;
                    local_s[1] = half - 1;
                    local_s[2] = half + 1;
                    local_s[3] = pa - 1;
                    nloc = 4;
                }
            } else {
                local_s[0] = 1;
                local_s[1] = pa - 1;
                nloc = 2;
            }

            // Combine via CRT
            int nnew = 0;
            for (int a = 0; a < nsols; a++) {
                for (int b = 0; b < nloc; b++) {
                    new_sols[nnew++] = crt2(sols[a], cur_mod, local_s[b], pa);
                }
            }

            cur_mod *= pa;
            nsols = nnew;
            memcpy(sols, new_sols, nsols * sizeof(long long));
        }

        // Find largest solution < n-1
        long long best = 1;
        for (int i = 0; i < nsols; i++) {
            long long s = sols[i];
            if (s > 1 && s < n - 1 && s > best) {
                best = s;
            }
        }

        total += best;
    }

    cout << total << endl;

    // Verification
    if (total == 153651073760956LL) {
        cout << "CORRECT!" << endl;
    } else {
        cout << "WRONG! Expected 153651073760956" << endl;
    }

    return 0;
}

Python project_euler/problem_451/solution.py

"""
Project Euler Problem 451: Modular Inverses

Find sum of I(n) for 3 <= n <= 2*10^7, where I(n) is the largest m < n-1
such that m^2 = 1 (mod n).

Approach: Sieve-based CRT. For each n, we build up the set of self-inverses
by processing one prime-power factor at a time. We store only the set of
known self-inverses so far (from partial factorization).

Answer: 153651073760956
"""

import time
from math import gcd

def compute_I_single(n):
    """Brute-force I(n) for verification: largest m < n-1 with m^2 = 1 mod n."""
    for m in range(n - 2, 0, -1):
        if (m * m) % n == 1:
            return m
    return 0  # only for n <= 2

def solve_sieve(N):
    """
    Compute I(n) for all 3 <= n <= N using a sieve approach.

    For each n, we maintain a list of all self-inverse residues found so far
    (from partial factorization). We process prime powers via a smallest-prime-
    factor sieve, and for each prime power q dividing n, we combine existing
    solutions with the solutions mod q using CRT.

    Returns array result where result[n] = I(n).
    """
    # Step 1: Smallest prime factor sieve
    spf = list(range(N + 1))  # spf[i] = smallest prime factor of i
    for i in range(2, int(N**0.5) + 1):
        if spf[i] == i:  # i is prime
            for j in range(i * i, N + 1, i):
                if spf[j] == j:
                    spf[j] = i

    # Step 2: For each n, compute I(n) by factoring via spf and using CRT
    result = [0] * (N + 1)

    for n in range(3, N + 1):
        # Factor n using spf
        temp = n
        prime_powers = []  # list of (p, p^a)
        while temp > 1:
            p = spf[temp]
            pa = 1
            while temp % p == 0:
                pa *= p
                temp //= p
            prime_powers.append((p, pa))

        # For each prime power, find solutions to x^2 = 1 mod p^a
        # Then combine via CRT
        # Solutions mod p^a for odd p: {1, p^a - 1}
        # Solutions mod 2^a: {1} if a=1, {1,3} if a=2, {1, 2^(a-1)-1, 2^(a-1)+1, 2^a-1} if a>=3

        all_solutions = [1]  # start with just x = 1 mod 1

        current_mod = 1
        for p, pa in prime_powers:
            if p == 2:
                if pa == 1:
                    local_sols = [1]
                elif pa == 2:
                    local_sols = [1, 3]
                else:
                    half = pa // 2
                    local_sols = [1, half - 1, half + 1, pa - 1]
            else:
                local_sols = [1, pa - 1]

            # Combine all_solutions (mod current_mod) with local_sols (mod pa) via CRT
            new_solutions = []
            for s1 in all_solutions:
                for s2 in local_sols:
                    # CRT: find x = s1 mod current_mod, x = s2 mod pa
                    combined = crt_two(s1, current_mod, s2, pa)
                    if combined is not None:
                        new_solutions.append(combined)

            current_mod *= pa
            all_solutions = new_solutions

        # I(n) = largest solution that is < n-1 and > 1
        # Solutions include 1 and n-1 always
        best = 1
        for s in all_solutions:
            if 1 < s < n - 1 and s > best:
                best = s

        # If no non-trivial solution exists, I(n) = 1
        result[n] = best

    return result

def extended_gcd(a, b):
    """Extended Euclidean algorithm. Returns (g, x, y) where a*x + b*y = g."""
    if a == 0:
        return b, 0, 1
    g, x, y = extended_gcd(b % a, a)
    return g, y - (b // a) * x, x

def crt_two(r1, m1, r2, m2):
    """
    Chinese Remainder Theorem for two congruences:
    x = r1 mod m1, x = r2 mod m2.
    Returns x mod (m1*m2), or None if no solution.
    Assumes gcd(m1, m2) = 1 (which holds since m1, m2 are coprime prime powers).
    """
    g, p, q = extended_gcd(m1, m2)
    if (r2 - r1) % g != 0:
        return None
    lcm = m1 * m2 // g
    x = (r1 + m1 * ((r2 - r1) // g) * p) % lcm
    return x

def solve_brute(N):
    """Brute force for small N, for verification."""
    result = [0] * (N + 1)
    for n in range(3, N + 1):
        result[n] = compute_I_single(n)
    return result

def create_visualization(N_viz=1000):
    """Create visualization of I(n) values and save to PNG."""
    result = solve_sieve(N_viz)

    ns = list(range(3, N_viz + 1))
    Is = [result[n] for n in ns]