求n组(sum_{i=1}^asum_{j=1}^b(gcd(i,j)==d)), 1≤n≤50 000,1≤d≤a,b≤50 000。
解
不难看出是约数组合计数问题，而解决该问题常用思路有容斥原理和Mobius反演。
法一：容斥原理
尽可能特殊化问题，令(a/=d,b/=d)，原问题等价于在这个范围里寻找互质的数的对数，而容斥原理关键在于构造不合理情况，即有别的约数，于是设(D(a,b,i))表示a，b范围里有公约数i的个数，不难得知(D(a,b,i)=[a/i][b/i]),于是由容斥原理有
[ans=sum_{i=1}^{min(a,b)}mu(i)D(a,b,i)=sum_{i=1}^{min(a,b)}mu(i)[a/i][b/i]
]
预处理(mu)函数前缀和，进行整出分块即可。
法二：Mobius反演
设
[f(d)=sum_{i=1}^asum_{j=1}^b(gcd(i,j)==d)
]
[F(d)=sum_{i=1}^asum_{j=1}^b(d|gcd(i,j))=[a/d][b/d]
]
由Mobius反演定理我们有
[f(d)=sum_{d|x}F(x)mu(x/d)=sum_{d|x}[a/x][b/x]mu(x/d)
]
于是
[ans=f(d)=sum_{d|x}[a/x][b/x]mu(x/d)=
]
[sum_{x=1}^{min([a/d],[b/d)]}[a/x][b/x]mu(x)
]
对式子进行整除分块即可。
参考代码
#include <iostream>
#include <cstdio>
#define il inline
#define ri register
#define swap(x,y) x^=y^=x^=y
using namespace std;
bool check[50001];
int prime[5137],pt,mb[50001];
il void read(int&);
il int min(int,int);
void pen(int),prepare(int);
int main(){
int lsy,a,b,n,i,j,ans;read(lsy);
prepare(50000);while(lsy--){
read(a),read(b),read(n),ans&=0;
a/=n,b/=n;if(a>b)swap(a,b);
for(i=1;i<=a;i=j+1)
j=min(a/(a/i),b/(b/i)),
ans+=(a/i)*(b/i)*(mb[j]-mb[i-1]);
pen(ans),putchar('
');
}
return 0;
}
void prepare(int n){
int i,j;check[1]|=mb[1]|=true;
for(i=2;i<=n;++i){
if(!check[i])mb[i]=-1,prime[++pt]=i;
for(j=1;j<=pt&&i*prime[j]<=n;++j){
check[i*prime[j]]|=true;
if(!(i%prime[j]))break;
mb[i*prime[j]]=-mb[i];
}
}for(i=1;i<=n;++i)mb[i]+=mb[i-1];
}
il int min(int a,int b){
return a<b?a:b;
}
void pen(int x){
if(x>9)pen(x/10);putchar(x%10+48);
}
il void read(int &x){
x&=0;ri char c;while(c=getchar(),c<'0'||c>'9');
while(c>='0'&&c<='9')x=(x<<1)+(x<<3)+(c^48),c=getchar();
}