51nod1237 最大公约数之和V3 杜教筛

首先肯定是求
$$\sum _{d=1}^{n}d\sum _{i=1}^n\sum _{j=1}^n[gcd(i,j)==d]$$
直接莫反一下
$$\sum _{d=1}^{n}d\sum _{i=1}^{[\frac{n}{d}]}\sum _{j=1}^{[\frac{n}{d}]}\sum _{T|gcd(i,j)}\mu [T]$$
把$T$提到外面，有
$$\sum _{d=1}^n\sum _{T=1}^{[\frac{n}{d}]}d\mu [T][\frac{n}{dT}{]}^2$$
我们可以直接对左边数论分块，那就变成求
$$\sum _{i=1}^n\mu [i][\frac{n}{i}{]}^2$$
这东西当然也可以直接数论分块套杜教筛就好啦。
然后我们算算复杂度，分块套分块是$O(n^{\frac{3}{4}})$的，杜教筛是$O(n^{2/3})$的。
我们对于第二个分块可以预处理前$O(n^{2/3})$个，这样总复杂度就是$O(n^{2/3})$了。
预处理的话考虑差分，设$S(n)={\sum }_{i=1}^n\mu [i][\frac{n}{i}{]}^2$
然后$$S(n)-S(n-1)=\sum _{d|n}\mu [d]\ast [(\frac{n}{d}{)}^2-(\frac{n}{d}-1{)}^2]$$
$$=2\ast \phi [n]-[n==1]$$
当然我们还可以对于式子继续化简。
化简的话考虑把$dT$换成$T$，那么有
$$\sum _{T=1}^n[\frac{n}{T}{]}^2\sum _{d|T}d\mu [\frac{T}{d}]$$
后面的实际上就是$id\ast \mu =\phi$
那么前面的数论分块，后面的杜教筛也行。

不化简考虑预处理的：

#include <iostream>
#include <cstdio>
#include <cmath>
#include <cstring>
#include <cstdlib>
#include <algorithm>
#include <unordered_map>
#include <ctime>

using namespace std;

using ll = long long;

int const mod = 1e9 + 7;
int const M = 1e7;
int const N = M + 5;

int f[N], pri[M / 3], tot, g[N];
bool isp[N];

void SAI() {
	f[1] = g[1] = 1;
	isp[1] = 1;
	for (int i = 2; i <= M; ++i) {
		if (isp[i] == 0) {
			pri[++tot] = i;
			f[i] = -1;
			g[i] = i - 1;
		}
		for (int j = 1; j <= tot && i * pri[j] <= M; ++j) {
			isp[i * pri[j]] = 1;
			if (i % pri[j]) {
				f[i * pri[j]] = -f[i];
				g[i * pri[j]] = g[i] * g[pri[j]];
			} else {
				f[i * pri[j]] = 0;
				g[i * pri[j]] = g[i] * pri[j];
				break;
			}
		}
	}
	for (int i = 1; i <= M; ++i) {
		f[i] = (1ll * f[i - 1] + 1ll * f[i] + mod * 2) % mod;
		g[i] = (1ll * g[i - 1] + 2ll * g[i]) % mod;
	}
	for (int i = 1; i <= M; ++i)
		g[i] = (g[i] + mod - 1) % mod;
}

unordered_map<ll, int> mp;

ll sum_(ll x) {
	if (x <= M)
		return f[x];
	if (mp.count(x))
		return mp[x];
	ll ret = 1;
	for (ll l = 2, r; l <= x; l = r + 1) {
		r = x / (x / l);
		ret += mod - sum_(x / l) * ((r - l + 1) % mod) % mod;
		if (ret >= mod)
			ret -= mod;
	}
	ret = (ret + mod) % mod;
	mp[x] = ret;
	return ret;
}

ll calc(ll n) {
	if (n <= M)
		return g[n];
	ll ret = 0;
	for (ll l = 1, r; l <= n; l = r + 1) {
		r = n / (n / l);
		ll nd = (n / l) % mod;
		ret += nd * nd % mod * (sum_(r) - sum_(l - 1)) % mod;
		if (ret < 0)
			ret += mod;
		if (ret >= mod)
			ret -= mod;
	}
	return (ret + mod) % mod;
}

int main() {
	ll n, ans = 0;
	cin >> n;
	SAI();
	for (ll l = 1, r; l <= n; l = r + 1) {
		r = n / (n / l);
		ll t = calc(n / l);
		ll t1 = l + r, t2 = r - l + 1;
		if (t1 & 1)
			t2 >>= 1;
		else
			t1 >>= 1;
		if (t1 >= mod)
			t1 %= mod;
		if (t2 >= mod)
			t2 %= mod;
		ans += t1 * t2 % mod * t % mod;
		if (ans >= mod)
			ans -= mod;
	}
	cout << ans << '\n';
}

化简后直接分块杜教筛的:：

#include <iostream>
#include <cstdio>
#include <cmath>
#include <cstring>
#include <cstdlib>
#include <algorithm>
#include <unordered_map>

using namespace std;

using ll = long long;

ll const mod = 1e9 + 7;
int const M = 1e7;
int const N = M + 5;

ll KSM(ll a, ll k) {
	ll ret = 1;
	for (; k; k >>= 1, a = a * a % mod)
		if (k & 1)
			ret = ret * a % mod;
	return ret;
}

int f[N], pri[M / 3], tot;
bool isp[N];
ll inv2;

void SAI() {
	f[1] = 1;
	isp[1] = 1;
	for (int i = 2; i <= M; ++i) {
		if (isp[i] == 0) {
			pri[++tot] = i;
			f[i] = i - 1;
		}
		for (int j = 1; j <= tot && 1ll * i * pri[j] <= M; ++j) {
			isp[i * pri[j]] = 1;
			if (i % pri[j]) 
				f[i * pri[j]] = f[i] * f[pri[j]];
			else {
				f[i * pri[j]] = f[i] * pri[j];
				break;
			}
		}
	}
	for (int i = 1; i <= M; ++i) 
		f[i] = (1ll * f[i - 1] + 1ll * f[i] + mod * 2) % mod;
}

unordered_map<ll, ll> mp;

ll sum_(ll x) {
	if (x <= M)
		return f[x];
	if (mp.count(x))
		return mp[x];
	ll ret = ((1 + x) % mod) * (x % mod) % mod * inv2 % mod;
	for (ll l = 2, r; l <= x; l = r + 1) {
		r = x / (x / l);
		ret = ((mod - sum_(x / l) * ((r - l + 1) % mod) % mod) + ret) % mod;
	}
	mp[x] = ret;
	return ret;
}

int main() {
	SAI();
	inv2 = KSM(2, mod - 2);
	ll n, ans = 0;
	cin >> n;
	for (ll l = 1, r; l <= n; l = r + 1) {
		r = n / (n / l);
		ll nd = (n / l) % mod;
		ans += nd * nd % mod * (sum_(r) - sum_(l - 1)) % mod;
		ans %= mod;
	}
	ans = (ans + mod) % mod;
	cout << ans << '\n';
}