Description：

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

$1\le N \le 10^{10}$

Solution：

反演枚举倍数求和 $(bzoj2694\ Lcm)$ $O(nlog_2n)$ 显然是不能承受的。

发现式子是对称的，于是

原式 $\begin{align}=2\times\sum_{i=1}^n\sum_{j=1}^{i-1}lcm(i,j)-\frac {n\times(n+1)}{2}\end{align}$

$\begin{align}=2\times\sum_{i=1}^n\sum_{j=1}^{i-1}\frac{ij}{gcd(i,j)}-\frac {n\times(n+1)}{2}\end{align}$

$\begin{align}=2\times\sum_{d=1}^n\sum_{i=1}^n\sum_{j=1}^{i-1}\frac{ij}d[gcd(i,j)=d]-\frac {n\times(n+1)}{2}\end{align}$

$\begin{align}=2\times\sum_{d=1}^nd\sum_{i=1}^{\lfloor \frac n d\rfloor}i\sum_{j=1}^{i-1}j[gcd(i,j)=1]-\frac {n\times(n+1)}{2}\end{align}$

这个时候不能再反演了，发现第三个求和号表示的就是与 $i$ 互质的数的和。

原式 $\begin{align}=2\times\sum_{d=1}^nd\sum_{i=1}^{\lfloor\frac n d\rfloor}i\frac{i \varphi(i)+\varepsilon(i)} 2 -\frac {n\times(n+1)}{2}\end{align}$

$\begin{align}=\sum_{d=1}^nd(\sum_{i=1}^{\lfloor\frac n d\rfloor}i^2\times\varphi(i)+1)-\frac {n\times(n+1)}{2}\end{align}$

$\begin{align}=\sum_{d=1}^nd\sum_{i=1}^{\lfloor\frac n d\rfloor}i^2\times\varphi(i)\end{align}$

设 $\begin{align}g(n)=\sum_{i=1}^ni^2\times\varphi(i)\end{align}$

原式 $\begin{align}=\sum_{d=1}^nd\times g(\lfloor\frac n d\rfloor)\end{align}$ ，可以除法分块。

问题在于怎么求 $g(n)$ ，即怎么在低于线性的时间复杂度内求出 $i^2\times\varphi(i)$ 的前缀和。

不难想到杜教筛，发现 $\begin{align}(id^2\times\varphi)\times(id^2)=\sum_{d|n}d^2\times\varphi(d)\times(\frac n d)^2=n^2\times\sum_{d|n}\varphi(d)=id^3\end{align}$

于是就可以筛了。

需要用到一些公式：

$\begin{align}\sum_{i=1}^ni=\frac {n\times (n+1)} 2\end{align}$

$\begin{align}\sum_{i=1}^ni^2=\frac {n\times (n+1)\times(2n+1)} 2\end{align}$

$\begin{align}\sum_{i=1}^ni^3=\frac {n^4+2n^3+n^2}{4}\end{align}$

模数是质数需要求逆元。

#include<algorithm>

#include<iostream>

#include<cstdlib>

#include<cstdio>

#include<cmath>

#include<map>

#include<cstring>

using namespace std;

#define MAXN 4641600

#define MOD 1000000007ll

typedef long long ll;

ll pr[MAXN],tot = 0;

ll phi[MAXN],g[MAXN];

bool isprime[MAXN];

map<ll,ll> _g;

void sieve()

{

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

	phi[1] = 1;g[1] = 1;

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

	{

		if(isprime[i])

		{

			pr[++tot] = i;

			phi[i] = i - 1;

		}

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

		{

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

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

			{

				phi[i * pr[j]] = phi[i] * pr[j];

				break;

			}

			else

			{

				phi[i * pr[j]] = phi[i] * (pr[j] - 1);

			}

		}

		g[i] = (i % MOD) * (i % MOD) % MOD * (phi[i] % MOD) % MOD;

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

	}

	return;

}

ll n;

ll sum1(ll n)

{

	n %= MOD;

	return n * (n + 1) % MOD * 500000004 % MOD;

}

ll sum2(ll n)

{

	n %= MOD;

	return n * (n + 1) % MOD * (2 * n + 1) % MOD * 166666668 % MOD;

}

ll sum3(ll n)

{

	n %= MOD;

	return ((n * n % MOD * n % MOD * n % MOD + 2 * n % MOD * n % MOD * n % MOD) % MOD + n * n % MOD) % MOD * 250000002 % MOD;

}

ll calc(ll n)

{

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

	if(_g[n])return _g[n];

	ll res = sum3(n);

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

	{

		r = n / (n / l);

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

	}

	return _g[n] = res;

}

int main()

{

	sieve();

	scanf("%lld",&n);

	ll ans = 0;

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

	{

		r = n / (n / l);

		ans = (ans + (sum1(r) - sum1(l - 1) + MOD) % MOD * calc(n / l) % MOD) % MOD;

	}

	cout << ans << endl;

	return 0;

}