51nod1238 最小公倍数之和 V3 莫比乌斯函数 杜教筛

题意：求\(\sum_{i = 1}^{n}\sum_{j = 1}^{n}lcm(i, j)\).
　　题解：虽然网上很多题解说用mu卡不过去，，，不过试了一下貌似时间还挺充足的，。。也许有时间用phi试试？
　　因为是用的莫比乌斯函数求的，所以推导比大部分题解多。。。而且我写式子一般都比较详细，所以可能看上去很多式子，实际上是因为每一步都写了，几乎没有跳过的。所以应该都可以看懂的。
　　末尾的\(e\)函数是指的\(e[1] = 1\),\(e[x] = 0(x != 1)\)这样一个函数
　　$$\sum_{i = 1}^{n}\sum_{j = 1}^{n}lcm(i, j)$$
　　$$\sum_{i = 1}^{n} \sum_{i = 1}^{n} \frac{ij}{gcd(i, j)}$$
　　枚举\(gcd\)
　　$$\sum_{d = 1}^{n} \sum_{i = 1}^{\lfloor {\frac{n}{d}} \rfloor} \sum_{j = 1}^{\lfloor {\frac{n}{d}} \rfloor}[gcd(i, j) == d] \frac{ij}{d}$$
　　因为\((\frac{ijd^2}{d} = ijd)\)，所以：
　　$$\sum_{d = 1}^{n} \sum_{i = 1}^{\lfloor {\frac{n}{d}} \rfloor} \sum_{j = 1}^{\lfloor {\frac{n}{d}} \rfloor}[gcd(i, j) == d] ijd$$
　　$$\sum_{d = 1}^{n}d\sum_{i = 1}^{\lfloor {\frac{n}{d}} \rfloor} \sum_{j = 1}^{\lfloor {\frac{n}{d}} \rfloor}ij[gcd(i, j) == 1]$$
　　$$\sum_{d = 1}^{n} d \sum_{i = 1}^{\lfloor {\frac{n}{d}} \rfloor} \sum_{j = 1}^{\lfloor {\frac{n}{d}} \rfloor} ij \sum_{k | gcd(i, j)} \mu(k)$$
　　枚举k,再枚举k的倍数。
　　$$\sum_{d = 1}^{n}d\sum_{k = 1}^{\lfloor {\frac{n}{d}} \rfloor}\mu(k) \sum_{i = 1}^{\lfloor {\frac{n}{dk}} \rfloor}ik \sum_{j = 1}^{\lfloor {\frac{n}{dk}} \rfloor}jk$$
　　设\(S(n) = \sum_{i = 1}^{n}i\)
　　$$\sum_{d = 1}^{n}d \sum_{k = 1}^{\lfloor {\frac{n}{d}} \rfloor} \mu(k) k ^ 2 S(\frac{n}{dk})$$
　　枚举\(T = dk\)
　　$$\sum_{T = 1}^{n} S(\frac{n}{T})^ 2 \sum_{k | T} \mu(k) k ^ 2 \frac{T}{k}$$
　　$$\sum_{T = 1}^{n} S(\frac{n}{T})^ 2 \sum_{k | T} \mu(k) kT$$
　　$$\sum_{T = 1}^{n} S(\frac{n}{k})^ 2 \cdot T \sum_{k | T} \mu(k)k$$
　　设\(f(T) = T\sum_{k | T} \mu(k) k\)，卷上\(id^2\)，因为\(S(\frac{n}{k})\)可以数论分块，所以我们只需要快速求出区间\([l, r]\)内的\(f\)之和即可,显然求出\(f\)的前缀和即可解决问题

　　$$(f * id^2)(n) = \sum_{i |n}f(i) \frac{n ^ 2}{i^2}=\sum_{i | n}i \sum_{k | i} \mu(k) k \frac{n ^ 2}{i ^ 2}$$
　　$$\sum_{i | n}\sum_{k | i} \mu(k)k\frac{n ^ 2}{i} = n \sum_{i | n}\sum_{k | i} \mu(k) k \frac{n}{i}$$
　　设$$h(i) = \sum_{k | i} \mu(k)k$$，则原式：
　　$$n \sum_{i | n} h(i) \frac{n}{i} = n (h * id)(n)$$
　　$$(f * id ^ 2)(n) = n (h * id)(n)$$
　　$$h(n) = \sum_{k | n}\mu(k)k = (\mu \cdot id) * 1$$
　　$$f * id ^ 2 = n [(\mu \cdot id) * 1 * id] = n[(\mu \cdot id) * id * 1]$$
　　其中$$(\mu \cdot id) * id = \sum_{i | n} \mu(i) i \frac{n}{i} = n \sum_{i | n}\mu(i) = e$$
　　所以
　　$$n[(\mu \cdot id) * id * 1] = n[e * 1] = n$$
　　带入杜教筛的式子：
　　$$g(1)S(n) = \sum_{i = 1}^{n} (f * g)(i) - \sum_{i = 2}^{n}g(i)S(\frac{n}{i})$$
　　$$= \sum_{i = 1}^{n}i - \sum_{i = 2}^{n}i ^ 2 S(\frac{n}{i})$$
　　然后直接上杜教就可以了.
　　其实还有一个问题。。。一开始预处理的前缀和怎么求?
　　要知道前缀和，首先要求出\(f\).
　　因为\(f(T) = T\sum_{k | T} \mu(k) k\)，所以如果我们可以快速求出\(\sum_{k | T}\mu(k)k\),然后就只需要再\(O(n)\)的乘上\(T\)就可以了.
　　我们先预处理出\(\mu(k)\)，然后对于每一个\(k\),枚举它的倍数，统计贡献。那么复杂度为 \(\frac{n}{1} + \frac{n}{2} + ... + \frac{n}{n} = nlogn\)(此处的\(n\)为原题面的\(\frac{2}{3}\)次方，即要预处理的\(f\)个数)

#include<bits/stdc++.h>
using namespace std;
#define R register int
#define LL long long
#define RL register LL
#define AC 3000
#define ac 5000000
#define p 1000000007LL
//#define h(x) ((x <= block) ? sum[x] : S[n / x])

LL n, ans, block; 
LL mu[ac], S[AC], sum[ac], inv[AC];
int pri[ac], tot;
bool z[ac], vis[AC];

inline LL h(LL x)
{
	return ((x <= block) ? sum[x] : S[n / x]);
}

inline void up(LL & a, LL b)
{
	a += b;
	if(a >= p) a -= p;
	if(a <= -p) a += p;
}

LL count(LL l, LL r){
	return (r - l + 1) % p * ((r + l) % p) % p * inv[2] % p;
}

void pre()
{
	scanf("%lld", &n), block = pow(n, 0.66666);
	mu[1] = 1;
	for(R i = 2; i <= block; i ++)
	{
		if(!z[i]) pri[++ tot] = i, mu[i] = -1;
		for(R j = 1; j <= tot; j ++)
		{
			int now = pri[j];
			if(i * now > block) break;
			z[i * now] = true;
			if(!(i % now)) break;
			mu[i * now] = - mu[i];
		}
	}
	inv[1] = 1;
	for(R i = 2; i <= 10; i ++) inv[i] = (p - p / i) * inv[p % i] % p;
	for(R i = 1; i <= block; i ++)//枚举mu(i)
		for(R j = 1; j; j ++)//枚举i的倍数
		{
			if(j * i > block) break;
			up(sum[i * j], mu[i] * i % p);
		}
	for(R i = 1; i <= block; i ++) sum[i] = sum[i] * i % p; 
	for(R i = 1; i <= block; i ++) up(sum[i], sum[i - 1]);//算出f数组后还要统计前缀和
}

LL get(LL x)
{
	x %= p;
	return x * (x + 1) % p * (2 * x + 1) % p * inv[6] % p;
}

void cal(LL x)
{
	if(x <= block || vis[n / x]) return ;
	LL rnt = count(1, x);
	for(RL i = 2, lim, now; i <= x; i = lim + 1)
	{
		lim = x / (x / i), now = x / i, cal(now);
		up(rnt, - ((get(lim) - get(i - 1)) % p * h(now) % p));
	}
	S[n / x] = rnt, vis[n / x] = true;
}

void work()
{
	for(RL i = 1, lim, now, x; i <= n; i = lim + 1)
	{
		lim = n / (n / i), now = n / i, x = count(1, now);
		up(ans, (x % p * x % p * ((h(lim) - h(i - 1)) % p) % p));
	}
	printf("%lld\n", (ans + p) % p);
}

int main()
{
	//freopen("in.in", "r", stdin);
	pre();
	cal(n);
	work();
//	fclose(stdin);
	return 0;
}