给出n个询问，询问(sum_{i=a}^bsum_{j=c}^d(gcd(i,j)==k)),100%的数据满足：1≤n≤50000，1≤a≤b≤50000，1≤c≤d≤50000，1≤k≤50000。
解
法一:Mobius反演
设
[bleq d,f(k)=sum_{i=a}^bsum_{j=c}^d(gcd(i,j)==k)
]
[F(k)=sum_{i=a}^bsum_{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
",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('
');
}
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);
}