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}} floor} sum_{j = 1}^{lfloor {frac{n}{d}} floor}[gcd(i, j) == d] frac{ij}{d}$$
　　因为((frac{ijd^2}{d} = ijd))，所以：
　　$$sum_{d = 1}^{n} sum_{i = 1}^{lfloor {frac{n}{d}} floor} sum_{j = 1}^{lfloor {frac{n}{d}} floor}[gcd(i, j) == d] ijd$$
　　$$sum_{d = 1}^{n}dsum_{i = 1}^{lfloor {frac{n}{d}} floor} sum_{j = 1}^{lfloor {frac{n}{d}} floor}ij[gcd(i, j) == 1]$$
　　$$sum_{d = 1}^{n} d sum_{i = 1}^{lfloor {frac{n}{d}} floor} sum_{j = 1}^{lfloor {frac{n}{d}} floor} ij sum_{k | gcd(i, j)} mu(k)$$
　　枚举k,再枚举k的倍数。
　　$$sum_{d = 1}^{n}dsum_{k = 1}^{lfloor {frac{n}{d}} floor}mu(k) sum_{i = 1}^{lfloor {frac{n}{dk}} floor}ik sum_{j = 1}^{lfloor {frac{n}{dk}} floor}jk$$
　　设(S(n) = sum_{i = 1}^{n}i)
　　$$sum_{d = 1}^{n}d sum_{k = 1}^{lfloor {frac{n}{d}} floor} 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) = Tsum_{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)kfrac{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) = Tsum_{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
", (ans + p) % p);
}

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