[HAOI2011]Problem b
给出n个询问，询问\(\sum_{i=a}^b\sum_{j=c}^d(gcd(i,j)==k)\),100%的数据满足：1≤n≤50000，1≤a≤b≤50000，1≤c≤d≤50000，1≤k≤50000。
解
法一:Mobius反演
设
\[b\leq d,f(k)=\sum_{i=a}^b\sum_{j=c}^d(gcd(i,j)==k)
\]
\[F(k)=\sum_{i=a}^b\sum_{j=c}^d(k|gcd(i,j))=
\]
\[([b/k]-[(a-1)/k])([d/k]-[(c-1)/k])
\]
由Mobius反演定理有
\[ans=f(k)=\sum_{k|x}F(x)\mu(x/k)=
\]
\[\sum_{k|x}([b/x]-[(a-1)/x])([d/x]-[(c-1)/x])\mu(x/k)=
\]
\[\sum_{x=1}^{[b/k]}([\frac{b}{xk}]-[\frac{a-1}{kx}])([\frac{d}{xk}]-[\frac{c-1}{xk}])\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[50001],pt,mb[50001];
il int min(int,int);
il void prepare(int);
int main(){
int a,b,c,d,k,i,j,ans(0),lsy;
prepare(50000),scanf("%d",&lsy);
while(lsy--){
scanf("%d%d%d%d%d",&a,&b,&c,&d,&k);
b/=k,d/=k,(--a)/=k,(--c)/=k,ans&=0;
if(b>d)swap(b,d),swap(a,c);
for(i=1;i<=b;i=j+1){
j=min(b/(b/i),d/(d/i));
if(a/i)j=min(j,a/(a/i));
if(c/i)j=min(j,c/(c/i));
ans+=(b/i-a/i)*(d/i-c/i)*(mb[j]-mb[i-1]);
}printf("%d\n",ans);
}
return 0;
}
il int min(int a,int b){
return a<b?a:b;
}
il void prepare(int n){
int i,j;check[1]|=mb[1]|=true;
for(i=2;i<=n;++i){
if(!check[i])prime[++pt]=i,mb[i]=-1;
for(j=1;j<=pt&&prime[j]<=n/i;++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];
}
法二：容斥原理
注意到下标的开始不为1，而我们能够做的是下标从1开始的问题，于是设\((a,b)\)表示a，b范围内，下界为1的满足题意的gcd限制的方案数，不难由容斥原理得
\[ans=(c,d)-(a-1,d)-(b,c-1)+(a-1,c-1)
\]
写出单个无下界限制的函数，套用公式即可。
参考代码：
#include <iostream>
#include <cstdio>
#define il inline
#define ri register
#define ll long long
#define swap(x,y) x^=y^=x^=y
using namespace std;
bool check[50001];
int prime[50001],pt,mb[50001];
il void read(int&);
il int min(int,int);
il ll mobius(int,int,int);
void prepare(int),pen(ll);
int main(){
int lsy,a,b,c,d,p;
read(lsy),prepare(50000);while(lsy--){
read(a),read(b),read(c),read(d),read(p);
pen(mobius(b,d,p)-mobius(a-1,d,p)-
mobius(b,c-1,p)+mobius(a-1,c-1,p)),putchar('\n');
}
return 0;
}
il ll mobius(int a,int b,int p){
int i,j;ll ans(0);a/=p,b/=p;
if(a>b)swap(a,b);
for(i=1;i<=a;i=j+1)
j=min(a/(a/i),b/(b/i)),
ans+=(ll)(a/i)*(b/i)*(mb[j]-mb[i-1]);
return ans;
}
void prepare(int n){
int i,j;check[1]|=mb[1]|=true;
for(i=2;i<=n;++i){
if(!check[i])prime[++pt]=i,mb[i]=-1;
for(j=1;i*prime[j]<=n&&j<=pt;++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;
}
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();
}
void pen(ll x){
if(x>9)pen(x/10);putchar(x%10+48);
}
