Problem 931: Totient Graph

Project Euler Problem 931 Solution

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

Problem Summary For each positive integer $n$, consider the graph whose vertices are the positive divisors of $n$. There is an upward edge from $b$ to $a$ exactly when $a/b$ is prime, so each edge multiplies a divisor by one prime factor and nothing else. The quantity attached to $n$ is the total edge weight $$T(n)=\sum_{\substack{b\mid n,\ a\mid n\\ a/b\text{ prime}}}\bigl(\varphi(a)-\varphi(b)\bigr).$$ The problem asks for the prefix sum $$S(N)=\sum_{n=1}^{N} T(n)\pmod{715827883},\qquad N=10^{12}.$$ A direct computation would have to build a divisor graph for every $n\le N$, which is far too expensive. The implementations succeed because the graph sum collapses to a clean prime-power formula, and the global prefix can then be reorganized into prime-power loops plus a prime-prefix query structure. Mathematical Approach The key is to understand one fixed $n$ exactly, and only then sum over all $n\le N$. The divisor graph splits by prime direction Write $$n=\prod_i p_i^{e_i}.$$ Every allowed edge in the graph is obtained by choosing one prime $p\mid n$ and moving from a divisor $d$ to $pd$, provided $pd\mid n$. So the whole graph decomposes into independent families of edges, one family for each prime dividing $n$....

Detailed mathematical approach

Problem Summary

For each positive integer $n$, consider the graph whose vertices are the positive divisors of $n$. There is an upward edge from $b$ to $a$ exactly when $a/b$ is prime, so each edge multiplies a divisor by one prime factor and nothing else. The quantity attached to $n$ is the total edge weight

$$T(n)=\sum_{\substack{b\mid n,\ a\mid n\\ a/b\text{ prime}}}\bigl(\varphi(a)-\varphi(b)\bigr).$$

The problem asks for the prefix sum

$$S(N)=\sum_{n=1}^{N} T(n)\pmod{715827883},\qquad N=10^{12}.$$

A direct computation would have to build a divisor graph for every $n\le N$, which is far too expensive. The implementations succeed because the graph sum collapses to a clean prime-power formula, and the global prefix can then be reorganized into prime-power loops plus a prime-prefix query structure.

Mathematical Approach

The key is to understand one fixed $n$ exactly, and only then sum over all $n\le N$.

The divisor graph splits by prime direction

Write

$$n=\prod_i p_i^{e_i}.$$

Every allowed edge in the graph is obtained by choosing one prime $p\mid n$ and moving from a divisor $d$ to $pd$, provided $pd\mid n$. So the whole graph decomposes into independent families of edges, one family for each prime dividing $n$.

Fix one prime $p$ with exact exponent $p^e\parallel n$, and write

$$n=p^e m,\qquad p\nmid m.$$

Then every edge whose ratio is $p$ has the form

$$p^r d \longrightarrow p^{r+1} d,$$

with $0\le r<e$ and $d\mid m$. So to find the total contribution of the $p$-direction, we only need to sum the totient differences on these edges.

Summing one prime direction exactly

Because $p\nmid d$, Euler's totient function satisfies

$$\varphi(p^r d)=\varphi(p^r)\varphi(d).$$

There are two cases.

If $r=0$, then

$$\varphi(pd)-\varphi(d)=(p-1)\varphi(d)-\varphi(d)=(p-2)\varphi(d).$$

If $r\ge 1$, then

$$\varphi(p^{r+1}d)-\varphi(p^r d)=p^r(p-1)\varphi(d)-p^{r-1}(p-1)\varphi(d)=p^{r-1}(p-1)^2\varphi(d).$$

Therefore the full contribution of all $p$-edges is

$$\sum_{d\mid m}\varphi(d)\left[(p-2)+\sum_{r=1}^{e-1}p^{r-1}(p-1)^2\right].$$

The bracket simplifies to

$$\begin{aligned} (p-2)+(p-1)^2\sum_{r=1}^{e-1}p^{r-1} &=(p-2)+(p-1)^2\frac{p^{e-1}-1}{p-1} \\ &=(p-2)+(p-1)(p^{e-1}-1) \\ &=p^e-p^{e-1}-1. \end{aligned}$$

Now use the classical identity

$$\sum_{d\mid m}\varphi(d)=m.$$

So the entire $p$-direction contributes

$$m\bigl(p^e-p^{e-1}-1\bigr)=n-\frac{n}{p}-\frac{n}{p^e}.$$

Summing over the distinct prime powers of $n$ gives the closed form used by the code:

$$\boxed{T(n)=\sum_{p^e\parallel n}\left(n-\frac{n}{p}-\frac{n}{p^e}\right).}$$

Worked example: $n=12$

The divisors of $12$ are $1,2,3,4,6,12$. The prime-ratio edges are

$$1\to2,\quad 1\to3,\quad 2\to4,\quad 2\to6,\quad 3\to6,\quad 4\to12,\quad 6\to12.$$

Their weights are

$$0,\ 1,\ 1,\ 1,\ 0,\ 2,\ 2,$$

because

$$\varphi(1)=1,\ \varphi(2)=1,\ \varphi(3)=2,\ \varphi(4)=2,\ \varphi(6)=2,\ \varphi(12)=4.$$

So $T(12)=7$.

The closed form gives exactly the same result. Since $12=2^2\cdot 3$,

$$T(12)=\left(12-\frac{12}{2}-\frac{12}{2^2}\right)+\left(12-\frac{12}{3}-\frac{12}{3}\right)=3+4=7.$$

This example shows what the formula is really doing: it adds the total contribution of each prime direction in the divisor graph.

Reordering the prefix sum by exact prime powers

To compute $S(N)$, rewrite every $n$ contributing to the term for $p^e$ as

$$n=p^e m,\qquad p\nmid m.$$

For such an $n$, the corresponding summand equals

$$n-\frac{n}{p}-\frac{n}{p^e}=m\bigl(p^e-p^{e-1}-1\bigr).$$

Hence

$$S(N)=\sum_{p^e\le N}\bigl(p^e-p^{e-1}-1\bigr)\sum_{\substack{m\le N/p^e\\ p\nmid m}} m.$$

The inner sum is simple once the multiples of $p$ are removed. Define the triangular sum

$$\operatorname{Tri}(x)=\frac{x(x+1)}{2}.$$

Then

$$\sum_{\substack{m\le x\\ p\nmid m}} m=\operatorname{Tri}(x)-p\,\operatorname{Tri}\!\left(\left\lfloor\frac{x}{p}\right\rfloor\right),$$

because the excluded multiples are $p,2p,\dots,p\lfloor x/p\rfloor$.

So the whole prefix sum becomes

$$\boxed{S(N)=\sum_{p^e\le N}\bigl(p^e-p^{e-1}-1\bigr)\left[\operatorname{Tri}\!\left(\left\lfloor\frac{N}{p^e}\right\rfloor\right)-p\,\operatorname{Tri}\!\left(\left\lfloor\frac{N}{p^{e+1}}\right\rfloor\right)\right].}$$

Why primes above $\sqrt N$ can be batched together

The implementations split the summation at $\sqrt N$.

If $p\le \sqrt N$, the code can iterate over $p,p^2,p^3,\dots$ directly.

If $p>\sqrt N$, then $p^2>N$, so only the exponent $e=1$ can occur. In that case

$$x=\left\lfloor\frac{N}{p}\right\rfloor<p,$$

which means there is no positive multiple of $p$ inside $\{1,\dots,x\}$. Therefore

$$\sum_{\substack{m\le x\\ p\nmid m}}m=\operatorname{Tri}(x),$$

and the contribution of a large prime is simply

$$\bigl(p-2\bigr)\operatorname{Tri}\!\left(\left\lfloor\frac{N}{p}\right\rfloor\right).$$

Now group large primes by the constant quotient

$$q=\left\lfloor\frac{N}{p}\right\rfloor.$$

All primes in the interval

$$L=\left\lfloor\frac{N}{q+1}\right\rfloor+1,\qquad R=\left\lfloor\frac{N}{q}\right\rfloor$$

share this same $q$, so their total contribution is

$$\operatorname{Tri}(q)\sum_{L\le p\le R}(p-2).$$

This reduces the large-prime side to interval queries for

$$\pi(x)=\#\{p\le x\},\qquad P(x)=\sum_{p\le x}p.$$

Once those two prime-prefix functions are available, every large-prime block can be evaluated as

$$\operatorname{Tri}(q)\Bigl(P(R)-P(L-1)-2\bigl(\pi(R)-\pi(L-1)\bigr)\Bigr).$$

How the Code Works

Validation comes first

The reference implementation does not assume the closed form blindly. It first builds the divisor graph directly for small $n$, evaluates the edge sum from the definition, and checks that this agrees with the prime-power formula for all $n$ in a small validation range. It also verifies several prefix values, including $S(10)=26$ and $S(100)=5282$, before moving on to the full $N=10^{12}$ computation.

Compressed prime-prefix preprocessing

To answer many prime-count and prime-sum queries quickly, the implementation stores the distinct floor-quotients

$$\left\lfloor\frac{N}{1}\right\rfloor,\left\lfloor\frac{N}{2}\right\rfloor,\left\lfloor\frac{N}{3}\right\rfloor,\dots$$

instead of every integer up to $N$. There are only $O(\sqrt N)$ distinct values of this form. On that compressed coordinate set it initializes counts and sums of all integers at least $2$, then removes composite contributions prime by prime. After that preprocessing, the code can query $\pi(x)$ and $P(x)$ for every relevant $x$ in constant time.

Accumulating the answer

The final sum is assembled in two phases. First, every prime $p\le\sqrt N$ is traversed through its powers $p,p^2,p^3,\dots$, and the exact prime-power formula is added. Second, primes above $\sqrt N$ are processed in quotient blocks with constant $q=\lfloor N/p\rfloor$, using prime counts and prime sums on intervals.

All arithmetic is reduced modulo $715827883$. The C++ implementation uses 128-bit intermediates before taking residues, the Python implementation relies on arbitrary-precision integers, and the Java entry point delegates the heavy computation to the compiled C++ solver and returns its result.

Complexity Analysis

The distinct values of $\lfloor N/t\rfloor$ form an $O(\sqrt N)$-sized state space, and both the large-prime grouping and the query structure are built around that compression. Memory usage is therefore $O(\sqrt N)$.

Time is dominated by the compressed prime-prefix sieve together with the sweep over primes up to $\sqrt N$ and their powers. The important point is qualitative: the algorithm is far below linear in $N$ and avoids any attempt to enumerate divisor graphs for all $n\le 10^{12}$. That is what makes the computation practical.

Footnotes and References

Problem page: https://projecteuler.net/problem=931
Euler's totient function: Wikipedia - Euler's totient function
Divisor: Wikipedia - Divisor
Prime-counting function: Wikipedia - Prime-counting function
Dirichlet hyperbola method: Wikipedia - Dirichlet hyperbola method

Problem 931 source code

C++

#include <algorithm>
#include <cassert>
#include <cmath>
#include <cstdint>
#include <iostream>
#include <vector>

namespace {

using i64 = std::int64_t;
using u64 = std::uint64_t;
using i128 = __int128_t;

constexpr u64 kMod = 715'827'883ULL;

i128 tri(u64 x) {
    return static_cast<i128>(x) * static_cast<i128>(x + 1) / 2;
}

bool is_prime_small(int x) {
    if (x < 2) {
        return false;
    }
    for (int d = 2; static_cast<i64>(d) * d <= x; ++d) {
        if (x % d == 0) {
            return false;
        }
    }
    return true;
}

int phi_small(int n) {
    int result = n;
    int x = n;
    for (int p = 2; static_cast<i64>(p) * p <= x; ++p) {
        if (x % p != 0) {
            continue;
        }
        while (x % p == 0) {
            x /= p;
        }
        result -= result / p;
    }
    if (x > 1) {
        result -= result / x;
    }
    return result;
}

i64 t_bruteforce_graph(int n) {
    std::vector<int> divisors;
    for (int d = 1; static_cast<i64>(d) * d <= n; ++d) {
        if (n % d != 0) {
            continue;
        }
        divisors.push_back(d);
        if (d * d != n) {
            divisors.push_back(n / d);
        }
    }

    i64 total = 0;
    for (int a : divisors) {
        for (int b : divisors) {
            if (a % b != 0) {
                continue;
            }
            const int q = a / b;
            if (!is_prime_small(q)) {
                continue;
            }
            total += static_cast<i64>(phi_small(a) - phi_small(b));
        }
    }
    return total;
}

i64 t_formula_small(int n) {
    int x = n;
    i64 total = 0;

    for (int p = 2; static_cast<i64>(p) * p <= x; ++p) {
        if (x % p != 0) {
            continue;
        }
        i64 pe = 1;
        while (x % p == 0) {
            x /= p;
            pe *= p;
        }
        total += static_cast<i64>(n) - n / p - n / pe;
    }
    if (x > 1) {
        total += static_cast<i64>(n) - 2LL * (n / x);
    }

    return total;
}

struct PrimePrefix {
    u64 n = 0;
    u64 root = 0;

    std::vector<u64> values;
    std::vector<int> idx_small;
    std::vector<int> idx_large;

    std::vector<i64> g_count;
    std::vector<i128> g_sum;

    std::vector<int> primes;

    int id(u64 x) const {
        if (x <= root) {
            return idx_small[x];
        }
        return idx_large[n / x];
    }

    explicit PrimePrefix(u64 limit) : n(limit), root(static_cast<u64>(std::sqrt(static_cast<long double>(limit)))) {
        for (u64 l = 1, r; l <= n; l = r + 1) {
            const u64 v = n / l;
            r = n / v;
            values.push_back(v);
        }

        idx_small.assign(root + 1, -1);
        idx_large.assign(root + 1, -1);

        const int m = static_cast<int>(values.size());
        g_count.resize(m);
        g_sum.resize(m);

        for (int i = 0; i < m; ++i) {
            const u64 v = values[i];
            if (v <= root) {
                idx_small[v] = i;
            } else {
                idx_large[n / v] = i;
            }
            g_count[i] = static_cast<i64>(v) - 1;
            g_sum[i] = tri(v) - 1;
        }

        std::vector<int> lp(root + 1, 0);
        for (u64 i = 2; i <= root; ++i) {
            if (lp[i] == 0) {
                lp[i] = static_cast<int>(i);
                primes.push_back(static_cast<int>(i));
            }
            for (int p : primes) {
                const u64 v = i * static_cast<u64>(p);
                if (v > root || p > lp[i]) {
                    break;
                }
                lp[v] = p;
            }
        }

        for (int p : primes) {
            const u64 p64 = static_cast<u64>(p);
            const u64 p2 = p64 * p64;
            if (p2 > n) {
                break;
            }

            const int idx_pm1 = id(p64 - 1);

            int i = 0;
            while (i < m && values[i] >= p2) {
                const int j = id(values[i] / p64);
                g_count[i] -= g_count[j] - g_count[idx_pm1];
                g_sum[i] -= static_cast<i128>(p64) * (g_sum[j] - g_sum[idx_pm1]);
                ++i;
            }
        }
    }

    i64 prime_count(u64 x) const {
        if (x < 2) {
            return 0;
        }
        return g_count[id(x)];
    }

    i128 prime_sum(u64 x) const {
        if (x < 2) {
            return 0;
        }
        return g_sum[id(x)];
    }
};

u64 solve(u64 N) {
    PrimePrefix prefix(N);
    const u64 root = prefix.root;

    u64 answer = 0;

    auto add_mod_i128 = [&](i128 x) {
        x %= static_cast<i128>(kMod);
        if (x < 0) {
            x += static_cast<i128>(kMod);
        }
        answer += static_cast<u64>(x);
        if (answer >= kMod) {
            answer -= kMod;
        }
    };

    for (int p : prefix.primes) {
        const u64 p64 = static_cast<u64>(p);
        u64 pe = p64;

        while (pe <= N) {
            const u64 x = N / pe;
            const i128 s = tri(x) - static_cast<i128>(p64) * tri(x / p64);
            const i128 coef = static_cast<i128>(pe) - static_cast<i128>(pe / p64) - 1;
            add_mod_i128(coef * s);

            if (pe > N / p64) {
                break;
            }
            pe *= p64;
        }
    }

    const u64 max_q = N / (root + 1);
    for (u64 q = 1; q <= max_q; ++q) {
        u64 L = N / (q + 1) + 1;
        u64 R = N / q;

        if (R <= root) {
            break;
        }
        if (L <= root) {
            L = root + 1;
        }
        if (L > R) {
            continue;
        }

        const i64 cnt = prefix.prime_count(R) - prefix.prime_count(L - 1);
        if (cnt == 0) {
            continue;
        }

        const i128 sp = prefix.prime_sum(R) - prefix.prime_sum(L - 1);
        const i128 group = tri(q) * (sp - 2 * static_cast<i128>(cnt));
        add_mod_i128(group);
    }

    return answer;
}

u64 brute_prefix_formula(int N) {
    i128 s = 0;
    for (int n = 1; n <= N; ++n) {
        s += t_formula_small(n);
    }
    s %= static_cast<i128>(kMod);
    if (s < 0) {
        s += static_cast<i128>(kMod);
    }
    return static_cast<u64>(s);
}

void run_validations() {
    assert(t_bruteforce_graph(45) == 52);
    for (int n = 1; n <= 200; ++n) {
        assert(t_formula_small(n) == t_bruteforce_graph(n));
    }

    assert(solve(10) == 26);
    assert(solve(100) == 5282);
    assert(solve(1000) == brute_prefix_formula(1000));
}

}  // namespace

int main() {
    run_validations();
    constexpr u64 kN = 1'000'000'000'000ULL;
    std::cout << solve(kN) << '\n';
    return 0;
}

Python

import math

def solve():
    MOD = 715827883; N = 1000000000000

    def tri(x): return x*(x+1)//2

    root = int(math.isqrt(N))
    while (root+1)*(root+1) <= N: root += 1
    while root*root > N: root -= 1

    # Min25 prime prefix sieve
    vals = []
    l = 1
    while l <= N:
        v = N // l; r = N // v; vals.append(v); l = r + 1

    idx_small = [0]*(root+1); idx_large = [0]*(root+1)
    m = len(vals)
    g_count = [0]*m; g_sum = [0]*m

    for i in range(m):
        v = vals[i]
        if v <= root: idx_small[v] = i
        else: idx_large[N//v] = i
        g_count[i] = v - 1; g_sum[i] = tri(v) - 1

    lp = [0]*(root+1); primes = []
    for i in range(2, root+1):
        if lp[i] == 0: lp[i] = i; primes.append(i)
        for p in primes:
            v = i*p
            if v > root or p > lp[i]: break
            lp[v] = p

    def vid(x):
        return idx_small[x] if x <= root else idx_large[N//x]

    for p in primes:
        p2 = p*p
        if p2 > N: break
        ipm1 = vid(p-1)
        i = 0
        while i < m and vals[i] >= p2:
            j = vid(vals[i]//p)
            g_count[i] -= g_count[j] - g_count[ipm1]
            g_sum[i] -= p * (g_sum[j] - g_sum[ipm1])
            i += 1

    answer = 0
    for p in primes:
        pe = p
        while pe <= N:
            x = N // pe
            s = tri(x) - p * tri(x // p)
            coef = pe - pe // p - 1
            v = coef * s % MOD
            answer = (answer + v) % MOD

            if pe > N // p: break
            pe *= p

    max_q = N // (root + 1)
    for q in range(1, max_q + 1):
        L = N // (q + 1) + 1; R = N // q
        if R <= root: break
        if L <= root: L = root + 1
        if L > R: continue
        cnt = g_count[vid(R)] - g_count[vid(L-1)]
        if cnt == 0: continue
        sp = g_sum[vid(R)] - g_sum[vid(L-1)]
        group = tri(q) * (sp - 2 * cnt)
        answer = (answer + group) % MOD

    return str(answer % MOD)

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

Java

import java.nio.file.*;
import java.util.*;
import java.util.regex.*;

public class Euler931 {
    private static final Pattern ANSWER_RE = Pattern.compile("answer\\s*:\\s*(.+)$", Pattern.CASE_INSENSITIVE);
    private static final Pattern EQUAL_RE = Pattern.compile("=\\s*(.+)$");

    private static String parseOutput(String stdout) {
        String[] lines = stdout.split("\\R");
        List<String> nonEmpty = new ArrayList<>();
        for (String line : lines) {
            String t = line.trim();
            if (!t.isEmpty()) {
                nonEmpty.add(t);
            }
        }
        if (nonEmpty.isEmpty()) {
            return "";
        }

        List<String> answers = new ArrayList<>();
        List<String> equals = new ArrayList<>();
        for (String line : nonEmpty) {
            Matcher m1 = ANSWER_RE.matcher(line);
            if (m1.find()) {
                answers.add(m1.group(1).trim());
            }
            Matcher m2 = EQUAL_RE.matcher(line);
            if (m2.find()) {
                equals.add(m2.group(1).trim());
            }
        }

        if (!answers.isEmpty()) {
            return answers.get(answers.size() - 1);
        }
        if (!equals.isEmpty()) {
            return equals.get(equals.size() - 1);
        }
        return nonEmpty.get(nonEmpty.size() - 1);
    }

    private static String pickCompiler() throws Exception {
        for (String compiler : List.of("clang++", "g++")) {
            Process probe = new ProcessBuilder("bash", "-lc", "command -v " + compiler)
                    .redirectErrorStream(true)
                    .start();
            String out = new String(probe.getInputStream().readAllBytes());
            int rc = probe.waitFor();
            if (rc == 0 && !out.trim().isEmpty()) {
                return compiler;
            }
        }
        throw new RuntimeException("No C++ compiler found (clang++/g++).");
    }

    private static Path cppSource(Path root) {
        return root.resolve("solutionsCpp").resolve("Euler931.cpp");
    }

    private static boolean shouldSkipCheckpoints(Path root) {
        Path src = cppSource(root);
        try {
            String text = Files.readString(src);
            return text.contains("--skip-checkpoints");
        } catch (Exception ex) {
            return false;
        }
    }

    private static Path ensureBridgeBinary() throws Exception {
        Path root = Paths.get(System.getProperty("user.dir"));
        Path src = cppSource(root);
        Path bin = root.resolve("solutionsCpp").resolve(".euler931_java_bridge");

        boolean rebuild = Files.notExists(bin)
                || Files.getLastModifiedTime(src).compareTo(Files.getLastModifiedTime(bin)) > 0;

        if (rebuild) {
            String compiler = pickCompiler();
            Process compile = new ProcessBuilder(
                    compiler,
                    "-std=c++17",
                    "-O2",
                    src.toString(),
                    "-o",
                    bin.toString())
                    .inheritIO()
                    .start();
            if (compile.waitFor() != 0) {
                throw new RuntimeException("Failed to compile Euler931 C++ bridge.");
            }
        }

        return bin;
    }

    private static String runBridge(Path bin, Path root, Path srcDir) throws Exception {
        List<String> cmd = new ArrayList<>();
        cmd.add(bin.toString());
        if (shouldSkipCheckpoints(root)) {
            cmd.add("--skip-checkpoints");
        }

        Process first = new ProcessBuilder(cmd)
                .directory(root.toFile())
                .redirectErrorStream(true)
                .start();
        String out = new String(first.getInputStream().readAllBytes());
        int rc = first.waitFor();
        if (rc == 0) {
            return out;
        }

        Process second = new ProcessBuilder(cmd)
                .directory(srcDir.toFile())
                .redirectErrorStream(true)
                .start();
        String out2 = new String(second.getInputStream().readAllBytes());
        int rc2 = second.waitFor();
        if (rc2 == 0) {
            return out2;
        }

        throw new RuntimeException("Euler931 C++ bridge failed.\n" + out + "\n" + out2);
    }

    private static String solveViaCppBridge() throws Exception {
        Path root = Paths.get(System.getProperty("user.dir"));
        Path src = cppSource(root);
        Path bin = ensureBridgeBinary();
        String out = runBridge(bin, root, src.getParent());
        String parsed = parseOutput(out);
        if (parsed.isEmpty()) {
            throw new RuntimeException("Euler931 C++ bridge produced empty output.");
        }
        return parsed;
    }

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