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

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