【BZOJ 2301】Problem B 莫比乌斯反演

题意

$n$个询问。
求$\sum_{i=a}^b\sum_{j=c}^d [\gcd(i,j)=k]$。
$1\leq n\leq 50000，≤a≤b≤50000，1≤c≤d≤50000，1≤k≤50000$

分析

用莫比乌斯反演求解这道题。
网上常见的分析有两种，一种是直接推一遍莫比乌斯反演，另一种是直接使用莫比乌斯反演。。

反演：
$F(n)=\sum_{n|d}f(d)\Rightarrow f(n)=\sum_{n|d}\mu({d\over n})F(d)$。

首先，根据容斥原理，把原问题转成4个子问题。
现在的问题相当于求：
$\sum_{i=1}^n \sum_{j=1}^m [\gcd(i,j)=k]$。

发现直接求解不容易，考虑使用反演。
设：$f(k)=\sum_{i=1}^n\sum_{j=1}^m [gcd(i,j)=k]$，即为所求的答案。
$F(k)=\sum_{i=1}^n \sum_{j=1}^m [k|\gcd(i,j)]=\lfloor {n\over k}\rfloor \lfloor {m\over k}\rfloor$。

则有$F(k)=\sum_{k|d}f(d)$。
$\therefore f(k)=\sum_{k|d}\mu({d\over k})F(d)=\sum_{k|d}\mu({d\over k})\lfloor {n\over d}\rfloor \lfloor {m\over d}\rfloor$

下一步在popoqqq的PPT里面没有讲清楚啊，然后就直接给代码了。不过这也是很容易想到的。
我们考虑枚举${d\over k}$的值$i$，则$d=ki$
$\therefore f(n)=\sum_{i=1}^n\mu(i)\lfloor {n\over {ki}}\rfloor \lfloor {m\over {ki}}\rfloor$。
到这里两种方法求出的结果已经相同了。

接下来的方法就不赘述了。
这里只是证明一下$\lfloor {n\over d}\rfloor$只有$O(\sqrt n)$个取值。
证明：①当$d\leq \sqrt n$时，因为只有$\sqrt n$个数，所以只有$\sqrt n$个取值。
②当$d>\sqrt n$时，${n\over d}<\sqrt n$，又因为是整数，所以最多只有$\sqrt n$个取值。
综上，有$O(\sqrt n)$个取值。

代码

#include <cstdio>
#include <cstring>
#include <cctype>
#include <algorithm>
using namespace std;

const int N=50001;

int vis[N],mu[N],pri[N],tot,suf[N];
int cas;
int a,b,c,d,k;
int s1,s2,s3,s4,ans;

inline int read(void)
{
    int x=0; char c=getchar();
    for (;!isdigit(c);c=getchar());
    for (;isdigit(c);c=getchar()) x=x*10+c-'0';
    return x;
}

int query(int n,int m)
{
    int a=n/k,b=m/k; int l,r; int calc=0;
    if (a>b) swap(a,b);
    for (l=1;l<=a;l=r+1)
    {
        r=min(a/(a/l),b/(b/l));
        calc=calc+(a/l)*(b/l)*(suf[r]-suf[l-1]);
    }
    return calc;
}

int main(void)
{
//  freopen("a.in","r",stdin);
//  freopen("a.out","w",stdout);

    vis[1]=mu[1]=1;
    for (int i=2;i<N;i++)
    {
        if (!vis[i])
        {
            mu[i]=-1;
            pri[++tot]=i;
        }
        for (int j=1;j<=tot;j++)
        {
            if (i*pri[j]>=N) break;
            vis[i*pri[j]]=1;
            if (i%pri[j]!=0)
                mu[i*pri[j]]=-mu[i];
            else
            {
                mu[i*pri[j]]=0;
                break;
            }
        }
    }
    for (int i=1;i<N;i++)
        suf[i]=suf[i-1]+mu[i];

    cas=read();
    for (int cc=1;cc<=cas;cc++)
    {
        a=read(),b=read(),c=read(),d=read(),k=read();
        s1=query(b,d);
        s2=query(a-1,d);
        s3=query(b,c-1);
        s4=query(a-1,c-1);
        ans=s1+s4-s2-s3;
        printf("%d\n",ans);
    }

    return 0;
}