Problem 429: Sum of Squares of Unitary Divisors

Mathematical Analysis

Unitary Divisors and Multiplicative Functions

If $n = p_1^{a_1} p_2^{a_2} \cdots p_k^{a_k}$, then a unitary divisor of $n$ must be of the form $d = p_1^{b_1} p_2^{b_2} \cdots p_k^{b_k}$ where each $b_i \in \{0, a_i\}$ (either take the full prime power or nothing, to ensure $\gcd(d, n/d) = 1$).

The sum of squares of unitary divisors is therefore:

Applying Legendre’s Formula

To find the prime factorization of $N! = (10^8)!$, we use Legendre’s formula for the exponent of prime $p$ in $N!$:

Editorial

Key insight: If n = p1^a1 * p2^a2 * … then sigma_2(n) = prod(1 + p_i^(2a_i)) Use Legendre’s formula for exponents in factorial, and sieve for primes. We generate all primes up to $N = 10^8$ using a sieve of Eratosthenes. We then iterate over each prime $p$, compute $a_p = v_p(N!)$ using Legendre’s formula. Finally, compute the product $\prod_p (1 + p^{2a_p}) \pmod{10^9 + 9}$ using modular exponentiation.

Pseudocode

Generate all primes up to $N = 10^8$ using a sieve of Eratosthenes
For each prime $p$, compute $a_p = v_p(N!)$ using Legendre's formula
Compute the product $\prod_p (1 + p^{2a_p}) \pmod{10^9 + 9}$ using modular exponentiation

Correctness

Theorem. The method described above computes exactly the quantity requested in the problem statement.

Proof. The preceding analysis identifies the admissible objects and derives the formula, recurrence, or exhaustive search carried out by the algorithm. The computation evaluates exactly that specification, so every valid contribution is included once and no invalid contribution is counted. Therefore the returned value is the required answer. $\square$

Complexity Analysis

• Time: $O(N)$ for the sieve, plus $O(\pi(N) \cdot \log N)$ for the Legendre computations and modular exponentiations, where $\pi(N) \approx N/\ln N$.
• Space: $O(N)$ for the prime sieve.

Answer

#include <bits/stdc++.h>
using namespace std;

const long long MOD = 1000000009;
const int N = 100000000;

long long power(long long base, long long exp, long long mod) {
    long long result = 1;
    base %= mod;
    while (exp > 0) {
        if (exp & 1) result = result * base % mod;
        base = base * base % mod;
        exp >>= 1;
    }
    return result;
}

int main() {
    // Sieve of Eratosthenes
    vector<bool> is_prime(N + 1, true);
    is_prime[0] = is_prime[1] = false;
    for (long long i = 2; i * i <= N; i++) {
        if (is_prime[i]) {
            for (long long j = i * i; j <= N; j += i)
                is_prime[j] = false;
        }
    }

    long long result = 1;
    for (long long p = 2; p <= N; p++) {
        if (!is_prime[p]) continue;
        // Legendre's formula: exponent of p in N!
        long long a = 0;
        long long pk = p;
        while (pk <= N) {
            a += N / pk;
            if (pk > N / p) break; // prevent overflow
            pk *= p;
        }
        // Multiply by (1 + p^(2a)) mod MOD
        result = result % MOD * ((1 + power(p, 2 * a, MOD)) % MOD) % MOD;
    }

    cout << result << endl;
    return 0;
}

"""
Project Euler Problem 429: Sum of Squares of Unitary Divisors

Find S(100_000_000!) mod 1_000_000_009 where S(n) is the sum of
squares of unitary divisors of n.

Key insight: If n = p1^a1 * p2^a2 * ... then
sigma*_2(n) = prod(1 + p_i^(2*a_i))

Use Legendre's formula for exponents in factorial,
and sieve for primes.
"""

def solve():
    N = 100_000_000
    MOD = 1_000_000_009

    # Sieve of Eratosthenes
    is_prime = bytearray(b'\x01') * (N + 1)
    is_prime[0] = is_prime[1] = 0
    for i in range(2, int(N**0.5) + 1):
        if is_prime[i]:
            is_prime[i*i::i] = bytearray(len(is_prime[i*i::i]))

    result = 1
    for p in range(2, N + 1):
        if not is_prime[p]:
            continue
        # Legendre's formula
        a = 0
        pk = p
        while pk <= N:
            a += N // pk
            pk *= p
        # Multiply (1 + p^(2a)) mod MOD
        result = result * (1 + pow(p, 2 * a, MOD)) % MOD

    print(result)

if __name__ == '__main__':
    solve()