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