Problem 445: Retractions A
View on Project EulerProject Euler Problem 445 Solution
EulerSolve provides an optimized solution for Project Euler Problem 445, Retractions A, with C++, Python, Java, and a step-by-step mathematical explanation.
Problem Summary For each integer \(n \gt 1\), consider the linear maps $$f(x)\equiv ax+b \pmod{n},\qquad 0 \lt a \lt n,\quad 0\le b \lt n.$$ The map is a retraction if it is idempotent on every residue class: $$f(f(x))\equiv f(x)\pmod{n}\qquad\text{for all }x.$$ Let \(R(n)\) be the number of such maps. The target of the problem is $$S(N)=\sum_{k=1}^{N-1} R\left(\binom{N}{k}\right),\qquad N=10^7,$$ and the final result is required modulo \(M=10^9+7\). The main challenge is that the binomial coefficients are enormous, so the implementations never factor them from scratch. Mathematical Approach Step 1: Translate the retraction condition into divisibility conditions Starting from \(f(x)=ax+b\), we compose once more: $$f(f(x))\equiv a(ax+b)+b \equiv a^2x+ab+b \pmod{n}.$$ For this to equal \(ax+b\) for every \(x\), the coefficient of \(x\) and the constant correction must both vanish modulo \(n\). Therefore $$n\mid a(a-1),\qquad n\mid ab.$$ So every retraction is determined by two simultaneous conditions: \(a\) must satisfy an idempotence congruence, and once \(a\) is fixed, \(b\) must solve a linear congruence. Step 2: Solve the condition on \(a\) prime power by prime power Write the modulus as $$n=\prod_{i=1}^{m} p_i^{e_i}.$$ Because consecutive integers are coprime, \(\gcd(a,a-1)=1\)....
Detailed mathematical approach
Problem Summary
For each integer \(n \gt 1\), consider the linear maps
$$f(x)\equiv ax+b \pmod{n},\qquad 0 \lt a \lt n,\quad 0\le b \lt n.$$
The map is a retraction if it is idempotent on every residue class:
$$f(f(x))\equiv f(x)\pmod{n}\qquad\text{for all }x.$$
Let \(R(n)\) be the number of such maps. The target of the problem is
$$S(N)=\sum_{k=1}^{N-1} R\left(\binom{N}{k}\right),\qquad N=10^7,$$
and the final result is required modulo \(M=10^9+7\). The main challenge is that the binomial coefficients are enormous, so the implementations never factor them from scratch.
Mathematical Approach
Step 1: Translate the retraction condition into divisibility conditions
Starting from \(f(x)=ax+b\), we compose once more:
$$f(f(x))\equiv a(ax+b)+b \equiv a^2x+ab+b \pmod{n}.$$
For this to equal \(ax+b\) for every \(x\), the coefficient of \(x\) and the constant correction must both vanish modulo \(n\). Therefore
$$n\mid a(a-1),\qquad n\mid ab.$$
So every retraction is determined by two simultaneous conditions: \(a\) must satisfy an idempotence congruence, and once \(a\) is fixed, \(b\) must solve a linear congruence.
Step 2: Solve the condition on \(a\) prime power by prime power
Write the modulus as
$$n=\prod_{i=1}^{m} p_i^{e_i}.$$
Because consecutive integers are coprime, \(\gcd(a,a-1)=1\). Hence if \(p_i^{e_i}\mid a(a-1)\), the full prime power \(p_i^{e_i}\) must divide exactly one of the two factors. For each \(p_i^{e_i}\parallel n\), we must have
$$a\equiv 0 \pmod{p_i^{e_i}}\qquad\text{or}\qquad a\equiv 1 \pmod{p_i^{e_i}}.$$
These choices are independent for distinct prime powers, and the Chinese remainder theorem combines them into exactly one residue class modulo \(n\).
If we define
$$d=\gcd(a,n),$$
then \(d\) is a unitary divisor of \(n\): it satisfies \(d\mid n\) and \(\gcd(d,n/d)=1\). Conversely, every unitary divisor \(d\lt n\) determines one admissible residue \(a\). The excluded case \(d=n\) corresponds to \(a\equiv 0\pmod{n}\), but the problem requires \(0 \lt a \lt n\).
Step 3: Count the admissible values of \(b\)
Fix one valid \(a\), and let \(d=\gcd(a,n)\). Write
$$a=d\,a_1,\qquad n=d\,n_1,\qquad \gcd(a_1,n_1)=1.$$
The second condition becomes
$$n\mid ab \iff d\,n_1 \mid d\,a_1 b \iff n_1 \mid a_1 b.$$
Since \(a_1\) is invertible modulo \(n_1\), this is equivalent to
$$n_1\mid b.$$
Among the residues \(0\le b\lt n=d\,n_1\), exactly \(d\) values are multiples of \(n_1\). Therefore each admissible \(a\) contributes exactly \(d=\gcd(a,n)\) choices for \(b\).
Summing over all unitary divisors gives
$$R(n)=\sum_{\substack{d\parallel n\\ d\lt n}} d = \left(\sum_{d\parallel n} d\right)-n.$$
Thus
$$\boxed{R(n)=\sigma^*(n)-n},$$
where \(\sigma^*(n)\) is the sum of unitary divisors. Because each prime power is either chosen completely or not chosen at all,
$$\sigma^*(n)=\sum_{d\parallel n} d=\prod_{p^e\parallel n}(1+p^e).$$
Step 4: Apply the formula to the binomial coefficients
Let
$$B_k=\binom{N}{k}.$$
If
$$B_k=\prod_p p^{e_p(k)},$$
then the closed form above becomes
$$R(B_k)=\prod_{e_p(k)\gt 0}\left(1+p^{e_p(k)}\right)-B_k.$$
So the entire problem reduces to maintaining the prime exponents of \(B_k\) while \(k\) moves across the row of Pascal's triangle.
Step 5: Update \(\binom{N}{k}\) incrementally
The standard recurrence
$$B_k=B_{k-1}\cdot \frac{N-k+1}{k}$$
implies, for every prime \(p\),
$$e_p(k)=e_p(k-1)+v_p(N-k+1)-v_p(k).$$
This means the implementations only factor the two ordinary integers \(N-k+1\) and \(k\) at each step. They do not factor \(\binom{N}{k}\) itself.
At the same time they maintain
$$Q_k=\prod_{e_p(k)\gt 0}\left(1+p^{e_p(k)}\right)\pmod{M},$$
and then compute
$$R(B_k)\equiv Q_k-B_k \pmod{M}.$$
Because \(M\) is prime and \(N \lt M\), the recurrence for \(B_k \bmod M\) can also use modular inverses:
$$B_k \equiv B_{k-1}(N-k+1)k^{-1}\pmod{M}.$$
Step 6: Use symmetry and check a small example
Binomial coefficients satisfy
$$\binom{N}{k}=\binom{N}{N-k},$$
so the sum is symmetric. The implementations only iterate up to \(\lfloor N/2\rfloor\), using weight \(2\) except for the middle term when \(N\) is even.
For \(N=10\), the distinct values are
$$\binom{10}{1}=10,\quad \binom{10}{2}=45,\quad \binom{10}{3}=120,\quad \binom{10}{4}=210,\quad \binom{10}{5}=252.$$
Applying \(R(n)=\sigma^*(n)-n\):
$$\begin{aligned} R(10)&=(1+2)(1+5)-10=8,\\ R(45)&=(1+9)(1+5)-45=15,\\ R(120)&=(1+8)(1+3)(1+5)-120=96,\\ R(210)&=(1+2)(1+3)(1+5)(1+7)-210=366,\\ R(252)&=(1+4)(1+9)(1+7)-252=148. \end{aligned}$$
Hence
$$S(10)=2(8+15+96+366)+148=1118,$$
which is the small checkpoint matched by the implementation.
How the Code Works
The C++, Python, and Java implementations first build a smallest-prime-factor table up to \(N\). They also precompute modular inverses \(1^{-1},2^{-1},\dots,N^{-1}\pmod M\) so that the current binomial coefficient can be updated in constant time once the numerator and denominator factors are known.
During the main loop, the implementation stores for each active prime its current exponent in \(B_k\) and the current value of \(p^{e_p(k)} \bmod M\). When a prime exponent changes, the old factor \((1+p^e)\) is removed from the running product \(Q_k\), the exponent is adjusted, and the new factor is inserted.
Division inside the modulus is handled by multiplying with inverses. The implementations cache inverses of previously seen nonzero factors, and they also keep track of how many factors are currently \(0 \pmod M\). If at least one such factor exists, then \(Q_k\equiv 0\pmod M\), so the running product remains correct even in that edge case.
Complexity Analysis
The smallest-prime-factor preprocessing uses \(O(N)\) memory. The total factor work over all updates is governed by the number of prime factors, counted with multiplicity, appearing in the integers \(1,2,\dots,N\), which is near \(O(N\log\log N)\) in aggregate. The remaining arithmetic per step is constant. Therefore the full method runs in near \(O(N\log\log N)\) time and \(O(N)\) memory, which is practical for \(N=10^7\).
Footnotes and References
- Problem page: https://projecteuler.net/problem=445
- Chinese remainder theorem: Wikipedia — Chinese remainder theorem
- Unitary divisor: Wikipedia — Unitary divisor
- Binomial coefficient: Wikipedia — Binomial coefficient
- Hardy and Wright, An Introduction to the Theory of Numbers, sections on multiplicative arithmetic functions and congruences.
Problem 445 source code
C++
#include <cstdint>
#include <iostream>
#include <limits>
#include <string>
#include <unordered_map>
#include <vector>
#include <algorithm>
#include <functional>
namespace {
using u32 = std::uint32_t;
using u64 = std::uint64_t;
constexpr u32 kMod = 1'000'000'007U;
struct Options {
u32 n = 10'000'000U;
bool run_checkpoints = true;
};
bool parse_u32_after_prefix(const std::string& arg, const std::string& prefix, u32& out) {
if (arg.rfind(prefix, 0U) != 0U) {
return false;
}
const std::string tail = arg.substr(prefix.size());
if (tail.empty()) {
return false;
}
try {
const unsigned long long parsed = std::stoull(tail);
if (parsed > static_cast<unsigned long long>(std::numeric_limits<u32>::max())) {
return false;
}
out = static_cast<u32>(parsed);
} 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_u32_after_prefix(arg, "--n=", options.n)) {
continue;
}
std::cerr << "Unknown argument: " << arg << '\n';
return false;
}
return options.n >= 2U;
}
u32 mod_pow(u64 base, u64 exp) {
u64 result = 1ULL;
base %= kMod;
while (exp > 0ULL) {
if ((exp & 1ULL) != 0ULL) {
result = (result * base) % kMod;
}
base = (base * base) % kMod;
exp >>= 1ULL;
}
return static_cast<u32>(result);
}
class TermProductTracker {
public:
void add(u32 value) {
if (value == 0U) {
++zero_terms_;
return;
}
product_non_zero_ = (product_non_zero_ * value) % kMod;
}
void remove(u32 value) {
if (value == 0U) {
--zero_terms_;
return;
}
product_non_zero_ = (product_non_zero_ * inverse(value)) % kMod;
}
u32 value() const {
return (zero_terms_ > 0) ? 0U : static_cast<u32>(product_non_zero_);
}
private:
u64 product_non_zero_ = 1ULL;
int zero_terms_ = 0;
std::unordered_map<u32, u32> inverse_cache_;
u32 inverse(const u32 x) {
const auto it = inverse_cache_.find(x);
if (it != inverse_cache_.end()) {
return it->second;
}
const u32 inv = mod_pow(x, static_cast<u64>(kMod) - 2ULL);
inverse_cache_.emplace(x, inv);
return inv;
}
};
u32 solve(const u32 n) {
std::vector<u32> spf(static_cast<std::size_t>(n) + 1U, 0U);
std::vector<u32> primes;
primes.reserve(static_cast<std::size_t>(n / 10U));
for (u32 i = 2U; i <= n; ++i) {
if (spf[i] == 0U) {
spf[i] = i;
primes.push_back(i);
}
for (const u32 p : primes) {
const u64 v = static_cast<u64>(i) * static_cast<u64>(p);
if (v > n || p > spf[i]) {
break;
}
spf[static_cast<std::size_t>(v)] = p;
}
}
std::vector<u32> inv(static_cast<std::size_t>(n) + 1U, 0U);
inv[1] = 1U;
for (u32 i = 2U; i <= n; ++i) {
inv[i] = static_cast<u32>(kMod - (static_cast<u64>(kMod / i) * inv[kMod % i]) % kMod);
}
std::vector<int> exponent(static_cast<std::size_t>(n) + 1U, 0);
std::vector<u32> power_mod(static_cast<std::size_t>(n) + 1U, 1U);
TermProductTracker tracker;
const auto adjust_prime = [&](const u32 p, const int delta) {
int old_exp = exponent[p];
u32 p_pow = power_mod[p];
if (old_exp > 0) {
tracker.remove((p_pow + 1U) % kMod);
}
if (delta > 0) {
for (int i = 0; i < delta; ++i) {
p_pow = static_cast<u32>((static_cast<u64>(p_pow) * p) % kMod);
}
} else if (delta < 0) {
const u32 inv_p = inv[p];
for (int i = 0; i < -delta; ++i) {
p_pow = static_cast<u32>((static_cast<u64>(p_pow) * inv_p) % kMod);
}
}
old_exp += delta;
exponent[p] = old_exp;
power_mod[p] = p_pow;
if (old_exp > 0) {
tracker.add((p_pow + 1U) % kMod);
}
};
const auto apply_factorization = [&](u32 x, const int sign) {
while (x > 1U) {
const u32 p = spf[x];
int cnt = 0;
while (x % p == 0U) {
x /= p;
++cnt;
}
adjust_prime(p, sign * cnt);
}
};
u64 binom_mod = 1ULL; // C(n,0)
u64 answer = 0ULL;
const u32 half = n / 2U;
for (u32 k = 1U; k <= half; ++k) {
apply_factorization(n - k + 1U, +1);
apply_factorization(k, -1);
binom_mod = (binom_mod * (n - k + 1ULL)) % kMod;
binom_mod = (binom_mod * inv[k]) % kMod;
const u32 q_mod = tracker.value();
const u32 r_mod = static_cast<u32>((q_mod + kMod - binom_mod) % kMod);
const u32 weight = (k == n - k) ? 1U : 2U;
answer += static_cast<u64>(weight) * r_mod;
answer %= kMod;
}
return static_cast<u32>(answer);
}
bool run_checkpoints() {
if (solve(10U) != 1118U) {
std::cerr << "Checkpoint failed: N=10\n";
return false;
}
if (solve(100'000U) != 628'701'600U) {
std::cerr << "Checkpoint failed: N=100000\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.n) << '\n';
return 0;
}
Python
import math
def solve():
MOD = 1000000007
n = 10000000
def mod_pow(base, exp):
r = 1; base %= MOD
while exp > 0:
if exp & 1: r = r * base % MOD
base = base * base % MOD
exp >>= 1
return r
# SPF sieve
spf = list(range(n + 1))
for p in range(2, int(n**0.5) + 1):
if spf[p] == p:
for q in range(p*p, n+1, p):
if spf[q] == q: spf[q] = p
inv = [0] * (n + 1)
inv[1] = 1
for i in range(2, n + 1):
inv[i] = (MOD - (MOD // i) * inv[MOD % i] % MOD) % MOD
exponent = [0] * (n + 1)
power_mod = [1] * (n + 1)
# Product of (p^e + 1) tracker
product = 1; zeros = 0
inv_cache = {}
def get_inv(x):
if x not in inv_cache: inv_cache[x] = mod_pow(x, MOD - 2)
return inv_cache[x]
def adjust(p, delta):
nonlocal product, zeros
old_exp = exponent[p]; pp = power_mod[p]
if old_exp > 0:
v = (pp + 1) % MOD
if v == 0: zeros -= 1
else: product = product * get_inv(v) % MOD
if delta > 0:
for _ in range(delta): pp = pp * p % MOD
else:
ip = inv[p]
for _ in range(-delta): pp = pp * ip % MOD
exponent[p] = old_exp + delta
power_mod[p] = pp
if old_exp + delta > 0:
v = (pp + 1) % MOD
if v == 0: zeros += 1
else: product = product * v % MOD
def apply_fact(x, sign):
while x > 1:
p = spf[x]; cnt = 0
while x % p == 0: x //= p; cnt += 1
adjust(p, sign * cnt)
binom_mod = 1; answer = 0; half = n // 2
for k in range(1, half + 1):
apply_fact(n - k + 1, 1)
apply_fact(k, -1)
binom_mod = binom_mod * (n - k + 1) % MOD * inv[k] % MOD
q_mod = 0 if zeros > 0 else product
r_mod = (q_mod + MOD - binom_mod) % MOD
weight = 1 if k == n - k else 2
answer = (answer + weight * r_mod) % MOD
return str(answer)
if __name__ == '__main__':
print(solve())
Java
import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;
public class Euler445 {
static final int MOD = 1000000007;
static long modPow(long base, long exp) {
long result = 1;
base %= MOD;
while (exp > 0) {
if ((exp & 1) != 0) {
result = (result * base) % MOD;
}
base = (base * base) % MOD;
exp >>= 1;
}
return result;
}
static class TermProductTracker {
long productNonZero = 1;
int zeroTerms = 0;
Map<Integer, Integer> inverseCache = new HashMap<>();
void add(int value) {
if (value == 0) {
zeroTerms++;
return;
}
productNonZero = (productNonZero * value) % MOD;
}
void remove(int value) {
if (value == 0) {
zeroTerms--;
return;
}
int inv = inverseCache.getOrDefault(value, 0);
if (inv == 0) {
inv = (int) modPow(value, MOD - 2);
inverseCache.put(value, inv);
}
productNonZero = (productNonZero * inv) % MOD;
}
int value() {
return zeroTerms > 0 ? 0 : (int) productNonZero;
}
}
public static String solve() {
int n = 10000000;
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 v = (long) i * p;
if (v > n || p > spf[i]) {
break;
}
spf[(int) v] = p;
}
}
int[] inv = new int[n + 1];
inv[1] = 1;
for (int i = 2; i <= n; ++i) {
inv[i] = (int) (MOD - ((long) (MOD / i) * inv[MOD % i]) % MOD);
}
int[] exponent = new int[n + 1];
int[] powerMod = new int[n + 1];
for (int i = 0; i <= n; i++)
powerMod[i] = 1;
TermProductTracker tracker = new TermProductTracker();
long binomMod = 1;
long answer = 0;
int half = n / 2;
for (int k = 1; k <= half; ++k) {
int x = n - k + 1;
while (x > 1) {
int p = spf[x];
int cnt = 0;
while (x % p == 0) {
x /= p;
cnt++;
}
int delta = cnt;
int oldExp = exponent[p];
long pPow = powerMod[p];
if (oldExp > 0) {
tracker.remove((int) ((pPow + 1) % MOD));
}
for (int i = 0; i < delta; ++i) {
pPow = (pPow * p) % MOD;
}
exponent[p] = oldExp + delta;
powerMod[p] = (int) pPow;
if (exponent[p] > 0) {
tracker.add((int) ((pPow + 1) % MOD));
}
}
x = k;
while (x > 1) {
int p = spf[x];
int cnt = 0;
while (x % p == 0) {
x /= p;
cnt++;
}
int delta = -cnt;
int oldExp = exponent[p];
long pPow = powerMod[p];
if (oldExp > 0) {
tracker.remove((int) ((pPow + 1) % MOD));
}
long invP = inv[p];
for (int i = 0; i < -delta; ++i) {
pPow = (pPow * invP) % MOD;
}
exponent[p] = oldExp + delta;
powerMod[p] = (int) pPow;
if (exponent[p] > 0) {
tracker.add((int) ((pPow + 1) % MOD));
}
}
binomMod = (binomMod * (n - k + 1)) % MOD;
binomMod = (binomMod * inv[k]) % MOD;
long qMod = tracker.value();
long rMod = (qMod + MOD - binomMod) % MOD;
long weight = (k == n - k) ? 1 : 2;
answer = (answer + weight * rMod) % MOD;
}
return Long.toString(answer);
}
public static void main(String[] args) {
System.out.println(solve());
}
}