Problem 451: Modular Inverses

Project Euler Problem 451 Solution

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

Problem Summary For each integer $n \ge 3$, define $$I(n)=\max\{m\in\mathbb{Z}_{>0}: m<n-1,\ m^{-1}\equiv m \pmod n\}.$$ The condition $m^{-1}\equiv m \pmod n$ means that $m$ is its own multiplicative inverse modulo $n$. Multiplying both sides by $m$ gives $$m^2\equiv 1 \pmod n.$$ So the task is to find, for every $n$, the largest solution of $x^2\equiv 1 \pmod n$ that is strictly smaller than $n-1$, and then sum those values for $3 \le n \le 2\cdot 10^7$. The checkpoints $I(7)=1$ and $I(100)=51$ confirm that we are looking for the largest nontrivial square root of $1$ modulo $n$, because $n-1\equiv -1\pmod n$ is always a solution but is excluded. Mathematical Approach Step 1: Self-inverse residues are roots of unity If $x^2\equiv 1 \pmod n$, then $n\mid (x-1)(x+1)$. Every solution is automatically coprime to $n$, because a common divisor of $x$ and $n$ would also divide $1$. Therefore the original inverse condition and the quadratic congruence are completely equivalent. This reduces the problem to understanding the solution set of $$x^2\equiv 1 \pmod n.$$ A brute-force scan over all $1\le x<n$ would be far too slow, so the implementation builds these roots from the prime-power factorization of $n$....

Detailed mathematical approach

Problem Summary

For each integer $n \ge 3$, define

$$I(n)=\max\{m\in\mathbb{Z}_{>0}: m<n-1,\ m^{-1}\equiv m \pmod n\}.$$

The condition $m^{-1}\equiv m \pmod n$ means that $m$ is its own multiplicative inverse modulo $n$. Multiplying both sides by $m$ gives

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

So the task is to find, for every $n$, the largest solution of $x^2\equiv 1 \pmod n$ that is strictly smaller than $n-1$, and then sum those values for $3 \le n \le 2\cdot 10^7$. The checkpoints $I(7)=1$ and $I(100)=51$ confirm that we are looking for the largest nontrivial square root of $1$ modulo $n$, because $n-1\equiv -1\pmod n$ is always a solution but is excluded.

Mathematical Approach

Step 1: Self-inverse residues are roots of unity

If $x^2\equiv 1 \pmod n$, then $n\mid (x-1)(x+1)$. Every solution is automatically coprime to $n$, because a common divisor of $x$ and $n$ would also divide $1$. Therefore the original inverse condition and the quadratic congruence are completely equivalent.

This reduces the problem to understanding the solution set of

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

A brute-force scan over all $1\le x<n$ would be far too slow, so the implementation builds these roots from the prime-power factorization of $n$.

Step 2: Solve the congruence on each prime power

Write the factorization of $n$ as

$$n=\prod_{j=1}^{t} q_j,$$

where each $q_j$ is a prime power and the factors are pairwise coprime.

For an odd prime power $p^a$, the congruence $x^2\equiv 1 \pmod{p^a}$ has exactly two solutions:

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

The reason is that $p^a\mid (x-1)(x+1)$, while $\gcd(x-1,x+1)\mid 2$. Since $p$ is odd, the two factors cannot both absorb a power of $p$, so one must contain all of $p^a$.

The power of $2$ is the only exceptional case:

$$x^2\equiv 1 \pmod{2^e}\quad\text{has}\quad \begin{cases} 1\text{ solution}, & e=1,\\ 2\text{ solutions}, & e=2,\\ 4\text{ solutions}, & e\ge 3. \end{cases}$$

For $e\ge 3$, the four residues are

$$x\equiv \pm 1,\qquad x\equiv 2^{e-1}\pm 1 \pmod{2^e}.$$

Indeed,

$$\left(2^{e-1}\pm 1\right)^2=2^{2e-2}\pm 2^e+1\equiv 1 \pmod{2^e}.$$

This special $2$-power behaviour is exactly why the implementation has two extra branches only when the factor $2^e$ satisfies $e\ge 3$.

Step 3: Recombine local roots with the Chinese Remainder Theorem

Because the prime-power factors $q_1,\dots,q_t$ are pairwise coprime, the Chinese Remainder Theorem says that a solution modulo $n$ is the same thing as choosing one local root modulo each $q_j$. Therefore the full solution set is the Cartesian product of the local solution sets.

If $r$ is the number of distinct odd prime divisors of $n$, then the total number of roots is

$$R(n)=2^r\cdot \begin{cases} 1, & 2\nmid n\text{ or }2\parallel n,\\ 2, & 4\mid n\text{ but }8\nmid n,\\ 4, & 8\mid n. \end{cases}$$

So the implementation never searches through all residues modulo $n$; it only enumerates these CRT combinations.

Step 4: The projector construction used by the implementation

For each prime-power factor $q_j$, define

$$M_j=\frac{n}{q_j}.$$

Since $\gcd(M_j,q_j)=1$, there exists an inverse $u_j$ such that

$$M_j u_j\equiv 1 \pmod{q_j}.$$

Now set

$$E_j\equiv M_j u_j \pmod n.$$

Then $E_j$ acts like a CRT projector:

$$E_j\equiv 1 \pmod{q_j},\qquad E_j\equiv 0 \pmod{q_i}\quad(i\ne j).$$

The implementations begin with the global residue

$$x_0=n-1,$$

which corresponds to choosing the local root $-1$ on every prime-power factor.

If a factor $q_j$ only has the two roots $\{-1,+1\}$, then switching the local choice from $-1$ to $+1$ changes that component by $2$. The lifted global correction is therefore

$$x\longmapsto x+2E_j \pmod n.$$

Nothing changes on the other factors, because $E_j\equiv 0$ there.

When $q_j=2^e$ with $e\ge 3$, the extra local roots are $2^{e-1}-1$ and $2^{e-1}+1$. Relative to the starting value $-1$, their offsets are

$$2^{e-1}\qquad\text{and}\qquad 2^{e-1}+2.$$

So the two additional lifted corrections are

$$x\longmapsto x+2^{e-1}E_j \pmod n,\qquad x\longmapsto x+\left(2^{e-1}+2\right)E_j \pmod n.$$

Processing the prime powers one by one and applying these corrections to every current candidate generates every root of $x^2\equiv 1 \pmod n$ exactly once.

Worked Example: $n=100$

This checkpoint illustrates the method cleanly. Since

$$100=4\cdot 25,$$

there are two prime-power factors.

For $q_1=4$, we have $M_1=25$, and $25\equiv 1 \pmod 4$, so

$$E_1\equiv 25 \pmod{100}.$$

For $q_2=25$, we have $M_2=4$. The inverse of $4$ modulo $25$ is $19$, hence

$$E_2\equiv 4\cdot 19=76 \pmod{100}.$$

Start from the trivial root

$$x_0=99.$$

Flipping the choice on the factor $4$ gives

$$99+2E_1\equiv 99+50\equiv 49 \pmod{100}.$$

Flipping the choice on the factor $25$ gives

$$99+2E_2\equiv 99+152\equiv 51 \pmod{100}.$$

Flipping both factors gives

$$99+2E_1+2E_2\equiv 1 \pmod{100}.$$

Thus the four roots are

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

The largest one below $99$ is $51$, so

$$I(100)=51.$$

How the Code Works

The C++, Python, and Java implementations first build a smallest-prime-factor sieve up to the overall limit. This allows each $n$ to be factored into prime powers quickly.

For a fixed $n$, the implementation starts from the residue $n-1$, factors $n$, and for each prime-power factor constructs the CRT projector described above. The modular inverse needed for that projector is computed with the extended Euclidean algorithm.

Each odd prime power, and also the factor $4$, contributes one extra branch obtained by replacing the local residue $-1$ with $+1$. A factor $2^e$ with $e\ge 3$ contributes two more branches coming from the additional roots $2^{e-1}\pm 1$. The candidate list is expanded incrementally, and the largest candidate strictly below $n-1$ is stored as $I(n)$.

Because the algorithm works with the full solution set of $x^2\equiv 1 \pmod n$, it is exact; there is no heuristic search and no dependence on trial residues.

Complexity Analysis

Let $N=2\cdot 10^7$. Building the smallest-prime-factor sieve costs $O(N\log\log N)$ time and $O(N)$ memory.

For one value of $n$, extracting the prime powers is close to logarithmic in practice, and the enumeration cost is proportional to the number of roots of $x^2\equiv 1 \pmod n$, not to $n$ itself. Since that root count is determined only by the distinct prime-power factors of $n$, it stays very small compared with $n$. This is what makes the full range up to $2\cdot 10^7$ feasible.

References

Problem page: https://projecteuler.net/problem=451
Chinese remainder theorem: Wikipedia — Chinese remainder theorem
Modular multiplicative inverse: Wikipedia — Modular multiplicative inverse
Kenneth H. Rosen, Elementary Number Theory and Its Applications, sections on linear congruences and the Chinese remainder theorem.

Problem 451 source code

C++

#include <pthread.h>
#include <algorithm>
#include <array>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <string>
#include <unistd.h>
#include <vector>

namespace {

using u32 = std::uint32_t;
using u64 = std::uint64_t;
using i64 = std::int64_t;

struct Options {
    int limit = 20'000'000;
    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;
    }
    int parsed = 0;
    for (char ch : tail) {
        if (ch < '0' || ch > '9') {
            return false;
        }
        parsed = parsed * 10 + static_cast<int>(ch - '0');
    }
    value = parsed;
    return true;
}

bool parse_arguments(int argc, char** argv, Options& options) {
    for (int i = 1; i < argc; ++i) {
        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 >= 3;
}

u32 mod_inverse(u32 a, u32 mod) {
    i64 x0 = 1;
    i64 y0 = 0;
    i64 x1 = 0;
    i64 y1 = 1;
    i64 aa = static_cast<i64>(a);
    i64 mm = static_cast<i64>(mod);

    while (mm != 0) {
        const i64 q = aa / mm;
        const i64 aa2 = aa - q * mm;
        aa = mm;
        mm = aa2;

        const i64 nx = x0 - q * x1;
        x0 = x1;
        x1 = nx;

        const i64 ny = y0 - q * y1;
        y0 = y1;
        y1 = ny;
    }

    if (aa != 1) {
        return 0U;
    }
    i64 v = x0 % static_cast<i64>(mod);
    if (v < 0) {
        v += static_cast<i64>(mod);
    }
    return static_cast<u32>(v);
}

std::vector<int> build_spf(const int limit) {
    std::vector<int> spf(static_cast<std::size_t>(limit + 1), 0);
    for (int i = 0; i <= limit; ++i) {
        spf[static_cast<std::size_t>(i)] = i;
    }
    for (int i = 2; i <= limit / i; ++i) {
        if (spf[static_cast<std::size_t>(i)] != i) {
            continue;
        }
        for (int j = i * i; j <= limit; j += i) {
            if (spf[static_cast<std::size_t>(j)] == j) {
                spf[static_cast<std::size_t>(j)] = i;
            }
        }
    }
    return spf;
}

u32 compute_I(const int n, const std::vector<int>& spf) {
    std::array<u32, 256> sums{};
    const u32 target = static_cast<u32>(n - 1);
    sums[0] = target;
    int count = 1;
    u32 best = 1U;

    auto push = [&](u32 s) {
        sums[static_cast<std::size_t>(count++)] = s;
        if (s > best && s < target) {
            best = s;
        }
    };

    int nn = n;
    while (nn != 1) {
        const int p = spf[static_cast<std::size_t>(nn)];
        u32 q = 1U;
        do {
            nn /= p;
            q *= static_cast<u32>(p);
        } while (nn % p == 0);

        if (q == 2U) {
            continue;
        }

        const u32 m = static_cast<u32>(n / static_cast<int>(q));
        const u32 inv = mod_inverse(m % q, q);
        const u32 b = static_cast<u32>((static_cast<u64>(m) * inv) % static_cast<u64>(n));
        const u32 delta2 = static_cast<u32>((2ULL * b) % static_cast<u64>(n));
        u32 delta_half = 0U;
        u32 delta_half_plus_two = 0U;
        if ((q & 1U) == 0U && q >= 8U) {
            const u32 half = q / 2U;
            delta_half = static_cast<u32>((static_cast<u64>(half) * b) % static_cast<u64>(n));
            delta_half_plus_two =
                static_cast<u32>((static_cast<u64>(half + 2U) * b) % static_cast<u64>(n));
        }

        const int current_count = count;
        for (int j = 0; j < current_count; ++j) {
            const u32 base = sums[static_cast<std::size_t>(j)];

            u32 s = base + delta2;
            if (s >= static_cast<u32>(n)) {
                s -= static_cast<u32>(n);
            }
            push(s);

            if ((q & 1U) == 0U && q >= 8U) {
                u32 s1 = base + delta_half;
                if (s1 >= static_cast<u32>(n)) {
                    s1 -= static_cast<u32>(n);
                }
                push(s1);

                u32 s2 = base + delta_half_plus_two;
                if (s2 >= static_cast<u32>(n)) {
                    s2 -= static_cast<u32>(n);
                }
                push(s2);
            }
        }
    }
    return best;
}

u64 solve_range(const std::vector<int>& spf, const int begin, const int end) {
    u64 sum = 0ULL;
    for (int n = begin; n <= end; ++n) {
        sum += compute_I(n, spf);
    }
    return sum;
}

struct WorkerTask {
    const std::vector<int>* spf = nullptr;
    int begin = 0;
    int end = -1;
    u64 partial = 0ULL;
};

void* worker_entry(void* raw) {
    auto* task = static_cast<WorkerTask*>(raw);
    task->partial = solve_range(*task->spf, task->begin, task->end);
    return nullptr;
}

int detect_thread_count(const int work_items) {
    long cores = ::sysconf(_SC_NPROCESSORS_ONLN);
    int threads = (cores > 0) ? static_cast<int>(cores) : 4;
    if (threads < 1) {
        threads = 1;
    }
    if (threads > work_items) {
        threads = work_items;
    }
    return threads;
}

u64 solve(const int limit) {
    const std::vector<int> spf = build_spf(limit);
    const int begin = 3;
    if (limit < begin) {
        return 0ULL;
    }
    const int work_items = limit - begin + 1;
    const int threads = detect_thread_count(work_items);
    if (threads == 1) {
        return solve_range(spf, begin, limit);
    }

    std::vector<pthread_t> handles(static_cast<std::size_t>(threads));
    std::vector<WorkerTask> tasks(static_cast<std::size_t>(threads));
    std::vector<char> launched(static_cast<std::size_t>(threads), 0);

    const int chunk = (work_items + threads - 1) / threads;
    int start = begin;
    for (int t = 0; t < threads; ++t) {
        if (start > limit) {
            tasks[static_cast<std::size_t>(t)].partial = 0ULL;
            continue;
        }
        const int end = std::min(limit, start + chunk - 1);
        auto& task = tasks[static_cast<std::size_t>(t)];
        task.spf = &spf;
        task.begin = start;
        task.end = end;
        task.partial = 0ULL;

        const int rc = pthread_create(&handles[static_cast<std::size_t>(t)], nullptr, worker_entry, &task);
        if (rc == 0) {
            launched[static_cast<std::size_t>(t)] = 1;
        } else {
            task.partial = solve_range(spf, start, end);
        }
        start = end + 1;
    }

    u64 sum = 0ULL;
    for (int t = 0; t < threads; ++t) {
        if (launched[static_cast<std::size_t>(t)] != 0) {
            pthread_join(handles[static_cast<std::size_t>(t)], nullptr);
        }
        sum += tasks[static_cast<std::size_t>(t)].partial;
    }
    return sum;
}

u32 brute_I(const int n) {
    u32 best = 1U;
    for (int m = 1; m < n - 1; ++m) {
        if (std::gcd(m, n) != 1) {
            continue;
        }
        if ((static_cast<u64>(m) * static_cast<u64>(m)) % static_cast<u64>(n) == 1ULL) {
            best = static_cast<u32>(m);
        }
    }
    return best;
}

bool run_checkpoints() {
    const std::vector<int> spf = build_spf(2'000);
    if (compute_I(7, spf) != 1U) {
        std::cerr << "Checkpoint failed: I(7)=1" << '\n';
        return false;
    }
    if (compute_I(100, spf) != 51U) {
        std::cerr << "Checkpoint failed: I(100)=51" << '\n';
        return false;
    }
    for (int n = 3; n <= 600; ++n) {
        if (compute_I(n, spf) != brute_I(n)) {
            std::cerr << "Checkpoint failed: brute-force cross-check at n=" << n << '\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 << solve(options.limit) << '\n';
    return 0;
}

Python

def solve():
    LIMIT = 20_000_000

    # Build SPF sieve
    spf = list(range(LIMIT + 1))
    for i in range(2, int(LIMIT**0.5) + 1):
        if spf[i] == i:
            for j in range(i*i, LIMIT+1, i):
                if spf[j] == j:
                    spf[j] = i

    def mod_inv(a, m):
        b, x0, x1 = m, 1, 0
        while b:
            q = a // b
            a, b = b, a - q * b
            x0, x1 = x1, x0 - q * x1
        return x0 % m

    def compute_I(n):
        target = n - 1
        sums = [target]
        best = 1

        nn = n
        while nn != 1:
            p = spf[nn]
            q = 1
            while nn % p == 0:
                nn //= p
                q *= p
            if q == 2:
                continue

            m = n // q
            inv = mod_inv(m % q, q)
            b = (m * inv) % n
            delta2 = (2 * b) % n

            delta_half = 0
            delta_half_plus_two = 0
            if q % 2 == 0 and q >= 8:
                half = q // 2
                delta_half = (half * b) % n
                delta_half_plus_two = ((half + 2) * b) % n

            current = list(sums)
            for base in current:
                s = (base + delta2) % n
                sums.append(s)
                if s > best and s < target: best = s

                if q % 2 == 0 and q >= 8:
                    s1 = (base + delta_half) % n
                    sums.append(s1)
                    if s1 > best and s1 < target: best = s1
                    s2 = (base + delta_half_plus_two) % n
                    sums.append(s2)
                    if s2 > best and s2 < target: best = s2

        return best

    total = sum(compute_I(n) for n in range(3, LIMIT + 1))
    return str(total)

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

Java

import java.util.stream.IntStream;

public class Euler451 {
    static final int LIMIT = 20000000;

    static int modInverse(int a, int mod) {
        long x0 = 1, y0 = 0;
        long x1 = 0, y1 = 1;
        long aa = a, mm = mod;

        while (mm != 0) {
            long q = aa / mm;
            long aa2 = aa - q * mm;
            aa = mm;
            mm = aa2;

            long nx = x0 - q * x1;
            x0 = x1;
            x1 = nx;
        }

        if (aa != 1)
            return 0;
        long v = x0 % mod;
        if (v < 0)
            v += mod;
        return (int) v;
    }

    static int[] buildSpf(int limit) {
        int[] spf = new int[limit + 1];
        for (int i = 0; i <= limit; i++)
            spf[i] = i;
        for (int i = 2; i * i <= limit; i++) {
            if (spf[i] == i) {
                for (int j = i * i; j <= limit; j += i) {
                    if (spf[j] == j) {
                        spf[j] = i;
                    }
                }
            }
        }
        return spf;
    }

    static int computeI(int n, int[] spf) {
        int[] sums = new int[256];
        int target = n - 1;
        sums[0] = target;
        int count = 1;
        int best = 1;

        int nn = n;
        while (nn != 1) {
            int p = spf[nn];
            long q = 1;
            do {
                nn /= p;
                q *= p;
            } while (nn % p == 0);

            if (q == 2)
                continue;

            long m = n / q;
            long inv = modInverse((int) (m % q), (int) q);
            long b = (m * inv) % n;
            long delta2 = (2L * b) % n;
            long deltaHalf = 0;
            long deltaHalfPlusTwo = 0;

            if ((q & 1) == 0 && q >= 8) {
                long half = q / 2;
                deltaHalf = (half * b) % n;
                deltaHalfPlusTwo = ((half + 2) * b) % n;
            }

            int currentCount = count;
            for (int j = 0; j < currentCount; j++) {
                int base = sums[j];

                long sLong = base + delta2;
                if (sLong >= n)
                    sLong -= n;
                int s = (int) sLong;
                sums[count++] = s;
                if (s > best && s < target)
                    best = s;

                if ((q & 1) == 0 && q >= 8) {
                    long s1Long = base + deltaHalf;
                    if (s1Long >= n)
                        s1Long -= n;
                    int s1 = (int) s1Long;
                    sums[count++] = s1;
                    if (s1 > best && s1 < target)
                        best = s1;

                    long s2Long = base + deltaHalfPlusTwo;
                    if (s2Long >= n)
                        s2Long -= n;
                    int s2 = (int) s2Long;
                    sums[count++] = s2;
                    if (s2 > best && s2 < target)
                        best = s2;
                }
            }
        }
        return best;
    }

    public static String solve() {
        int[] spf = buildSpf(LIMIT);
        int chunkSize = 16384;
        int totalChunks = (LIMIT - 3 + 1 + chunkSize - 1) / chunkSize;

        long totalSum = IntStream.range(0, totalChunks).parallel().mapToLong(chunkIdx -> {
            int start = 3 + chunkIdx * chunkSize;
            int end = Math.min(LIMIT, start + chunkSize - 1);
            long localSum = 0;
            for (int n = start; n <= end; n++) {
                localSum += computeI(n, spf);
            }
            return localSum;
        }).sum();

        return Long.toString(totalSum);
    }

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