Problem 625: Gcd Sum

View on Project Euler

Project Euler Problem 625 Solution

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

Problem Summary Define $$G(N)=\sum_{j=1}^{N}\sum_{i=1}^{j}\gcd(i,j).$$ The problem asks for \(G(10^{11}) \bmod 998244353\). A direct double loop is impossible at that scale, so the solution rewrites the gcd contribution as a divisor sum, introduces the summatory totient function, and then evaluates the remaining expressions by floor-division blocks. Mathematical Approach It is convenient to isolate the inner sum $$T(j)=\sum_{i=1}^{j}\gcd(i,j),$$ and also to define the totient prefix $$\Phi(x)=\sum_{k=1}^{x}\varphi(k).$$ The entire method is built around turning \(T(j)\) into something that depends on divisors of \(j\), and then exploiting the fact that floor quotients stay constant on long intervals. Step 1: Rewrite the inner gcd sum as a divisor sum Fix \(j\). Group the indices \(i\) by the value \(d=\gcd(i,j)\). If \(d\) is fixed, then we can write $$i=d\,a,\qquad j=d\,b,\qquad \gcd(a,b)=1.$$ Here \(b=j/d\), and the admissible values of \(a\) are exactly the integers \(1\le a\le b\) that are coprime to \(b\). Their number is \(\varphi(b)\). So each divisor \(d\mid j\) contributes the value \(d\) exactly \(\varphi(j/d)\) times, which gives $$T(j)=\sum_{d\mid j} d\,\varphi\!\left(\frac{j}{d}\right).$$ This identity removes the explicit gcd from the inner loop and replaces it with a clean divisor formula....

Detailed mathematical approach

Problem Summary

Define

$$G(N)=\sum_{j=1}^{N}\sum_{i=1}^{j}\gcd(i,j).$$

The problem asks for \(G(10^{11}) \bmod 998244353\). A direct double loop is impossible at that scale, so the solution rewrites the gcd contribution as a divisor sum, introduces the summatory totient function, and then evaluates the remaining expressions by floor-division blocks.

Mathematical Approach

It is convenient to isolate the inner sum

$$T(j)=\sum_{i=1}^{j}\gcd(i,j),$$

and also to define the totient prefix

$$\Phi(x)=\sum_{k=1}^{x}\varphi(k).$$

The entire method is built around turning \(T(j)\) into something that depends on divisors of \(j\), and then exploiting the fact that floor quotients stay constant on long intervals.

Step 1: Rewrite the inner gcd sum as a divisor sum

Fix \(j\). Group the indices \(i\) by the value \(d=\gcd(i,j)\). If \(d\) is fixed, then we can write

$$i=d\,a,\qquad j=d\,b,\qquad \gcd(a,b)=1.$$

Here \(b=j/d\), and the admissible values of \(a\) are exactly the integers \(1\le a\le b\) that are coprime to \(b\). Their number is \(\varphi(b)\).

So each divisor \(d\mid j\) contributes the value \(d\) exactly \(\varphi(j/d)\) times, which gives

$$T(j)=\sum_{d\mid j} d\,\varphi\!\left(\frac{j}{d}\right).$$

This identity removes the explicit gcd from the inner loop and replaces it with a clean divisor formula.

Step 2: Reverse the order of summation

Now sum \(T(j)\) over all \(j\le N\):

$$G(N)=\sum_{j=1}^{N}\sum_{d\mid j} d\,\varphi\!\left(\frac{j}{d}\right).$$

Write \(j=d\,k\). Every pair \((d,k)\) with \(d\,k\le N\) appears exactly once, so

$$G(N)=\sum_{d=1}^{N} d\sum_{k\le N/d}\varphi(k)=\sum_{d=1}^{N} d\,\Phi\!\left(\left\lfloor\frac{N}{d}\right\rfloor\right).$$

At this point the whole problem is reduced to answering many queries for the prefix sum \(\Phi(x)\).

Step 3: Derive a recurrence for the totient prefix

The classical identity

$$\sum_{d\mid n}\varphi(d)=n$$

holds for every positive integer \(n\). Summing it over \(1\le n\le x\) gives

$$\sum_{n=1}^{x} n=\sum_{n=1}^{x}\sum_{d\mid n}\varphi(d).$$

Swap the order of summation. A fixed \(d\) divides exactly \(\left\lfloor x/d\right\rfloor\) integers up to \(x\), so

$$\frac{x(x+1)}{2}=\sum_{d=1}^{x}\varphi(d)\left\lfloor\frac{x}{d}\right\rfloor.$$

Rewrite the right-hand side by grouping equal quotients. This yields

$$\frac{x(x+1)}{2}=\sum_{m=1}^{x}\Phi\!\left(\left\lfloor\frac{x}{m}\right\rfloor\right).$$

Isolating the \(m=1\) term gives the recurrence

$$\Phi(x)=\frac{x(x+1)}{2}-\sum_{m=2}^{x}\Phi\!\left(\left\lfloor\frac{x}{m}\right\rfloor\right).$$

This is exactly the large-argument formula used by the implementation.

Step 4: Compress the recurrence with floor-division blocks

The quantity \(\left\lfloor x/m\right\rfloor\) does not change at every index. If

$$q=\left\lfloor\frac{x}{\ell}\right\rfloor,\qquad r=\left\lfloor\frac{x}{q}\right\rfloor,$$

then every \(m\in[\ell,r]\) has the same quotient \(q\). Therefore the recurrence becomes

$$\Phi(x)=\frac{x(x+1)}{2}-\sum_{\text{blocks }[\ell,r]} (r-\ell+1)\,\Phi(q).$$

Instead of iterating over all \(m\), we only iterate over the distinct quotient blocks. For large \(x\), there are only \(O(\sqrt{x})\) such blocks.

The C++, Python, and Java implementations precompute \(\varphi(n)\) and its prefix values for all \(n\le 5\times 10^6\), and only use this recursive block formula beyond that cutoff.

Step 5: Apply the same block idea to the outer sum

The transformed target expression

$$G(N)=\sum_{d=1}^{N} d\,\Phi\!\left(\left\lfloor\frac{N}{d}\right\rfloor\right)$$

has the same floor structure. If \(\left\lfloor N/d\right\rfloor=q\) on a block \(d\in[\ell,r]\), then the entire block contributes

$$\left(\sum_{d=\ell}^{r} d\right)\Phi(q).$$

The arithmetic progression sum is

$$\sum_{d=\ell}^{r} d=\frac{(\ell+r)(r-\ell+1)}{2}.$$

So the final evaluation is again a loop over quotient blocks rather than a loop over all \(d\le N\).

Worked Example: \(N=10\)

The implementations verify the small checkpoint \(G(10)=122\). The block formula reproduces it directly.

First compute the needed totient prefixes:

$$\Phi(1)=1,\qquad \Phi(2)=2,\qquad \Phi(3)=4,\qquad \Phi(5)=10,\qquad \Phi(10)=32.$$

Now write

$$G(10)=\sum_{d=1}^{10} d\,\Phi\!\left(\left\lfloor\frac{10}{d}\right\rfloor\right).$$

The quotients \(\left\lfloor 10/d\right\rfloor\) are constant on the blocks

$$[1,1],\qquad [2,2],\qquad [3,3],\qquad [4,5],\qquad [6,10],$$

with quotient values \(10,5,3,2,1\), respectively. Hence

$$\begin{aligned} G(10)&=1\cdot \Phi(10)+2\cdot \Phi(5)+3\cdot \Phi(3)+(4+5)\Phi(2)+(6+7+8+9+10)\Phi(1)\\ &=1\cdot 32+2\cdot 10+3\cdot 4+9\cdot 2+40\cdot 1\\ &=32+20+12+18+40=122. \end{aligned}$$

This small case shows exactly why grouping by equal floor quotients removes the need for a linear scan.

How the Code Works

The C++, Python, and Java implementations all follow the same pipeline. First they build \(\varphi(n)\) up to \(5\times 10^6\) with a linear sieve and turn it into a prefix table, so every small \(\Phi(x)\) query becomes a constant-time lookup.

For larger arguments, the implementation evaluates \(\Phi(x)\) from the recurrence above, grouping equal values of \(\left\lfloor x/m\right\rfloor\) into blocks and caching each large result the first time it is computed. The same large quotient can appear repeatedly, so memoization is essential.

After that, the implementation computes \(G(N)\) with another quotient-block loop over \(d\). On each block it evaluates the arithmetic progression sum, multiplies it by the already known value of \(\Phi(q)\), and accumulates everything modulo \(998244353\). Because the modulus is odd, every division by \(2\) is handled as multiplication by the modular inverse of \(2\).

Complexity Analysis

Let \(L=5\times 10^6\). The totient sieve and prefix construction cost \(O(L)\) time and \(O(L)\) memory. Each large \(\Phi(x)\) query is memoized once and processed by quotient blocks, so its work is proportional to the number of distinct values of \(\left\lfloor x/k\right\rfloor\), namely \(O(\sqrt{x})\). The outer sum is also traversed by quotient blocks, so the full computation is far below \(O(N)\) and is practical for \(N=10^{11}\).

Footnotes and References

  1. Problem page: https://projecteuler.net/problem=625
  2. Greatest common divisor: Wikipedia - Greatest common divisor
  3. Euler's totient function: Wikipedia - Euler's totient function
  4. Dirichlet convolution: Wikipedia - Dirichlet convolution
  5. Linear sieve: cp-algorithms - Linear Sieve

Problem 625 source code

C++

#include <algorithm>
#include <cassert>
#include <cstdint>
#include <iostream>
#include <numeric>
#include <unordered_map>
#include <utility>
#include <vector>

using u64 = unsigned long long;
using u128 = __uint128_t;

static constexpr int MOD = 998244353;
static constexpr int INV2 = (MOD + 1) / 2;

static inline int mod_add(int a, int b) {
    int s = a + b;
    if (s >= MOD) s -= MOD;
    return s;
}

static inline int mod_sub(int a, int b) {
    int s = a - b;
    if (s < 0) s += MOD;
    return s;
}

static inline int mod_mul(u64 a, u64 b) { return (int)((u128)(a % MOD) * (b % MOD) % MOD); }

static inline int tri_mod(u64 l, u64 r) {
    u64 cnt = r - l + 1;
    int a = (int)((l + r) % MOD);
    int b = (int)(cnt % MOD);
    return mod_mul((u64)a * b, INV2);
}

struct PhiPrefix {
    int limit;
    std::vector<int> pref;
};

static PhiPrefix build_phi_prefix(int limit) {
    std::vector<int> phi(limit + 1, 0);
    std::vector<int> primes;
    std::vector<unsigned char> is_comp(limit + 1, 0);
    primes.reserve((size_t)limit / 10);

    phi[0] = 0;
    if (limit >= 1) phi[1] = 1;

    for (int i = 2; i <= limit; ++i) {
        if (!is_comp[i]) {
            primes.push_back(i);
            phi[i] = i - 1;
        }
        for (int p : primes) {
            long long v = 1LL * p * i;
            if (v > limit) break;
            is_comp[(int)v] = 1;
            if (i % p == 0) {
                phi[(int)v] = phi[i] * p;
                break;
            }
            phi[(int)v] = phi[i] * (p - 1);
        }
    }

    std::vector<int> pref(limit + 1, 0);
    for (int i = 1; i <= limit; ++i) {
        pref[i] = mod_add(pref[i - 1], phi[i] % MOD);
    }
    return PhiPrefix{limit, std::move(pref)};
}

static int sum_phi(u64 n, const PhiPrefix &pre, std::unordered_map<u64, int> &memo) {
    if (n <= (u64)pre.limit) return pre.pref[(size_t)n];
    auto it = memo.find(n);
    if (it != memo.end()) return it->second;

    int res = mod_mul(n, n + 1);
    res = mod_mul(res, INV2);

    for (u64 l = 2; l <= n;) {
        u64 q = n / l;
        u64 r = n / q;
        int cnt = (int)((r - l + 1) % MOD);
        int sub = mod_mul(cnt, (u64)sum_phi(q, pre, memo));
        res = mod_sub(res, sub);
        l = r + 1;
    }

    memo.emplace(n, res);
    return res;
}

static int G(u64 N, const PhiPrefix &pre, std::unordered_map<u64, int> &memo) {
    (void)sum_phi(N, pre, memo);
    int ans = 0;
    for (u64 l = 1; l <= N;) {
        u64 q = N / l;
        u64 r = N / q;
        int sum_d = tri_mod(l, r);
        int sphi = sum_phi(q, pre, memo);
        ans = mod_add(ans, mod_mul((u64)sum_d, (u64)sphi));
        l = r + 1;
    }
    return ans;
}

static u64 G_bruteforce(u64 N) {
    u64 ans = 0;
    for (u64 j = 1; j <= N; ++j) {
        for (u64 i = 1; i <= j; ++i) ans += std::gcd(i, j);
    }
    return ans;
}

int main() {
    const int LIMIT = 5'000'000;
    const PhiPrefix pre = build_phi_prefix(LIMIT);
    std::unordered_map<u64, int> memo;
    memo.reserve(1 << 18);

    assert(G_bruteforce(10) == 122);
    assert(G(10, pre, memo) == 122);
    std::cout << G(100'000'000'000ULL, pre, memo) << "\n";
    return 0;
}

Python

def solve():
    MOD = 998244353
    INV2 = (MOD + 1) // 2
    N = 100_000_000_000
    SIEVE_LIMIT = 5_000_000

    def mod_mul(a, b):
        return a % MOD * (b % MOD) % MOD

    def tri_mod(l, r):
        a = (l + r) % MOD
        b = (r - l + 1) % MOD
        return a * b % MOD * INV2 % MOD

    # Totient sieve
    phi = list(range(SIEVE_LIMIT + 1))
    is_comp = bytearray(SIEVE_LIMIT + 1)
    primes = []
    phi[0] = 0
    phi[1] = 1
    for i in range(2, SIEVE_LIMIT + 1):
        if not is_comp[i]:
            primes.append(i)
            phi[i] = i - 1
        for p in primes:
            if i * p > SIEVE_LIMIT:
                break
            is_comp[i * p] = 1
            if i % p == 0:
                phi[i * p] = phi[i] * p
                break
            phi[i * p] = phi[i] * (p - 1)

    pref = [0] * (SIEVE_LIMIT + 1)
    for i in range(1, SIEVE_LIMIT + 1):
        pref[i] = (pref[i-1] + phi[i]) % MOD

    memo = {}

    def sum_phi(n):
        if n <= SIEVE_LIMIT:
            return pref[n]
        if n in memo:
            return memo[n]
        res = n % MOD * ((n + 1) % MOD) % MOD * INV2 % MOD
        l = 2
        while l <= n:
            q = n // l
            r = n // q
            cnt = (r - l + 1) % MOD
            sub = cnt * sum_phi(q) % MOD
            res = (res - sub) % MOD
            l = r + 1
        memo[n] = res
        return res

    sum_phi(N)

    ans = 0
    l = 1
    while l <= N:
        q = N // l
        r = N // q
        sd = tri_mod(l, r)
        sp = sum_phi(q)
        ans = (ans + sd * sp) % MOD
        l = r + 1

    return str(ans % MOD)

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

Java

import java.util.ArrayList;
import java.util.HashMap;
import java.util.List;
import java.util.Map;

public class Euler625 {
    static final int MOD = 998244353;
    static final int INV2 = (MOD + 1) / 2;

    static int triMod(long l, long r) {
        long cnt = r - l + 1;
        long a = (l + r) % MOD;
        long b = cnt % MOD;
        return (int) ((a * b % MOD) * INV2 % MOD);
    }

    static class PhiPrefix {
        int limit;
        int[] pref;

        PhiPrefix(int limit, int[] pref) {
            this.limit = limit;
            this.pref = pref;
        }
    }

    static PhiPrefix buildPhiPrefix(int limit) {
        int[] phi = new int[limit + 1];
        List<Integer> primes = new ArrayList<>();
        byte[] isComp = new byte[limit + 1];

        if (limit >= 1)
            phi[1] = 1;

        for (int i = 2; i <= limit; i++) {
            if (isComp[i] == 0) {
                primes.add(i);
                phi[i] = i - 1;
            }
            for (int p : primes) {
                long v = (long) p * i;
                if (v > limit)
                    break;
                isComp[(int) v] = 1;
                if (i % p == 0) {
                    phi[(int) v] = phi[i] * p;
                    break;
                }
                phi[(int) v] = phi[i] * (p - 1);
            }
        }

        int[] pref = new int[limit + 1];
        for (int i = 1; i <= limit; i++) {
            pref[i] = (pref[i - 1] + phi[i]) % MOD;
        }
        return new PhiPrefix(limit, pref);
    }

    static int sumPhi(long n, PhiPrefix pre, Map<Long, Integer> memo) {
        if (n <= pre.limit)
            return pre.pref[(int) n];
        if (memo.containsKey(n))
            return memo.get(n);

        long res = (n % MOD) * ((n + 1) % MOD) % MOD;
        res = (res * INV2) % MOD;

        for (long l = 2; l <= n;) {
            long q = n / l;
            long r = n / q;
            long cnt = (r - l + 1) % MOD;
            long sub = (cnt * sumPhi(q, pre, memo)) % MOD;
            res = (res - sub + MOD) % MOD;
            l = r + 1;
        }

        int finalRes = (int) res;
        memo.put(n, finalRes);
        return finalRes;
    }

    static int G(long N, PhiPrefix pre, Map<Long, Integer> memo) {
        sumPhi(N, pre, memo);
        long ans = 0;
        for (long l = 1; l <= N;) {
            long q = N / l;
            long r = N / q;
            long sumD = triMod(l, r);
            long sphi = sumPhi(q, pre, memo);
            ans = (ans + sumD * sphi) % MOD;
            l = r + 1;
        }
        return (int) ans;
    }

    public static String solve() {
        int LIMIT = 5000000;
        PhiPrefix pre = buildPhiPrefix(LIMIT);
        Map<Long, Integer> memo = new HashMap<>();

        long ans = G(100000000000L, pre, memo);
        return Long.toString(ans);
    }

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