P2568 GCD（莫比乌斯反演）

题目传送门
本题题意转化成为：
$$\sum _{i=1}^n\sum _{j=1}^n[gcd(i,j)==p]$$
这道题用莫比乌斯式子推其实也挺简单的说…………欧拉函数的那个还理解了半天我真菜

$step:1$提取公因式
$$\sum _{p=1}^{p\in prime,p<=n}\sum _{i=1}^{[n/p]}\sum _{j=1}^{[n/p]}\epsilon (gcd(i,j))$$

$step:2$利用迪利克雷卷积进行反演
$$\sum _{p=1}^{p\in prime,p<=n}\sum _{i=1}^{[n/p]}\sum _{j=1}^{[n/p]}\sum _{k|i,k|j}\mu (k)$$
$step:3$交换枚举顺序
$$\sum _{p=1}^{p\in prime,p<=n}\sum _{k=1}^{[n/p]}\sum _{i=1}^{[n/(p\ast k)]}\sum _{j=1}^{[n/(p\ast k)]}\mu (k)$$
$$\sum _{p=1}^{p\in prime,p<=n}\sum _{k=1}^{[n/p]}\mu (k)\ast (n/(p\ast k){)}^2$$
然后利用我们的整除分块就可以愉快的实现我们的代码了

#include<iostream>
#include<vector> 
using namespace std;
typedef long long ll;
const ll MAXN=1e7+10;
//莫比乌斯函数
ll mu[MAXN];
bool isnp[MAXN];
vector<ll> primes;
void init(ll n)
{
    mu[1] = 1;
    for (ll i = 2; i <= n; i++)
    {
        if (!isnp[i])
            primes.push_back(i), mu[i] = -1; // 质数为-1
        for (ll p : primes)
        {
            if (p * i > n)
                break;
            isnp[p * i] = 1;
            if (i % p == 0)
            {
                mu[p * i] = 0; // 有平方因数为0
                break;
            }
            else
                mu[p * i] = mu[p] * mu[i]; // 互质，利用积性函数性质
        }
    }
    for(ll i=1;i<=n;i++){
    	mu[i]+=mu[i-1];
	}

}
int main(){
	ll n,ans=0;
	init(10000000);
	cin>>n;
	for(ll i=0;primes[i]<=n;i++){
		//cout<<primes[i]<<endl;
		ll u=n/primes[i];
		ll l=1,r;
		for(ll l=1,r;l<=u;l=r+1){
			r=u/(u/l);
			ans+=(mu[r]-mu[l-1])*(u/l)*(u/l);
		}
	}
	cout<<ans<<endl;
}