Problem statement

Problem Statement

This archive keeps the full statement, math, and original media on the page.

For a positive integer construct a graph using all the divisors of as the vertices. An edge is drawn between and if is divisible by and is prime, and is given weight , where is the Euler totient function.
Define to be the total weight of this graph.
The example below shows that

Let . You are given and .

Find . Give your answer modulo .

Problem 931 content is attributed to Project Euler and included here under the CC BY-NC-SA 4.0 license.

Open official problem License details

Problem 931: Sum of Radical Expressions

Mathematical Foundation

Definition. For $n \geq 1$, the radical of $n$ is $\operatorname{rad}(n) = \prod_{p \mid n} p$, with $\operatorname{rad}(1) = 1$.

Lemma 1. The radical function $\operatorname{rad}$ is multiplicative: if $\gcd(a, b) = 1$, then $\operatorname{rad}(ab) = \operatorname{rad}(a) \cdot \operatorname{rad}(b)$.

Proof. Let $a = \prod p_i^{a_i}$ and $b = \prod q_j^{b_j}$ where $\{p_i\}$ and $\{q_j\}$ are disjoint (since $\gcd(a,b) = 1$). Then $ab = \prod p_i^{a_i} \prod q_j^{b_j}$ and the set of distinct prime factors of $ab$ is $\{p_i\} \cup \{q_j\}$. Hence $\operatorname{rad}(ab) = \prod p_i \cdot \prod q_j = \operatorname{rad}(a) \cdot \operatorname{rad}(b)$. $\square$

Theorem 1. The function $R(n) = \sum_{d \mid n} \operatorname{rad}(d)$ is multiplicative, and for a prime power $p^a$ with $a \geq 1$:

Proof. Since $R = \operatorname{rad} * \mathbf{1}$ is the Dirichlet convolution of two multiplicative functions, $R$ is multiplicative (a standard result in multiplicative number theory). For the prime power formula, the divisors of $p^a$ are $1, p, p^2, \ldots, p^a$. We have $\operatorname{rad}(1) = 1$ and $\operatorname{rad}(p^k) = p$ for all $k \geq 1$. Therefore:

$\square$

Theorem 2 (Closed-form product). For $n = p_1^{a_1} \cdots p_r^{a_r}$:

Proof. Immediate from Theorem 1 and the multiplicativity of $R$. $\square$

Theorem 3 (Summation by Dirichlet interchange). For any $N \geq 1$:

Proof. We have:

where the second equality follows by swapping the order of summation (each pair $(d, n)$ with $d \mid n$ and $1 \leq n \leq N$ is counted once on each side). $\square$

Editorial

Alternative (direct product approach using Theorem 2):*. We linear sieve to compute rad(d) for all d <= N. We then iterate over p from 2 to N. Finally, compute the sum via Dirichlet interchange (Theorem 3).

Pseudocode

Linear sieve to compute rad(d) for all d <= N
for p from 2 to N
Compute the sum via Dirichlet interchange (Theorem 3)
for d from 1 to N
Sieve smallest prime factor spf[n] for n <= N
for p from 2 to N
For each n, compute R(n) = prod(1 + a_i * p_i)
for n from 1 to N

Complexity Analysis

• Time: $O(N \log \log N)$ for the sieve; $O(N)$ for the summation in the Dirichlet interchange approach. The direct factorization approach takes $O(N \log N)$ in the worst case (due to repeated division), but $O(N)$ amortized. Overall: $O(N \log \log N)$.
• Space: $O(N)$ for storing $\operatorname{rad}(d)$ or the smallest prime factor array.

Answer

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

/*
 * Problem 931: Sum of Radical Expressions
 *
 * R(n) = sum_{d|n} rad(d). Find sum_{n=1}^{10^7} R(n) mod 10^9+7.
 *
 * rad(n) = product of distinct prime factors.
 * R is multiplicative: R(p^a) = 1 + a*p.
 * R(n) = prod_{p^a || n} (1 + a*p).
 *
 * Alternative: sum R(n) = sum_{d=1}^N rad(d) * floor(N/d).
 *
 * Algorithm: sieve to compute rad(d) for all d, then compute sum.
 */

int main() {
    const int N = 10000000;
    const long long MOD = 1000000007;
    
    // Compute rad(n) via sieve
    vector<long long> rad(N + 1, 1);
    for (int p = 2; p <= N; p++) {
        if (rad[p] == 1) { // p is prime
            for (int m = p; m <= N; m += p) {
                rad[m] *= p;
            }
        }
    }
    
    // sum R(n) = sum_{d=1}^N rad(d) * floor(N/d)
    long long total = 0;
    for (int d = 1; d <= N; d++) {
        total = (total + rad[d] % MOD * (N / d) % MOD) % MOD;
    }
    cout << total << endl;

    // Verify: R(6) = rad(1)+rad(2)+rad(3)+rad(6) = 1+2+3+6 = 12
    assert(rad[1] == 1 && rad[2] == 2 && rad[3] == 3 && rad[6] == 6);

    return 0;
}

"""
Problem 931: Sum of Radical Expressions

For rad(d) = product of distinct prime factors of d, define
R(n) = sum_{d | n} rad(d). Compute sum_{n=1}^{N} R(n) mod (10^9 + 7).

Key insight:
  sum_{n=1}^{N} R(n) = sum_{n=1}^{N} sum_{d|n} rad(d)
                     = sum_{d=1}^{N} rad(d) * floor(N/d)
  by swapping summation order (each d divides floor(N/d) values of n).

Methods:
  - solve: Compute rad via modified sieve, then sum rad(d)*floor(N/d) mod p.
  - solve_small: Same logic for smaller N, also returns rad array for viz.
  - rad_brute: Brute-force rad(n) by trial division for verification.

Complexity: O(N log log N) for rad sieve, O(N) for the final sum.
"""

from collections import Counter

def solve(N=10**7, MOD=10**9 + 7):
    """Compute sum_{n=1}^{N} R(n) mod MOD via rad sieve."""
    rad = [1] * (N + 1)
    is_prime = [True] * (N + 1)
    for p in range(2, N + 1):
        if is_prime[p]:
            for m in range(p, N + 1, p):
                rad[m] *= p
                if m > p:
                    is_prime[m] = False
    S = 0
    for d in range(1, N + 1):
        S = (S + rad[d] * (N // d)) % MOD
    return S

def rad_brute(n):
    """Compute rad(n) by trial division."""
    if n == 1:
        return 1
    result = 1
    temp = n
    d = 2
    while d * d <= temp:
        if temp % d == 0:
            result *= d
            while temp % d == 0:
                temp //= d
        d += 1
    if temp > 1:
        result *= temp
    return result

def solve_small(N, MOD=10**9 + 7):
    """Solve for small N, return (answer, rad_array)."""
    rad = [1] * (N + 1)
    is_prime = [True] * (N + 1)
    for p in range(2, N + 1):
        if is_prime[p]:
            for m in range(p, N + 1, p):
                rad[m] *= p
                if m > p:
                    is_prime[m] = False
    S = 0
    for d in range(1, N + 1):
        S = (S + rad[d] * (N // d)) % MOD
    return S, rad

# Verification
# Check rad values
assert rad_brute(1) == 1
assert rad_brute(12) == 6  # 12 = 2^2 * 3 -> rad = 2*3 = 6
assert rad_brute(30) == 30  # 30 = 2*3*5 -> rad = 30
assert rad_brute(8) == 2   # 8 = 2^3 -> rad = 2

# Cross-check sieve rad vs brute force
_, rad_sieve = solve_small(100)
for n in range(1, 101):
    assert rad_sieve[n] == rad_brute(n), f"rad({n}): sieve={rad_sieve[n]}, brute={rad_brute(n)}"

# Verify small answer: sum R(n) for n=1..10
# R(1)=1, R(2)=3, R(3)=4, R(4)=5, R(5)=6,
# R(6)=12, R(7)=8, R(8)=7, R(9)=7, R(10)=18
# sum = 1+3+4+5+6+12+8+7+7+18 = 71
assert solve_small(10)[0] == 71, f"Expected 71, got {solve_small(10)[0]}"

# Compute answer
ans_small, rad_arr = solve_small(10**5)
print(f"Answer for N=10^5: {ans_small}")