题目描述
题解
人生第一次正儿八经的莫比乌斯反演?
题目要求
∑i=1b∑j=1d[(i,j)=k]
变形
∑i=1b∑j=1d[k|i][k|j][(ik,jk)=1]
令i=ik,j=jk得
∑i=1bk∑j=1dk[(i,j)=1]
然后我们把式子化成这样
∑i=1bk∑j=1dk∑t|(i,j)μ(t)=∑i=1bk[t|i]∑j=1dk[t|j]μ(t)=∑t=1max(bk,dk)⌊bkt⌋⌊dkt⌋μ(t)
就可以用分块解决问题了?
不不不,我们写这道题是为了练反演
我们设f(n)为最大公约数为n的数对个数,F(n)为公约数为n的数对个数,即最大公约数为n的倍数的数的个数,则f(1)就是我们要求的答案。
存在关系
F(n)=∑n|df(d)
利用莫比乌斯反演公式可以得到
f(n)=∑n|dμ(dn)F(d)
那么接下来的问题就是求解F(n)
由F(n)的定义可以得到公式
F(n)=∑i=1b∑j=1d[n|(i,j)]
其实也就是[1,b][1,d]范围内可以同时整除n的数对个数
设[1,b]内可以整除n的数个数为p,[1,d]内可以整除n的数的个数为q,那么由乘法原理可得F(n)=p*q
那么利用结论
p=⌊bn⌋,q=⌊dn⌋
可以得到
F(n)=⌊bn⌋∗⌊dn⌋
那么我们现在就可以O(1)得到F(n)的函数值
那么我们回到原来的公式
f(n)=∑n|dμ(dn)F(d)
我们要求解f(1),那么公式变成了
f(1)=∑μ(d)F(d)
那么我们只要枚举d就可以了
看上去d好像是没有上限的,但是仔细观察公式 F(n)=⌊bn⌋∗⌊dn⌋ 可以知道,当枚举的d>min(b,d)(注意前面的d是枚举的莫比乌斯反演公式中的d,后面的b,d是题目中给的)时,F(n)的函数值为0,就不会再计算了。所以枚举d的上限为min(b,d)
这里的 μ 是可以线性筛出来的,F(n)是可以O(1)算出来的,算法的思路就很清晰了。
不过这道题还牵扯到一个“去重”问题,假设b < <script type="math/tex" id="MathJax-Element-140"><</script>d,显然根据以上的计算步骤,只有[1,b]范围内的数对被重复计算了。所以我们需要计算出[1,b]范围内不重复的满足条件的数对的个数,再用刚才算出的f(1)减去它。而这就相当于d和b相等,再用同样的方法计算一遍f(1)即可。因为两个上界的范围相同,都是b,所以所有数对都被计算了两次,可以直接除以2来获得不重复的答案。时间复杂度是O(n)的。
感谢ATP神犇的帮助!顺便%一发!
代码
#include<iostream>
#include<cstring>
#include<cstdio>
using namespace std;
#define LL long long
const int N=1e5;
int T,a,b,c,d,k;
LL ans,re;
int mu[N+5],p[N+5],prime[N+5];
inline void get_mu(){
mu[1]=1;
for (int i=2;i<=N;++i){
if (!p[i]){
prime[++prime[0]]=i;
mu[i]=-1;
}
for (int j=1;j<=prime[0]&&i*prime[j]<=N;++j){
p[i*prime[j]]=1;
if (i%prime[j]==0){
mu[i*prime[j]]=0;
break;
}
else mu[i*prime[j]]=-mu[i];
}
}
}
inline LL F(int n){
return (LL)(b/n)*(d/n);
}
int main(){
get_mu();
scanf("%d",&T);
for (int t=1;t<=T;++t){
scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
if (k==0){printf("Case %d: %d\n",t,0);continue;}
b/=k; d/=k; if (b>d) swap(b,d);
ans=re=0;
for (int i=1;i<=b;++i) ans+=mu[i]*F(i);
d=b;
for (int i=1;i<=b;++i) re+=mu[i]*F(i);
ans-=re/2;
printf("Case %d: %I64d\n",t,ans);
}
}