【51nod】Gcd and Phi（反演）

基准时间限制：2 秒 空间限制：131072 KB 分值: 320
$F(n)=∑_{i=1}^n∑_{j=1}^nϕ(ϕ_i,ϕ_j)$
其中 ϕ 表示欧拉函数。欧拉函数ϕn 是不超过n的数中与n互质的数的数目。
ϕ(ϕi,ϕj) 表示i,j欧拉函数值的最大公约数的欧拉函数值.

给出n，求F(n)的值。
Input
一行，包含正整数T，表示数据组数。T<=5
接下来T行每行一个正整数n。n<=2*10^6
Output
共T行，每行包含一个答案。
Input示例
1
1
Output示例
1
思路：发现对于枚举的i是没有直接关系的，而是一个映射即$\phi(i)$,所以应该直接去枚举这个，所以我们设$g(x)为有多少\phi(i)$等于$x$,这个可以线性统计一下。
则转化为：

$$ans=\sum_{i=1}^n\sum_{j=1}^ng(i)g(j)\phi(gcd(i,j))$$
接下来显然我们肯定要去枚举gcd的，所以：

$$=\sum_{d=1}^n\phi(d)\sum_{d|i}^n\sum_{d|j}^ng(i)g(j)[gcd(i,j)=d]\\ =\sum_{d=1}^n\phi(d)\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}g(id)g(jd)[gcd(i,j)=1]$$
这里就有可以引入莫比乌斯函数了：
$$=\sum_{d=1}^n\phi(d)\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}g(id)g(jd)\sum_{k|gcd(i,j)}\mu(k)\\ =\sum_{d=1}^n\phi(d)\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{n}{d}}g(id)g(jd)\sum_{k|i,k|j}\mu(k)\\ =\sum_{d=1}^n\phi(d)\sum_{k=1}^{\frac{n}{d}}\mu(k)\sum_{i=1}^{\frac{n}{dk}}\sum_{j=1}^{\frac{n}{dk}}g(idk)g(jdk)$$
然后就可以搞了。

#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#define ull  unsigned long long
#define ll long long
#define ul unsigned int
#define maxn 2000005
#define mod 1000000007
using namespace std;
bool isP[maxn];
int prime[maxn];
int mu[maxn];
int phi[maxn];
int cnt;
void init()
{
    phi[1]=mu[1]=1;
    for(int i=2;i<maxn;i++)
    {
        if(!isP[i])
        {
            prime[cnt++]=i;
            phi[i]=i-1;
            mu[i]=-1;
        }
        for(int j=0;j<cnt&&(ll)i*prime[j]<maxn;j++)
        {
            int now=i*prime[j];
            isP[now]=true;
            if(i%prime[j])
            {
                phi[now]=phi[i]*(prime[j]-1);
                mu[now]=-mu[i];
            }
            else
            {
                phi[now]=phi[i]*prime[j];
                mu[now]=0;
                break;
            }
        }
    }
}
int g[maxn];
ll G[maxn];
int n;
void _init()
{
    memset(g,0,sizeof(g));
    memset(G,0,sizeof(G));
    for(int i=1;i<=n;++i)
        ++g[phi[i]];
    for(int i=1;i<=n;++i)
        for(int j=i;j<=n;j+=i)
            G[i]+=g[j];
    for(int i=1;i<=n;i++)
        G[i]*=G[i];
}
int main()
{
    init();
    int t;
    cin>>t;
    while(t--)
    {
        cin>>n;
        _init();
        ll ans=0;
        for(int d=1;d<=n;d++)
        {
            int limit=n/d;
            for(int k=1;k<=limit;k++)
                if(mu[k])
                    ans+=mu[k]*G[k*d]*phi[d];
        }
        cout<<ans<<endl;
    }
    return 0;
}