Description：

求 $\begin{align}\sum_{i=1}^n\sum_{j=1}^md(ij)\end{align}$

$1\le N,M\le 10^9$

Solution：

$\begin{align}\sum_{i=1}^n\sum_{j=1}^nd(ij)\end{align}$

$\begin{align}=\sum_{i=1}^n\sum_{j=1}^n\sum_{i'|i}\sum_{j'|j}[gcd(i',j')==1]\end{align}$

$\begin{align}=\sum_{i=1}^n\sum_{j=1}^n\sum_{i'|i}\sum_{j'|j}\sum_{k|i',k|j'}\mu(k)\end{align}$

$\begin{align}=\sum_{k=1}^n\mu(k)\sum_{i=1}^n\sum_{j=1}^n\sum_{k|i',i'|i}\sum_{k|j',j'|j}1\end{align}$

$\begin{align}=\sum_{k=1}^n\mu(k)\sum_{i=1}^{\lfloor \frac n k\rfloor}\sum_{j=1}^{\lfloor \frac n k\rfloor}\sum_{i'|i}\sum_{j'|j}1\end{align}$

$\begin{align}=\sum_{k=1}^n\mu(k)\sum_{i'}^{\lfloor \frac n k\rfloor}\sum_{i'|i}^{\lfloor \frac n k\rfloor}\sum_j^{\lfloor \frac n k\rfloor}\sum_{j'|j}^{\lfloor \frac n k\rfloor}1\end{align}$

$\begin{align}=\sum_{k=1}^n\mu(k)\sum_{i'}^{\lfloor \frac n k\rfloor}\sum_{i=1}^{\lfloor \frac n {ki'}\rfloor}\sum_j^{\lfloor \frac n k\rfloor}\sum_{j=1}^{\lfloor \frac n {kj'}\rfloor}1\end{align}$

$\begin{align}=\sum_{k=1}^n\mu(k)\times\sum_{i'}^{\lfloor \frac n k\rfloor}{\lfloor \frac n {ki'}\rfloor}\times\sum_j^{\lfloor \frac n k\rfloor}{\lfloor \frac n {kj'}\rfloor}\end{align}$

设 $\begin{align}g(n)=\sum_{i=1}^n\lfloor\frac n i\rfloor\end{align}$ ， $g$ 可以除法分块求

则原式 $=\begin{align}\sum_{k=1}^n\mu(k)g(\lfloor\frac n k\rfloor)^2\end{align}$ 也可以除法分块。

数据范围太大没办法线性筛，用杜教筛优化求前缀和的过程。

杜教筛 $O(n^\frac 2 3log_2n)$ ，两层分块求 $g$ 复杂度 $O(n^\frac 3 4)$ 。

Code：

#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<map>

#include<cstring>

using namespace std;

#define MOD 1000000007ll

int n;

#define MAXN 1000010

int pr[MAXN],tot = 0,mu[MAXN];

bool isprime[MAXN];

typedef long long ll;

map<int,ll> _mu,_g;

void sieve()

{

	for(int i = 2;i < MAXN;++i)isprime[i] = true;

	mu[1] = 1;

	for(int i = 2;i < MAXN;++i)

	{

		if(isprime[i])

		{

			pr[++tot] = i;

			mu[i] = -1;

		}

		for(int j = 1;j <= tot && i * pr[j] < MAXN;++j)

		{

			isprime[i * pr[j]] = false;

			if(i % pr[j] == 0)

			{

				mu[i * pr[j]] = 0;

				break;

			}

			else

			{

				mu[i * pr[j]] = -1 * mu[i];

			}

		}

	}

	for(int i = 1;i < MAXN;++i)

	{

		mu[i] = (mu[i] + mu[i - 1]) % MOD;

	}

	return;

}

ll sum(int n)

{

	if(n < MAXN)return mu[n];

	map<int,ll>::iterator it;

	if((it = _mu.find(n)) != _mu.end())return (it -> second);

	ll res = 1;

	for(int l = 2,r = 0;l <= n;l = r + 1)

	{

		r = n / (n / l);

		res = (res - (r - l + 1) * sum(n / l) % MOD + MOD) % MOD;

	}

	return _mu[n] = res;

}

ll g(int n)

{

	map<int,ll>::iterator it;

	if((it = _g.find(n)) != _g.end())return (it -> second);

	ll res = 0;

	for(int l = 1,r = 0;l <= n;l = r + 1)

	{

		r = n / (n / l);

		res = (res + (ll)(r - l + 1) * (ll)(n / l) % MOD) % MOD;

	}

	return _g[n] = (res * res % MOD);

}

int main()

{

	sieve();

	scanf("%d",&n);

	ll ans = 0;

	for(int l = 1,r = 0;l <= n;l = r + 1)

	{

		r = n / (n / l);

		ans = (ans + (sum(r) - sum(l - 1) + MOD) % MOD * g(n / l) % MOD) % MOD;

	}

	cout << ans << endl;

	return 0;

}