P3455 [POI2007]ZAP-Queries
斯德哥尔摩
2018-05-13 18:27:02
[P3455 [POI2007]ZAP-Queries](https://www.luogu.org/problemnew/show/P3455)
题目要求：$\sum_{i=1}^n\sum_{i=1}^m[gcd(i,j)==d]$
先保证：$n<=m$
设：$$f(d)=\sum_{i=1}^n\sum_{i=1}^m[gcd(i,j)==d]$$
$$F(d)=\sum_{i=1}^n\sum_{i=1}^m[d|gcd(i,j)]$$
一个数对要满足它们的$gcd$是$d$的倍数，则$i$和$j$中都必须包含$d$这个因子，所以有：
$$F(d)=\lfloor\frac{n}{d}\rfloor\lfloor\frac{m}{d}\rfloor$$
而$$F(d)=\sum_{id=1}^nf(i\times d)=\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}f(id)$$
然后反演一下：
$$f(d)=\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\mu(i)F(i\times d)$$
把$F(d)$带入：
$$f(d)=\sum_{i=1}^{\lfloor\frac{n}{d}\rfloor}\mu(i)\lfloor\frac{n}{id}\rfloor\lfloor\frac{m}{id}\rfloor$$
那个$id$很烦人，设$D=id$，则：
$$f(d)=\sum_{D=1}^n\mu(\frac{D}{d})\lfloor\frac{n}{D}\rfloor\lfloor\frac{m}{D}\rfloor$$
这个式子显然可以$O(n)$求解，由于是多组询问，我们用一下数论分块就可以做到$O(\sqrt n)$求解了。
附代码：
```cpp
#include<iostream>
#include<algorithm>
#include<cstdio>
#define MAXN 50010
using namespace std;
int n,m,d,k=0,prime[MAXN],mu[MAXN],sum[MAXN];
bool np[MAXN];
inline int read(){
int date=0,w=1;char c=0;
while(c<'0'||c>'9'){if(c=='-')w=-1;c=getchar();}
while(c>='0'&&c<='9'){date=date*10+c-'0';c=getchar();}
return date*w;
}
void make(){//莫比乌斯专属线性筛
int m=MAXN-10;
mu[1]=1;
for(int i=2;i<=m;i++){
if(!np[i]){
mu[i]=-1;
prime[++k]=i;
}
for(int j=1;j<=k&&prime[j]*i<=m;j++){
np[prime[j]*i]=true;
if(i%prime[j]==0)break;
else mu[prime[j]*i]=-mu[i];
}
}
for(int i=1;i<=m;i++)sum[i]=sum[i-1]+mu[i];
}
void work(){//求解
long long ans=0;
n=read();m=read();d=read();
if(n>m)swap(n,m);
n/=d;m/=d;
for(int i=1,last=1;i<=n;i=last+1){
last=min(n/(n/i),m/(m/i));
ans+=(long long)(sum[last]-sum[i-1])*(n/i)*(m/i);
}
printf("%lld\n",ans);
}
int main(){
int t=read();
make();
while(t--)work();//多组询问
return 0;
}
```