Problem 745: Sum of Squares II

Mathematical Foundation

Theorem 1 (Unique Squarefree Decomposition). Every positive integer $n$ can be written uniquely as $n = a^2 b$ where $b$ is squarefree. In this decomposition, $g(n) = a^2$.

Proof. Write $n = \prod p_i^{e_i}$. Set $a = \prod p_i^{\lfloor e_i / 2 \rfloor}$ and $b = \prod p_i^{e_i \bmod 2}$. Then $n = a^2 b$, $b$ is squarefree, and $a^2$ is the largest perfect square dividing $n$. Uniqueness follows from the uniqueness of prime factorization. $\square$

Theorem 2 (Jordan Totient Reformulation). Define the Jordan totient function of order $2$:

Then:

Proof. From Theorem 1:

where we used the standard Mobius inversion formula $Q(m) = \sum_{d=1}^{\lfloor \sqrt{m} \rfloor} \mu(d) \lfloor m/d^2 \rfloor$ for the squarefree counting function. Substituting $e = ad$:

The identity $\sum_{a \mid e} a^2 \mu(e/a) = J_2(e)$ is the definition of the Jordan totient as the Dirichlet convolution $\mathrm{id}_2 * \mu$. The multiplicative formula $J_2(e) = e^2 \prod_{p \mid e}(1 - 1/p^2)$ follows from multiplicativity and evaluation at prime powers: $J_2(p^k) = p^{2k} - p^{2(k-1)}$. $\square$

Lemma 1 (Sieve for $J_2$). The values $J_2(e)$ for $1 \le e \le L$ can be computed by a multiplicative sieve analogous to the Euler totient sieve:

1. Initialize $J_2(e) = e^2$ for all $e$.
2. For each prime $p \le L$: for each multiple $m$ of $p$, set $J_2(m) \leftarrow J_2(m) \cdot (p^2 - 1) / p^2$.

This runs in $O(L \log \log L)$ time.

Proof. The sieve correctly computes $J_2(e) = e^2 \prod_{p \mid e}(1 - 1/p^2)$ by multiplying in the factor $(1 - 1/p^2) = (p^2-1)/p^2$ for each prime divisor. The time complexity is the same as the Euler sieve. The division by $p^2$ is exact in modular arithmetic since we track the factor $(p^2-1)/p^2$ as a modular inverse multiplication. $\square$

Editorial

For a positive integer n, g(n) = max perfect square dividing n. S(N) = sum of g(n) for n=1..N. Find S(10^14) mod 10^9+7. Uses Jordan totient function J_2 sieve approach. We sieve: for each prime, multiply in (1 - 1/p^2) factor. Finally, compute S(N) = sum_{e=1}^{L} J2[e] * floor(N / e^2) mod p.

Pseudocode

Sieve J_2(e) mod p for e = 1..L
Sieve: for each prime, multiply in (1 - 1/p^2) factor
Compute S(N) = sum_{e=1}^{L} J2[e] * floor(N / e^2) mod p

Complexity Analysis

• Time: $O(\sqrt{N} \log \log \sqrt{N})$ for the sieve, $O(\sqrt{N})$ for the final summation. Total: $O(\sqrt{N} \log \log \sqrt{N})$. For $N = 10^{14}$, $\sqrt{N} = 10^7$.
• Space: $O(\sqrt{N})$ for the $J_2$ array and the prime sieve.

Answer

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

const long long MOD = 1000000007;
const long long N = 100000000000000LL; // 10^14
const int L = 10000000; // sqrt(N) = 10^7

long long j2[L + 1]; // Jordan totient J_2
bool is_composite[L + 1];
vector<int> primes;

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() {
    // Initialize J_2(e) = e^2 mod MOD
    for (int i = 1; i <= L; i++) {
        j2[i] = (1LL * i % MOD) * (i % MOD) % MOD;
    }

    // Sieve: for each prime p, multiply J_2(m) by (1 - 1/p^2) = (p^2-1)/p^2
    // We need modular inverse of p^2
    for (int p = 2; p <= L; p++) {
        if (!is_composite[p]) {
            primes.push_back(p);
            long long p2 = (1LL * p * p) % MOD;
            long long factor = (p2 - 1 + MOD) % MOD * power(p2, MOD - 2, MOD) % MOD;
            for (long long m = p; m <= L; m += p) {
                is_composite[m] = (m != p);
                j2[m] = j2[m] * factor % MOD;
            }
        } else {
            // Mark composites for sieve
        }
    }

    // Actually we need a proper sieve. Let me redo using smallest prime factor.
    // Reset and use proper Euler sieve approach
    // The above loop already handles it: for each prime p, we iterate through
    // all multiples and multiply by (p^2-1)/p^2. This correctly computes J_2
    // since J_2 is multiplicative and J_2(n) = n^2 * prod_{p|n} (1 - 1/p^2).

    // Compute S(N) = sum_{e=1}^{L} J_2(e) * floor(N/e^2) mod MOD
    long long ans = 0;
    for (int e = 1; e <= L; e++) {
        long long q = (N / ((long long)e * e)) % MOD;
        ans = (ans + j2[e] * q) % MOD;
    }

    cout << ans << endl;
    return 0;
}

"""
Project Euler Problem 745: Sum of Squares II

For a positive integer n, g(n) = max perfect square dividing n.
S(N) = sum of g(n) for n=1..N.
Find S(10^14) mod 10^9+7.

Uses Jordan totient function J_2 sieve approach.
"""

def solve():
    MOD = 1000000007
    N = 10**14
    L = int(N**0.5)  # 10^7

    # Sieve J_2 values
    # J_2(n) = n^2 * prod_{p|n} (1 - 1/p^2)
    # We work modulo MOD

    j2 = [0] * (L + 1)
    for i in range(1, L + 1):
        j2[i] = (i * i) % MOD

    # Sieve of Eratosthenes style
    is_prime = bytearray(b'\x01') * (L + 1)
    is_prime[0] = is_prime[1] = 0

    for p in range(2, L + 1):
        if is_prime[p]:
            p2 = (p * p) % MOD
            inv_p2 = pow(p2, MOD - 2, MOD)
            factor = (p2 - 1) * inv_p2 % MOD
            for m in range(p, L + 1, p):
                if m != p:
                    is_prime[m] = 0
                j2[m] = j2[m] * factor % MOD

    # Compute S(N)
    ans = 0
    for e in range(1, L + 1):
        q = (N // (e * e)) % MOD
        ans = (ans + j2[e] * q) % MOD

    print(ans)

if __name__ == "__main__":
    solve()