Problem 407: Idempotents
View on Project EulerProject Euler Problem 407 Solution
EulerSolve provides an optimized solution for Project Euler Problem 407, Idempotents, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary For each integer \(n \le N\), define \(M(n)\) as the largest residue \(a \lt n\) satisfying $$a^2 \equiv a \pmod{n}.$$ The goal is to compute $$\sum_{n=1}^{N} M(n)$$ for the large bound \(N=10^7\). The key observation is that these residues are exactly the idempotents modulo \(n\), and their structure is controlled by the prime-power factorization of \(n\). Mathematical Approach Step 1: Rewrite the Congruence The defining condition is equivalent to $$a^2 \equiv a \pmod{n}\iff a(a-1)\equiv 0 \pmod{n}.$$ This form is more informative because the two factors differ by \(1\), so $$\gcd(a,a-1)=1.$$ Therefore every prime power dividing \(n\) must divide one of these two coprime factors completely. Step 2: Prime-Power Blocks Force Local Choices Write the factorization of \(n\) as $$n=\prod_{i=1}^{k} q_i,\qquad q_i=p_i^{e_i},$$ where the numbers \(q_1,\dots,q_k\) are pairwise coprime prime powers. If \(q_i\mid a(a-1)\), then the prime underlying \(q_i\) cannot divide both \(a\) and \(a-1\), so the whole block \(q_i\) must divide exactly one of them....
Detailed mathematical approach
Problem Summary
For each integer \(n \le N\), define \(M(n)\) as the largest residue \(a \lt n\) satisfying
$$a^2 \equiv a \pmod{n}.$$
The goal is to compute
$$\sum_{n=1}^{N} M(n)$$
for the large bound \(N=10^7\). The key observation is that these residues are exactly the idempotents modulo \(n\), and their structure is controlled by the prime-power factorization of \(n\).
Mathematical Approach
Step 1: Rewrite the Congruence
The defining condition is equivalent to
$$a^2 \equiv a \pmod{n}\iff a(a-1)\equiv 0 \pmod{n}.$$
This form is more informative because the two factors differ by \(1\), so
$$\gcd(a,a-1)=1.$$
Therefore every prime power dividing \(n\) must divide one of these two coprime factors completely.
Step 2: Prime-Power Blocks Force Local Choices
Write the factorization of \(n\) as
$$n=\prod_{i=1}^{k} q_i,\qquad q_i=p_i^{e_i},$$
where the numbers \(q_1,\dots,q_k\) are pairwise coprime prime powers. If \(q_i\mid a(a-1)\), then the prime underlying \(q_i\) cannot divide both \(a\) and \(a-1\), so the whole block \(q_i\) must divide exactly one of them. Hence for each \(i\),
$$a\equiv 0 \pmod{q_i}\qquad\text{or}\qquad a\equiv 1 \pmod{q_i}.$$
Conversely, any residue that is \(0\) or \(1\) modulo every \(q_i\) is automatically idempotent modulo \(n\), because locally \(0^2=0\) and \(1^2=1\), and the Chinese remainder theorem glues the local conditions into one residue modulo \(n\).
So idempotents modulo \(n\) are in one-to-one correspondence with the \(2^k\) ways to choose, for each prime-power block, whether the residue is \(0\) or \(1\).
Step 3: Encode an Idempotent by a Subset
Choose a subset \(U\subseteq\{1,\dots,k\}\). Let
$$u=\prod_{i\in U} q_i,\qquad v=\frac{n}{u}.$$
The subset \(U\) means: force the residue to be \(0\) on the blocks inside \(U\), and \(1\) on the complementary blocks. That is, solve
$$a\equiv 0 \pmod{u},\qquad a\equiv 1 \pmod{v}.$$
Because \(u\) and \(v\) are coprime, the Chinese remainder theorem gives a unique solution modulo \(n\).
Step 4: Closed Form for Each Candidate
Since \(a\equiv 0\pmod{u}\), write \(a=ut\). The second congruence becomes
$$ut\equiv 1 \pmod{v}.$$
As \(\gcd(u,v)=1\), a unique inverse exists modulo \(v\). Let \(t_U\) be the integer with
$$0\le t_U\lt v,\qquad ut_U\equiv 1 \pmod{v}.$$
Then the idempotent attached to \(U\) is
$$a_U = u\,t_U.$$
When \(U=\emptyset\), we get \(u=1\), \(a=1\). When \(U=\{1,\dots,k\}\), we get \(u=n\), which corresponds to \(a=0\). These two trivial idempotents can never beat a nontrivial larger residue, so the maximum only needs nonempty proper subsets. Therefore
$$M(n)=\begin{cases} 0, & n=1,\\ 1, & n>1 \text{ and } k=1,\\ \max\limits_{\emptyset\neq U\neq\{1,\dots,k\}} a_U, & k\ge 2. \end{cases}$$
The case \(k=1\) means \(n\) is a prime power. Then only \(0\) and \(1\) are idempotent, so the largest residue below \(n\) is \(1\).
Step 5: Worked Example
Take \(n=12\). Its prime-power blocks are
$$12=4\cdot 3.$$
There are four idempotent patterns:
$$\begin{aligned} (0 \bmod 4,\ 0 \bmod 3) &\Rightarrow a=0,\\ (1 \bmod 4,\ 1 \bmod 3) &\Rightarrow a=1,\\ (0 \bmod 4,\ 1 \bmod 3) &\Rightarrow a=4,\\ (1 \bmod 4,\ 0 \bmod 3) &\Rightarrow a=9. \end{aligned}$$
So the idempotents modulo \(12\) are \(\{0,1,4,9\}\), and therefore
$$M(12)=9.$$
The same construction for \(n=6=2\cdot 3\) gives the nontrivial idempotents \(3\) and \(4\), hence \(M(6)=4\), matching the standard checkpoint for this problem.
How the Code Works
The C++, Python, and Java implementations begin by building a smallest-prime-factor table up to the limit. That lets them factor each \(n\) quickly and compress repeated prime factors into the prime-power blocks \(q_i\) used in the derivation above.
If an \(n\) has only one distinct prime factor, the implementation immediately adds \(1\). Otherwise it enumerates every nonempty proper subset of the prime-power blocks. To make that fast, it precomputes the product attached to each subset, so each candidate only needs the subset product \(u\), the complementary factor \(v=n/u\), one modular inverse of \(u\) modulo \(v\), and one multiplication to recover \(a_U\). The maximum candidate is then added to the global sum. The Java implementation uses the same per-\(n\) mathematics inside a parallel outer summation.
Complexity Analysis
Building the smallest-prime-factor table with a linear sieve costs \(O(N)\) time and \(O(N)\) memory. For a fixed \(n\), the factorization step is cheap, and the dominant work is enumerating all subsets of the \(k=\omega(n)\) distinct prime factors. That contributes \(O(2^k)\) subset products and \(O(2^k \log n)\) time if the modular inverse is computed with the extended Euclidean algorithm.
A precise overall bound is therefore
$$O\left(N+\sum_{n=2}^{N} 2^{\omega(n)}\log n\right)$$
time and \(O(N)\) memory. For the actual bound \(N=10^7\), one has \(\omega(n)\le 8\), so no number ever needs more than \(2^8-2=254\) nontrivial subset checks. That is why the method is practical even though it explicitly searches all idempotent patterns for each \(n\).
Footnotes and References
- Problem page: https://projecteuler.net/problem=407
- Chinese remainder theorem: Wikipedia — Chinese remainder theorem
- Idempotent elements: Wikipedia — Idempotent (ring theory)
- Modular multiplicative inverse: Wikipedia — Modular multiplicative inverse
- Fundamental theorem of arithmetic: Wikipedia — Fundamental theorem of arithmetic
Problem 407 source code
C++
#include <cstdint>
#include <iostream>
#include <string>
#include <vector>
namespace {
using u64 = std::uint64_t;
using i64 = long long;
struct Options {
int limit = 10000000;
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;
}
try {
value = std::stoi(tail);
} catch (...) {
return false;
}
return true;
}
bool parse_arguments(int argc, char** argv, Options& options) {
for (int i = 1; i < argc; ++i) {
const 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 >= 1;
}
std::vector<int> build_spf(const int n) {
std::vector<int> spf(static_cast<std::size_t>(n + 1), 0);
std::vector<int> primes;
primes.reserve(n / 10);
for (int i = 2; i <= n; ++i) {
if (spf[static_cast<std::size_t>(i)] == 0) {
spf[static_cast<std::size_t>(i)] = i;
primes.push_back(i);
}
for (const int p : primes) {
const i64 x = 1LL * i * p;
if (x > n) {
break;
}
spf[static_cast<std::size_t>(x)] = p;
if (p == spf[static_cast<std::size_t>(i)]) {
break;
}
}
}
return spf;
}
int inv_mod(int a, int mod) {
int b = mod;
int u = 1;
int v = 0;
while (b != 0) {
const int t = a / b;
a -= t * b;
std::swap(a, b);
u -= t * v;
std::swap(u, v);
}
if (u < 0) {
u += mod;
}
return u;
}
u64 sum_M(const int limit) {
const std::vector<int> spf = build_spf(limit);
std::vector<int> prime_powers(10, 0);
std::vector<int> subset_prod(1 << 10, 1);
u64 ans = 0ULL;
for (int n = 1; n <= limit; ++n) {
if (n == 1) {
continue;
}
int x = n;
int k = 0;
while (x > 1) {
const int p = spf[static_cast<std::size_t>(x)];
int pe = 1;
while (x % p == 0) {
x /= p;
pe *= p;
}
prime_powers[static_cast<std::size_t>(k++)] = pe;
}
if (k == 1) {
ans += 1ULL;
continue;
}
const int full = (1 << k);
subset_prod[0] = 1;
for (int mask = 1; mask < full; ++mask) {
const int lsb = mask & -mask;
const int bit = __builtin_ctz(static_cast<unsigned>(lsb));
subset_prod[mask] = subset_prod[mask ^ lsb] * prime_powers[static_cast<std::size_t>(bit)];
}
int best = 1;
for (int mask = 1; mask < full - 1; ++mask) {
const int u = subset_prod[mask];
const int v = n / u;
const int inv = inv_mod(u % v, v);
const int a = static_cast<int>((1LL * u * inv) % n);
if (a > best) {
best = a;
}
}
ans += static_cast<u64>(best);
}
return ans;
}
bool run_checkpoints() {
// From statement: M(6)=4.
const u64 s6 = sum_M(6) - sum_M(5);
if (s6 != 4ULL) {
std::cerr << "Checkpoint failed: M(6)\n";
return false;
}
// Small consistency: M(1)=0 and M(prime)=1.
const u64 s2 = sum_M(2);
if (s2 != 1ULL) {
std::cerr << "Checkpoint failed: sum up to 2\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 << sum_M(options.limit) << '\n';
return 0;
}
Python
def solve():
LIMIT = 10_000_000
# Build SPF via linear sieve
spf = [0] * (LIMIT + 1)
primes = []
for i in range(2, LIMIT + 1):
if spf[i] == 0:
spf[i] = i
primes.append(i)
for p in primes:
if i * p > LIMIT: break
spf[i * p] = p
if p == spf[i]: break
def inv_mod(a, m):
b, u, v = m, 1, 0
while b:
t = a // b
a, b = b, a - t * b
u, v = v, u - t * v
return u % m
ans = 0
for n in range(2, LIMIT + 1):
x = n
k = 0
pp = []
while x > 1:
p = spf[x]
pe = 1
while x % p == 0:
x //= p
pe *= p
pp.append(pe)
k += 1
if k == 1:
ans += 1
continue
full = 1 << k
subset_prod = [0] * full
subset_prod[0] = 1
for mask in range(1, full):
lsb = mask & -mask
bit = lsb.bit_length() - 1
subset_prod[mask] = subset_prod[mask ^ lsb] * pp[bit]
best = 1
for mask in range(1, full - 1):
u = subset_prod[mask]
v = n // u
inv = inv_mod(u % v, v)
a = (u * inv) % n
if a > best:
best = a
ans += best
return str(ans)
if __name__ == '__main__':
print(solve())
Java
import java.util.ArrayList;
import java.util.List;
import java.util.stream.IntStream;
public class Euler407 {
static int[] buildSpf(int n) {
int[] spf = new int[n + 1];
List<Integer> primes = new ArrayList<>(n / 10);
for (int i = 2; i <= n; i++) {
if (spf[i] == 0) {
spf[i] = i;
primes.add(i);
}
for (int p : primes) {
long x = (long) i * p;
if (x > n)
break;
spf[(int) x] = p;
if (p == spf[i])
break;
}
}
return spf;
}
static int invMod(int a, int mod) {
int b = mod;
int u = 1;
int v = 0;
while (b != 0) {
int t = a / b;
int tempA = a - t * b;
a = b;
b = tempA;
int tempU = u - t * v;
u = v;
v = tempU;
}
if (u < 0)
u += mod;
return u;
}
static String solve() {
int limit = 10000000;
int[] spf = buildSpf(limit);
long total = IntStream.rangeClosed(1, limit)
.parallel()
.mapToLong(n -> {
if (n == 1)
return 0;
int[] primePowers = new int[10];
int x = n;
int k = 0;
while (x > 1) {
int p = spf[x];
int pe = 1;
while (x % p == 0) {
x /= p;
pe *= p;
}
primePowers[k++] = pe;
}
if (k == 1)
return 1;
int full = 1 << k;
int[] subsetProd = new int[full];
subsetProd[0] = 1;
for (int mask = 1; mask < full; mask++) {
int lsb = mask & -mask;
int bit = Integer.numberOfTrailingZeros(lsb);
subsetProd[mask] = subsetProd[mask ^ lsb] * primePowers[bit];
}
int best = 1;
for (int mask = 1; mask < full - 1; mask++) {
int u = subsetProd[mask];
int v = n / u;
int inv = invMod(u % v, v);
int a = (int) (((long) u * inv) % n);
if (a > best)
best = a;
}
return best;
})
.sum();
return Long.toString(total);
}
public static void main(String[] args) {
System.out.println(solve());
}
}