Detailed mathematical approach

Problem Summary

For each integer \(n \gt 1\), consider 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 when it is idempotent on every residue class:

$f(f(x))\equiv f(x)\pmod{n}\qquad \forall x.$

Let \(R(n)\) be the number of such maps. Problem 447 asks for

$F(N)=\sum_{n=2}^{N}R(n)\pmod{10^9+7},$

with the checkpoint \(F(10^7)\equiv 638042271 \pmod{10^9+7}\) and target \(N=10^{14}\). Direct enumeration is impossible at this scale, so the implementation rewrites \(R(n)\) as an arithmetic function with a fast summatory formula.

Mathematical Approach

Step 1: Express \(R(n)\) through unitary divisors

Composing \(f(x)=ax+b\) once more gives

$f(f(x))\equiv a^2x+ab+b \pmod{n}.$

For this to equal \(ax+b\) for every residue class, both corrections must vanish modulo \(n\):

$n \mid a(a-1),\qquad n \mid ab.$

If we write \(d=\gcd(a,n)\), then \(d\mid n\) and \(\gcd(d,n/d)=1\), so \(d\) is a unitary divisor of \(n\). Conversely, each unitary divisor \(d \lt n\) determines one admissible residue class for \(a\), and once \(a\) is fixed, the number of valid \(b\) values is exactly \(d\). Therefore

$R(n)=\sum_{\substack{d \mid n\\ \gcd(d,n/d)=1\\ d \lt n}} d=\sigma^*(n)-n,$

where \(\sigma^*(n)\) is the sum of unitary divisors of \(n\). Since \(\sigma^*(1)=1\), the missing term at \(n=1\) contributes \(1-1=0\), so

$F(N)=\sum_{n=1}^{N}\bigl(\sigma^*(n)-n\bigr)=\left(\sum_{n=1}^{N}\sigma^*(n)\right)-\frac{N(N+1)}{2}.$

Step 2: Replace \(\sigma^*(n)\) by a Möbius square decomposition

The unitary condition can be encoded by Möbius inversion:

$\mathbf{1}_{\gcd(d,n/d)=1}=\sum_{t \mid \gcd(d,n/d)} \mu(t).$

Insert this into the divisor sum:

$\sigma^*(n)=\sum_{d \mid n} d \sum_{t \mid \gcd(d,n/d)} \mu(t).$

If \(t \mid \gcd(d,n/d)\), then \(t^2 \mid n\), and we may write \(d=t\,m\) with \(m \mid n/t^2\). Hence

$\sigma^*(n)=\sum_{t^2 \mid n} \mu(t)\, t \sum_{m \mid n/t^2} m=\sum_{t^2 \mid n} \mu(t)\, t\, \sigma\!\left(\frac{n}{t^2}\right).$

This identity is the bridge from unitary divisors to the ordinary sum-of-divisors function \(\sigma\), which is much easier to sum over long intervals.

Step 3: Convert the prefix sum into a summatory \(\sigma\) problem

Define

$T(x)=\sum_{m=1}^{x}\sigma(m).$

Summing the previous identity for \(n \le N\) and swapping the order of summation gives

$\sum_{n=1}^{N}\sigma^*(n)=\sum_{t \le \sqrt{N}} \mu(t)\, t\, T\!\left(\left\lfloor\frac{N}{t^2}\right\rfloor\right).$

Now use the classical divisor-summatory identity

$T(x)=\sum_{m=1}^{x}\sum_{d \mid m} d=\sum_{d=1}^{x} d\left\lfloor\frac{x}{d}\right\rfloor.$

So the original problem has been reduced to evaluating many instances of \(T(x)\) quickly.

Step 4: Evaluate \(T(x)\) by quotient blocks

The value \(\left\lfloor x/i \right\rfloor\) stays constant on intervals. If \(l\) is the first index of a block, define

$q=\left\lfloor\frac{x}{l}\right\rfloor,\qquad r=\left\lfloor\frac{x}{q}\right\rfloor.$

Then every \(i\in[l,r]\) has the same quotient \(q\), so the whole block contributes

$q\sum_{i=l}^{r} i=q\cdot \frac{(l+r)(r-l+1)}{2}.$

Advancing directly from \(l\) to \(r+1\) skips an entire constant-quotient range, reducing the cost of \(T(x)\) from \(O(x)\) to about \(O(\sqrt{x})\). This quotient-block decomposition is the central speedup in the implementation.

Step 5: Final formula and a worked check

Putting everything together,

$\boxed{F(N)\equiv \sum_{t \le \sqrt{N}} \mu(t)\, t\, T\!\left(\left\lfloor\frac{N}{t^2}\right\rfloor\right)-\frac{N(N+1)}{2}\pmod{10^9+7}.}$

For a small check, take \(N=10\). Quotient blocks give

$T(10)=1\cdot 10+2\cdot 5+3\cdot 3+(4+5)\cdot 2+(6+7+8+9+10)\cdot 1=87.$

Also \(T(2)=4\), \(T(1)=1\), and the nonzero Möbius values up to \(\sqrt{10}\) are \(\mu(1)=1\), \(\mu(2)=-1\), \(\mu(3)=-1\). Therefore

$\sum_{n=1}^{10}\sigma^*(n)=1\cdot 87-2\cdot 4-3\cdot 1=76,$

so

$F(10)=76-\frac{10\cdot 11}{2}=21.$

This agrees with direct evaluation of \(R(2),\dots,R(10)\), so the derivation is consistent before moving to \(10^{14}\).

How the Code Works

The C++, Python, and Java implementations first sieve the Möbius function up to \(\lfloor\sqrt{N}\rfloor\). The outer sum then skips every \(t\) with \(\mu(t)=0\), because such terms contribute nothing. For each remaining \(t\), the implementation evaluates one inner query \(T\!\left(\lfloor N/t^2\rfloor\right)\) using quotient blocks instead of a linear loop.

After the summatory unitary-divisor part has been accumulated, the implementation subtracts the triangular number \(N(N+1)/2\) modulo \(10^9+7\). The C++ and Java implementations split the outer Möbius sum cyclically across worker threads so that expensive small \(t\) values are balanced more evenly; the Python implementation uses the same mathematics sequentially. In every language, modular subtraction is normalized back into the range \(0,\dots,10^9+6\).

Complexity Analysis

The Möbius sieve up to \(\sqrt{N}\) needs \(O(\sqrt{N})\) memory and near-linear preprocessing in that range. The dominant work is

$\sum_{t \le \sqrt{N}} O\!\left(\sqrt{\frac{N}{t^2}}\right)=\sum_{t \le \sqrt{N}} O\!\left(\frac{\sqrt{N}}{t}\right)=O(\sqrt{N}\log N).$

So the overall method runs in \(O(\sqrt{N}\log N)\) time and \(O(\sqrt{N})\) memory. Parallel execution improves wall-clock time in the compiled implementations, but the asymptotic bound is unchanged.

Footnotes and References

1. Problem page: https://projecteuler.net/problem=447
2. Unitary divisor: Wikipedia - Unitary divisor
3. Möbius function and inversion: Wikipedia - Möbius function
4. Divisor summatory function: Wikipedia - Divisor summatory function
5. T. M. Apostol, Introduction to Analytic Number Theory, chapters on multiplicative arithmetic functions and Möbius inversion.