Problem 447: Retractions C

View on Project Euler

Project Euler Problem 447 Solution

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

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\)....

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.

Problem 447 source code

C++

#include <iostream>
#include <vector>
#include <cmath>
#include <future>
#include <thread>
#include <numeric>

// Type definitions for handling large numbers
using int128 = __int128;
using ll = long long;

const ll MOD = 1000000007;

// Modular arithmetic helper to handle negative results
ll safe_mod(int128 val) {
    val %= MOD;
    if (val < 0) val += MOD;
    return (ll)val;
}

// Modular inverse of 2 for division
const ll INV2 = 500000004; // (MOD + 1) / 2

// Class to handle the Sieve of Eratosthenes for Mobius function
class Sieve {
public:
    std::vector<int> mu;
    
    Sieve(int n) {
        mu.resize(n + 1);
        std::vector<bool> is_prime(n + 1, true);
        std::vector<int> primes;
        
        mu[1] = 1;
        for (int i = 2; i <= n; ++i) {
            if (is_prime[i]) {
                primes.push_back(i);
                mu[i] = -1;
            }
            for (int p : primes) {
                if (i * p > n) break;
                is_prime[i * p] = false;
                if (i % p == 0) {
                    mu[i * p] = 0;
                    break;
                }
                mu[i * p] = -mu[i];
            }
        }
    }
};

// Calculate sum of sigma_1(k) for k=1 to x
// T(x) = sum_{i=1}^x i * floor(x/i)
ll calc_T(ll x) {
    ll total_sum = 0;
    for (ll l = 1, r; l <= x; l = r + 1) {
        r = x / (x / l);
        // Sum of arithmetic progression from l to r: (l+r)*(r-l+1)/2
        // We do calculation in int128 to prevent overflow before modulo
        int128 count = (r - l + 1);
        int128 sum_l_r = (l + r);
        
        // (count * sum_l_r) / 2 % MOD
        ll sum_seq = safe_mod((count * sum_l_r) % (2 * MOD) / 2);
        
        int128 term = (int128)sum_seq * (x / l);
        total_sum = safe_mod(total_sum + term);
    }
    return total_sum;
}

// Solver function
ll solve(ll N, const Sieve& sieve) {
    ll limit = sqrt(N);
    
    // We need to compute sum_{d=1}^{limit} d * mu[d] * T(N/d^2)
    // We will parallelize this loop.
    
    int num_threads = std::thread::hardware_concurrency();
    if (num_threads == 0) num_threads = 4;
    
    std::vector<std::future<ll>> futures;
    
    auto task = [&](int thread_id) -> ll {
        ll local_sum = 0;
        // Cyclic distribution for load balancing (small d is expensive, large d is cheap)
        for (ll d = thread_id + 1; d <= limit; d += num_threads) {
            if (sieve.mu[d] == 0) continue;
            
            ll inner = calc_T(N / (d * d));
            int128 term = (int128)d * inner;
            term = safe_mod(term);
            
            if (sieve.mu[d] == 1) {
                local_sum = safe_mod(local_sum + term);
            } else {
                local_sum = safe_mod(local_sum - term);
            }
        }
        return local_sum;
    };
    
    for (int i = 0; i < num_threads; ++i) {
        futures.push_back(std::async(std::launch::async, task, i));
    }
    
    ll S_N = 0;
    for (auto& f : futures) {
        S_N = safe_mod(S_N + f.get());
    }
    
    // F(N) = S(N) - N(N+1)/2
    int128 n_sum = (int128)N * (N + 1);
    // Be careful with modulo division
    ll n_sum_mod = safe_mod(n_sum % (2 * MOD) / 2);
    
    return safe_mod(S_N - n_sum_mod);
}

int main() {
    // 1. Validation Check Point
    std::cout << "Running validation for N = 10^7..." << std::endl;
    
    // Precompute sieve up to sqrt(10^14) = 10^7
    // This covers the requirement for both 10^7 and 10^14 cases.
    Sieve sieve(10000000); 
    
    ll val_N = 10000000;
    ll val_result = solve(val_N, sieve);
    ll expected = 638042271;
    
    std::cout << "F(10^7) = " << val_result << std::endl;
    
    if (val_result == expected) {
        std::cout << "Validation PASSED. Proceeding to solve for 10^14..." << std::endl;
    } else {
        std::cout << "Validation FAILED. Expected " << expected << ". Aborting." << std::endl;
        return 1;
    }
    
    std::cout << "------------------------------------------------" << std::endl;
    
    // 2. Main Problem
    ll target_N = 100000000000000LL; // 10^14
    ll result = solve(target_N, sieve);
    
    std::cout << "F(10^14) modulo 10^9 + 7 is:" << std::endl;
    std::cout << result << std::endl;
    std::cout << "Answer: " << result << std::endl;
    
    return 0;
}

Python

import math

def solve():
    MOD = 1000000007
    N = 100000000000000
    INV2 = 500000004

    limit = int(math.isqrt(N))
    # Sieve Mobius function
    mu = [0] * (limit + 1)
    mu[1] = 1
    is_prime = bytearray(b'\x01') * (limit + 1)
    primes = []
    for i in range(2, limit + 1):
        if is_prime[i]:
            primes.append(i)
            mu[i] = -1
        for p in primes:
            if i * p > limit: break
            is_prime[i * p] = 0
            if i % p == 0:
                mu[i * p] = 0
                break
            mu[i * p] = -mu[i]

    def calc_T(x):
        total = 0
        l = 1
        while l <= x:
            r = x // (x // l)
            count = r - l + 1
            sum_lr = l + r
            sum_seq = (count * sum_lr % (2 * MOD) // 2) % MOD
            total = (total + sum_seq * (x // l)) % MOD
            l = r + 1
        return total

    S_N = 0
    for d in range(1, limit + 1):
        if mu[d] == 0: continue
        inner = calc_T(N // (d * d))
        term = d * inner % MOD
        if mu[d] == 1:
            S_N = (S_N + term) % MOD
        else:
            S_N = (S_N - term) % MOD

    n_sum = N % MOD * ((N + 1) % MOD) % MOD * INV2 % MOD
    return str((S_N - n_sum) % MOD)

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

Java

import java.util.ArrayList;
import java.util.List;
import java.util.stream.IntStream;

public class Euler447 {
    static final long MOD = 1000000007;

    static long safeMod(long val) {
        val %= MOD;
        if (val < 0)
            val += MOD;
        return val;
    }

    static class Sieve {
        byte[] mu;

        Sieve(int n) {
            mu = new byte[n + 1];
            boolean[] isPrime = new boolean[n + 1];
            for (int i = 2; i <= n; i++)
                isPrime[i] = true;

            List<Integer> primes = new ArrayList<>(n / 10);

            mu[1] = 1;
            for (int i = 2; i <= n; i++) {
                if (isPrime[i]) {
                    primes.add(i);
                    mu[i] = -1;
                }
                for (int p : primes) {
                    if ((long) i * p > n)
                        break;
                    isPrime[i * p] = false;
                    if (i % p == 0) {
                        mu[i * p] = 0;
                        break;
                    }
                    mu[i * p] = (byte) -mu[i];
                }
            }
        }
    }

    static long calcT(long x) {
        long totalSum = 0;
        long l = 1;
        while (l <= x) {
            long q = x / l;
            long r = x / q;

            long count = r - l + 1;
            long sumLR = l + r;

            long sumSeq = safeMod((count % (2 * MOD)) * (sumLR % (2 * MOD)) % (2 * MOD) / 2);

            long term = safeMod(sumSeq * (q % MOD));
            totalSum = safeMod(totalSum + term);

            l = r + 1;
        }
        return totalSum;
    }

    public static String solve() {
        long N = 100000000000000L;
        int limit = (int) Math.sqrt(N);

        Sieve sieve = new Sieve(limit);

        int threads = Math.max(1, Runtime.getRuntime().availableProcessors());

        long sN = IntStream.range(0, threads).parallel().mapToLong(threadId -> {
            long localSum = 0;
            for (int d = threadId + 1; d <= limit; d += threads) {
                if (sieve.mu[d] == 0)
                    continue;

                long inner = calcT(N / ((long) d * d));
                long term = safeMod((d % MOD) * inner);

                if (sieve.mu[d] == 1) {
                    localSum = safeMod(localSum + term);
                } else {
                    localSum = safeMod(localSum - term);
                }
            }
            return localSum;
        }).reduce(0, (a, b) -> safeMod(a + b));

        long nSumMod = safeMod((N % (2 * MOD)) * ((N + 1) % (2 * MOD)) % (2 * MOD) / 2);
        long result = safeMod(sN - nSumMod);

        return Long.toString(result);
    }

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