Project Euler Problem 625

Solution

Answer: 551614306

Let

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

We use the standard identity

$$\sum_{i=1}^n \gcd(i,n)=\sum_{d\mid n} d,\varphi!\left(\frac nd\right),$$

where $\varphi$ is Euler's totient function.

Hence

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

Write $n=dk$. Then

$$G(N) =\sum_{k\le N}\varphi(k)\sum_{d\le N/k} d.$$

Using

$$\sum_{d\le m} d = \frac{m(m+1)}2,$$

we get

$$G(N)=\sum_{k=1}^N \varphi(k),\frac{t(t+1)}2, \qquad t=\left\lfloor \frac Nk \right\rfloor.$$

To evaluate this efficiently for $N=10^{11}$, we group equal values of $\left\lfloor N/k\right\rfloor$ and use the summatory totient function

$$\Phi(x)=\sum_{n\le x}\varphi(n).$$

The standard recursion is

$$\Phi(n)=\frac{n(n+1)}2-\sum_{l=2}^{n}(r-l+1)\Phi!\left(\left\lfloor\frac nl\right\rfloor\right),$$

with quotient grouping $r=\left\lfloor n/\lfloor n/l\rfloor\right\rfloor$.

The following Python program computes the answer:

from functools import lru_cache

MOD = 998244353
N = 10**11

# Precompute phi and prefix sums up to a cutoff
M = 3_000_000

phi = [0] * (M + 1)
primes = []
is_comp = bytearray(M + 1)

phi[1] = 1

# Linear sieve for Euler phi
for i in range(2, M + 1):
    if not is_comp[i]:
        primes.append(i)
        phi[i] = i - 1

    for p in primes:
        x = i * p
        if x > M:
            break

        is_comp[x] = 1

        if i % p == 0:
            phi[x] = phi[i] * p
            break
        else:
            phi[x] = phi[i] * (p - 1)

# Prefix sums: phi[n] = sum_{k<=n} phi(k)
for i in range(2, M + 1):
    phi[i] += phi[i - 1]

@lru_cache(maxsize=None)
def Phi(n):
    """Summatory totient function."""
    if n <= M:
        return phi[n]

    res = n * (n + 1) // 2

    l = 2
    while l <= n:
        v = n // l
        r = n // v
        res -= (r - l + 1) * Phi(v)
        l = r + 1

    return res

def G(N):
    ans = 0
    l = 1

    while l <= N:
        v = N // l
        r = N // v

        # sum phi(k) over [l, r]
        block_phi = Phi(r) - Phi(l - 1)

        tri = v * (v + 1) // 2

        ans = (ans + (block_phi % MOD) * (tri % MOD)) % MOD

        l = r + 1

    return ans

print(G(10))        # 122
print(G(10**11))

The program reproduces the given check:

$$G(10)=122.$$

Running it for $10^{11}$ gives

$$G(10^{11}) \bmod 998244353 = 551614306.$$

Answer: 551614306