2301: [HAOI2011]Problem b/1101: [POI2007]Zap 莫比乌斯反演

今天决定入坑莫比乌斯反演qwq。
《多年的心头大恨终于解决了系列》
更多关于莫比乌斯反演的问题我打算一会儿写到[总结]数学专题里去
我们用f(k)表示1≤x≤n,1≤y≤m且gcd(x,y)=k的数对(x,y)的个数，并假设n≤m。
我们发现直接计算f(k)是比较困难的。
引入一个g(k)表示1≤x≤n,1≤y≤m且k∣gcd(x,y)的数对(x,y)的个数，则有g(k)=∑⌊nk⌋d=1f(k×d)，g(k)是显然的，g(k)=⌊nk⌋×⌊mk⌋。
根据莫比乌斯反演，
g(k)=∑⌊nk⌋d=1f(k×d)
=>f(k)=∑⌊nk⌋d=1μ(d)×g(k×d)=∑⌊nk⌋d=1μ(d)×⌊nk×d⌋×⌊mk×d⌋。
然后我们分块求解f(k)，然后用前缀和求出答案。

#include<iostream>#include<cstdio>#define MAXN 50005using namespace std;int sum[MAXN],mobius[MAXN],prime[MAXN];bool flag[MAXN];inline int read(){    int a=0,f=1; char c=getchar();    while (c<'0'||c>'9') {if (c=='-') f=-1; c=getchar();}    while (c>='0'&&c<='9') {a=a*10+c-'0'; c=getchar();}    return a*f;}inline void prepare(){    mobius[1]=1;    for (int i=2;i<MAXN;i++)    {        if (!flag[i]) prime[++prime[0]]=i,mobius[i]=-1;        for (int j=1;j<=prime[0]&&i*prime[j]<MAXN;j++)        {            flag[i*prime[j]]=1;            if (i%prime[j]==0) {mobius[i*prime[j]]=0; break;}            mobius[i*prime[j]]=-mobius[i];        }    }    for (int i=1;i<MAXN;i++) sum[i]=sum[i-1]+mobius[i];}inline int calc(int n,int m){    if (n>m) swap(n,m);    int ans=0,pos;    for (int i=1;i<=n;i=pos+1)    {        pos=min(n/(n/i),m/(m/i));        ans+=(sum[pos]-sum[i-1])*(n/i)*(m/i);    }    return ans;}int main(){    prepare();    int testcase=read();    while (testcase--)    {        int a=read()-1,b=read(),c=read()-1,d=read(),k=read();        a/=k; b/=k; c/=k; d/=k;        int ans=calc(a,c)+calc(b,d)-calc(a,d)-calc(b,c);        printf("%d\n",ans);    }    return 0;}