I was reading a lot of stuff about Dirichlet Convolution and how to solve sums over arithmetic functions. There are a lot of resources and I don't consider that it's an easy topic. Even when I don't use Dirichlet here, It's something that I'll try to practice more because it's really useful when you're working with Mobius function. Maybe I need to review my past codes to see if I can improve them.

Statement

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

Find $G(10^{11})$. Give your answer modulo $998244353$.

Solution

Naive solution

Let's try with the simplest and naive solution, just doing a nested for:

long long G = 0;
for(long long j=1; i<=N; i++){
    for(long long i=1; i<=j; i++) G += __gcd(i, j);
}

Of course, the time complexity of this solution is: $\mathcal{O}(N^{2})$ and not feasible to compute since the value for $N$ is $10^{11}$.

Linear solution

From my previous solution in a past blog, using the Mobius function it is possible to reduce the solution to $\mathcal{O}(n\log(n))$ and by applying the Mertens function optimization, it is possible to reduce it to $\mathcal{O}(n)$.

Even when this time complexity is acceptable (it could take a few hours), the memory utilization is close to the order of terabytes, so it is not efficient.

long long G = 0; 
vector <int> values;
vector <pair <int,int> > range;
fill(n, values, range);

for(int i=0; i<values.size(); i++){
    int v = values[i];
    long long aux = 0;
    vector <int> n_v; 
    vector <pair <int,int> > n_r;
    fill(v, n_v,n_r);
    for(int j=0; j<n_v.size(); j++){
        int k = n_v[j];
        long long s = 1LL*k*k-(1LL*k*(k+1))/2;
        aux += s*(s_mobius[n_r[j].first] - s_mobius[n_r[j].second-1]);
    }
    long long s = (1LL*range[i].first*(range[i].first+1))/2;
    s -= (1LL*range[i].second*(range[i].second-1))/2;
    G += aux*s;
}

Sub-linear solution

Since even the linear solution is not acceptable for this problem. We need to find a sub-linear approach, in this case, I will try to formulate a $\mathcal{O}(n^{2/3})$ solution, which can be computed in few seconds.

First, I'll calculate the sum of $\gcd(i, n)$:

$$g(n) = \sum_{i=1}^{n} \gcd(i, n)\\ = \sum_{k\vert n} k \phi\Big( \frac{n}{k} \Big)\\ = (n * \phi)(n)$$

The last line is the dirichlet notation, it could be useful.

Now, we can use $g(n)$ to calculate $G(N)$:

$$G(N) = \sum_{i=1}^{N} g(i)\\ = \sum_{i=1}^{N} \sum_{k \vert i} \phi(k) \frac{i}{k}\\ = \sum_{k=1}^{N} \sum_{p=1}^{\frac{N}{k}} \phi(k) p\\ = \sum_{k=1}^{N} \sum_{p=1}^{\frac{N}{k}} \phi(k) p\\ = \sum_{k=1}^{N} \phi(k) \frac{\lfloor \frac{N}{k}\rfloor (\lfloor \frac{N}{k}\rfloor + 1)}{2}$$

As you can see, we want to calculate sums of $\phi(k)$, so let's to reduce it:

$$S(n) = \sum_{i=1}^{n} \phi(i) \\ S(n) = \Phi(n) \\ = \sum_{b=1}^{n} \sum_{a=1}^{b} \big\vert \gcd(a, b) = 1 \big\vert \\ = \frac{n(n+1)}{2} - \sum_{m=2}^{n} \sum_{y=1}^{\lfloor n/m \rfloor} \sum_{x=1}^{\lfloor b/m \rfloor} \big\vert \gcd(x, y) = 1 \big\vert \\ = \frac{n(n+1)}{2} - \sum_{m=2}^{n} S(\lfloor n/m \rfloor) \\ = \frac{n(n+1)}{2} - \sum_{m=2}^{\lfloor \sqrt{n} \rfloor} S(\lfloor n/m \rfloor) - \sum_{m=1}^{\lfloor \frac{n}{\lfloor \sqrt{n} \rfloor} \rfloor - 1} \Big(\Big\lfloor\frac{n}{m} \Big\rfloor - \Big\lfloor\frac{n}{m+1} \Big\rfloor \Big) S(m)$$

This is a recursive definition. It's not hard to see that the time complexity for this is $\mathcal{O}(n^{2/3})$.

Similar to the explained in a previous post, the list of possible values of $\lfloor N/k \rfloor$ has a size of $2\sqrt{N}$. I will use the hyperbola method to calculate it:

$$G(N) = \sum_{k=1}^{N} \phi(k) \frac{\lfloor \frac{N}{k}\rfloor (\lfloor \frac{N}{k}\rfloor + 1)}{2}\\ = \sum_{k=1}^{N} \phi(k) \sum_{ab=k} a \phi(b) \\ = \sum_{a=1}^{\lfloor \sqrt{N} \rfloor} \sum_{b=1}^{\lfloor N/a\rfloor} a\phi(b) + \sum_{b=1}^{\lfloor \sqrt{N}\rfloor} \sum_{a=1}^{\lfloor N/b\rfloor} a\phi(b) - \sum_{a=1}^{\lfloor\sqrt{N}\rfloor} \sum_{b=1}^{\lfloor\sqrt{N}\rfloor} a\phi(b)\\ = \sum_{a=1}^{\lfloor\sqrt{N}\rfloor} a \Phi(\lfloor N/a \rfloor) + \sum_{b=1}^{\lfloor \sqrt{N} \rfloor} \phi(b) \frac{\lfloor N/b\rfloor \times (\lfloor N/b\rfloor + 1)}{2} - \frac{\lfloor \sqrt{N}\rfloor (\lfloor \sqrt{N} \rfloor + 1)}{2}\Phi(\lfloor \sqrt{N} \rfloor)\\$$

To optimize the complexity, suppose that we will do a sieve until $N/k$, so the overall complexity will be:

$$\mathcal{O}(G(N)) = \mathcal{O}\Big(\sum_{a=1}^{k} \mathcal{O}\big(\Phi(\lfloor N/a\rfloor)\big)\Big) + \mathcal{O}(\sqrt{N}) + \mathcal{O}(N/k \log(N/k))\\ = \mathcal{O}\Big(\sum_{a=1}^{k} 1\Big) + \mathcal{O}(N^{2/3}) + \mathcal{O}(\sqrt{N}) + \mathcal{O}(N/k \log(N/k))\\ = \mathcal{O}(k) + \mathcal{O}(N^{2/3}) + \mathcal{O}(\sqrt{N}) + \mathcal{O}(N/k \log(N/k))\\ \text{Optimal value is when:}\\ k \approx N^{69/125}\\ \mathcal{O}(G(N)) = \mathcal{O}(k + N^{2/3} + \sqrt{N} + N/k\log(N/k))\\ = \mathcal{O}(N^{2/3})$$

Code

Here is the code, it's not so fast due to mod operations.

#include <bits/stdc++.h>
using namespace std;
using ll = long long;
#define N 10000000
#define MOD 998244353

int phi[N];
int c_Phi[N];
unordered_map <ll, int> m_Phi;
void cphi(){
    phi[1] = 1;
    for(int i=2; i<N; i++){
        if(!phi[i]){
            phi[i] = i-1;
            for(ll j=2*i; j<N; j+=i){
                if(phi[j] == 0) phi[j] = j;
                phi[j] = phi[j]/i*(i-1);
            }
        }
    }
    c_Phi[1] = phi[1];
    for(int i=2; i<N; i++) c_Phi[i] = (c_Phi[i-1] + phi[i])%MOD;
}

int T(ll n){
    n %= MOD;
    ll ans = n*(n+1);
    ans %= MOD;
    ans *= 499122177;
    return ans%MOD;
}

ll Phi(ll n){
    if(n < N) return c_Phi[n];
    if(m_Phi.find(n) != m_Phi.end()) return m_Phi[n];
    ll ans = T(n);
    ll sqrt_n = floor(sqrt(n));
    for(int m=1; m < n/sqrt_n; m++){
        ans -= 1LL*(n/m - n/(m+1))*c_Phi[m];
        ans %= MOD;
        }
    for(int m=2; m<=sqrt_n; m++){
        ans -= Phi(n/m);
        ans %= MOD;
    }
    ans += MOD; ans %= MOD;
    m_Phi[n] = ans;
    return ans;
}

ll G(ll n){
    ll ans = 0;
    int sqrt_n = floor(sqrt(n));
    ans -= 1LL*T(sqrt_n)*Phi(sqrt_n);
    ans %= MOD;
    for(int i=1; i<=sqrt_n; i++){
        ans += Phi(n/i)*i;
        ans %= MOD;
    }

    for(int i=1; i<=sqrt_n; i++){
        ans += 1LL*phi[i]*T(n/i);
        ans %= MOD;
    }
    ans += MOD; ans %= MOD;
    return ans;
}

int main(){
    cphi();
    ll n; cin >> n;
    cout << G(n) << '\n';
    return 0;
}

Resource:

I want to list them because I would like to go deeper into them and should review them further at my leisure.

• griff's blog
• Dirichlet Convolution
• Beni's blog
• codeforces blog