[bzoj] 1101 Zap || 莫比乌斯反演

原题

求x<=a，y<=b，并且gcd(x,y)=d的x，y对数

---

\(\sum^{n}_{x=1}\sum^{m}_{y=1}[gcd(x,y)==d]\)
\(=\sum^{\lfloor \frac{n}{d} \rfloor}_{x=1}\sum^{\lfloor \frac{m}{d} \rfloor}_{y=1}[gcd(x,y)==1]\)
根据莫比乌斯反演，可以把式子化为
\(=\sum^{\lfloor \frac{n}{d} \rfloor}_{x=1}\sum^{\lfloor \frac{m}{d} \rfloor}_{y=1}\sum_{d|gcd(x,y)}\mu(d)\)
因为$d|gcd(x,y) $ => \(d|x,d|y\)
所以式子变为 \(\sum^{n/d}_{i=1}\mu(d)*\lfloor \frac{\lfloor \frac{n}{d} \rfloor}{d} \rfloor*\lfloor \frac{\lfloor \frac{m}{d} \rfloor}{d} \rfloor\)

用n/d的取值分一个段，然后对mu求一个前缀和即可解决……

分段：
因为我们知道\(\lfloor \frac{n}{d} \rfloor\)的取值有\(\sqrt n\)种取值，所以对于\(\lfloor \frac{n}{d} \rfloor\)和\(\lfloor \frac{m}{d} \rfloor\)一共有\(\sqrt n + \sqrt m\)种情况
经典枚举方法：

for (int i=1,last=0;i<=a;i=last+1)
	{
	    last=min(a/(a/i),b/(b/i));
	    ans+=(ll)(a/i)*(b/i)*(sum[last]-sum[i-1]);
	}

对mu求前缀和：
利用欧拉筛法：

void init()
{
    miu[1]=sum[1]=1;
    for (int i=2;i<=N;i++)
    {
	if (!notprime[i]) miu[i]=-1,prime[++tot]=i;
	for (int j=1;j<=tot && i*prime[j]<=N;j++)
	{
	    notprime[i*prime[j]]=1;
	    if (i%prime[j]) miu[i*prime[j]]=-miu[i];
	    else
	    {
		miu[i*prime[j]]=0;
		break;
	    }
	}
	sum[i]=sum[i-1]+miu[i];
    }
}

AC代码：

#include<cstdio>
#include<algorithm>
#define N 50000
typedef long long ll;
using namespace std;
int t,a,b,d,miu[N+10],prime[N+10],tot,sum[N];
bool notprime[N+10];
ll ans;

int read()
{
    int ans=0,fu=1;
    char j=getchar();
    for (;j<'0' || j>'9';j=getchar()) if (j=='-') fu=-1;
    for (;j>='0' && j<='9';j=getchar()) ans*=10,ans+=j-'0';
    return ans*fu;
}

void init()
{
    miu[1]=sum[1]=1;
    for (int i=2;i<=N;i++)
    {
	if (!notprime[i]) miu[i]=-1,prime[++tot]=i;
	for (int j=1;j<=tot && i*prime[j]<=N;j++)
	{
	    notprime[i*prime[j]]=1;
	    if (i%prime[j]) miu[i*prime[j]]=-miu[i];
	    else
	    {
		miu[i*prime[j]]=0;
		break;
	    }
	}
	sum[i]=sum[i-1]+miu[i];
    }
}

int main()
{
    init();
    t=read();
    while (t--)
    {
	a=read();b=read();d=read();
	a/=d;
	b/=d;
	ans=0;
	if (a>b) swap(a,b);
	for (int i=1,last=0;i<=a;i=last+1)
	{
	    last=min(a/(a/i),b/(b/i));
	    ans+=(ll)(a/i)*(b/i)*(sum[last]-sum[i-1]);
	}
	printf("%lld\n",ans);
    }
    return 0;
}