【BZOJ2820】YY的GCD-莫比乌斯反演+数论分块+分类讨论

测试地址：YY的GCD
题目大意：求∑p=prime∑ni=1∑mj=1[gcd(i,j)=p]$\sum _{p=prime}\sum _{i=1}^n\sum _{j=1}^m[gcd(i,j)=p]$，T(≤104)$T(\le {10}^4)$组询问，n,m≤107$n,m\le {10}^7$。
做法：本题需要用到莫比乌斯反演+数论分块+分类讨论。
首先按照套路推式子（不妨设n≤m$n\le m$）：
ans=∑p=prime∑ni=1∑mj=1[gcd(i,j)=p]$ans=\sum _{p=prime}\sum _{i=1}^n\sum _{j=1}^m[gcd(i,j)=p]$
=∑p=prime∑n/pi=1∑m/pj=1[gcd(i,j)=1]$=\sum _{p=prime}\sum _{i=1}^{n/p}\sum _{j=1}^{m/p}[gcd(i,j)=1]$
=∑p=prime∑n/pi=1∑m/pj=1∑d|i,d|jμ(d)$=\sum _{p=prime}\sum _{i=1}^{n/p}\sum _{j=1}^{m/p}\sum _{d|i,d|j}\mu (d)$
=∑p=prime∑n/pd=1μ(d)⌊npd⌋⌊mpd⌋$=\sum _{p=prime}\sum _{d=1}^{n/p}\mu (d)\lfloor \frac{n}{pd}\rfloor \lfloor \frac{m}{pd}\rfloor$
令k=pd$k=pd$，有：
ans=∑p=prime∑p|kμ(kp)⌊nk⌋⌊mk⌋$ans=\sum _{p=prime}\sum _{p|k}\mu (\frac{k}{p})\lfloor \frac{n}{k}\rfloor \lfloor \frac{m}{k}\rfloor$
互换p,k$p,k$的位置，有：
ans=∑nk=1⌊nk⌋⌊mk⌋∑p=primeμ(kp)$ans=\sum _{k=1}^n\lfloor \frac{n}{k}\rfloor \lfloor \frac{m}{k}\rfloor \sum _{p=prime}\mu (\frac{k}{p})$
这是一个数论分块的形式，如果我们能预处理出g(k)=∑p=primeμ(kp)$g(k)=\sum _{p=prime}\mu (\frac{k}{p})$，我们就可以以O(Tn−−√)$O(T\sqrt{n})$的复杂度完成本题了。
直接枚举是O(nlogn)$O(n\log n)$的，无法承受，而g$g$显然并不是一个积性函数，这是不是就意味着我们无法用线性筛求了呢？不是的，只是我们需要进一步探求g$g$的性质。
令k=p0x$k=p_{0}x$，其中p0$p_0$为k$k$的最小质因子，我们尝试用和x$x$有关的式子O(1)$O(1)$算出g(k)$g(k)$。
情况一：当p0|x$p_0|x$时，考虑在g(k)$g(k)$式子中枚举的p$p$：
当p0=p$p_0=p$时，显然μ(kp)=μ(x)$\mu (\frac{k}{p})=\mu (x)$；
当p0≠p$p_0\ne p$时，kp$\frac{k}{p}$包含因数p20$p_0^2$，所以μ(kp)=0$\mu (\frac{k}{p})=0$。
综上所述，g(k)=μ(x)$g(k)=\mu (x)$。
情况二：当p0∤x$p_0\nmid x$时，考虑在g(k)$g(k)$式子中枚举的p$p$：
当p0=p$p_0=p$时，显然μ(kp)=μ(x)$\mu (\frac{k}{p})=\mu (x)$；
当p0≠p$p_0\ne p$时，因为p0$p_0$是质数，且p0∤x$p_0\nmid x$，所以p0$p_0$和x$x$互质，那么p0$p_0$和xp$\frac{x}{p}$显然也互质。又因为μ$\mu$是积性函数，所以μ(kp)=μ(p0)⋅μ(xp)$\mu (\frac{k}{p})=\mu (p_0)\cdot \mu (\frac{x}{p})$。将μ(p0)$\mu (p_0)$提出和式后，发现剩下的就是g(x)$g(x)$的表达式。
综上所述，g(k)=μ(x)−g(x)$g(k)=\mu (x)-g(x)$。
于是我们就可以用线性筛O(n)$O(n)$求出g$g$了，这样我们就解决了这一题。
以下是本人代码：

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
int T,prime[10000010];
ll mxn=0,n[10010],m[10010];
ll mu[10000010],g[10000010];
bool vis[10000010]={0};

void calc(int n)
{
    mu[1]=1,g[1]=0;
    prime[0]=0;
    for(int i=2;i<=n;i++)
    {
        if (!vis[i])
        {
            prime[++prime[0]]=i;
            mu[i]=-1;
            g[i]=1;
        }
        for(int j=1;j<=prime[0]&&i*prime[j]<=n;j++)
        {
            vis[i*prime[j]]=1;
            if (i%prime[j]==0)
            {
                mu[i*prime[j]]=0;
                g[i*prime[j]]=mu[i];
                break;
            }
            mu[i*prime[j]]=-mu[i];
            g[i*prime[j]]=mu[i]-g[i];
        }
    }
    for(int i=1;i<=n;i++)
        g[i]+=g[i-1];
}

int main()
{
    scanf("%d",&T);
    for(int i=1;i<=T;i++)
    {
        scanf("%lld%lld",&n[i],&m[i]); 
        if (n[i]>m[i]) swap(n[i],m[i]);
        mxn=max(mxn,n[i]);
    }

    calc(mxn);

    for(int i=1;i<=T;i++)
    {
        ll ans=0;
        for(ll j=n[i];j;j=max(n[i]/(n[i]/j+1),m[i]/(m[i]/j+1)))
        {
            ll l=max(n[i]/(n[i]/j+1),m[i]/(m[i]/j+1)),r=j;
            ans+=(n[i]/j)*(m[i]/j)*(g[r]-g[l]);
        }
        printf("%lld\n",ans);
    }

    return 0;
}