Description
对于给出的n个询问,每次求有多少个数对(x,y),满足a≤x≤b,c≤y≤d,且gcd(x,y) = k,gcd(x,y)函数为x和y的最大公约数。
Input
第一行一个整数n,接下来n行每行五个整数,分别表示a、b、c、d、k
Output
共n行,每行一个整数表示满足要求的数对(x,y)的个数
Sample Input
2
2 5 1 5 1
1 5 1 5 2
Sample Output
14
3
HINT
100%的数据满足:1≤n≤50000,1≤a≤b≤50000,1≤c≤d≤50000,1≤k≤50000
算是写的第一道反演题。题目要求的是一段区间,通常的做法是求出1到右边界,然后再用容斥来求解区间。
我们考虑如何求1到右边界,设右边界为 n n 和,则要求的式子是 ∑ni=1∑mj=1[gcd(i,j)=k] ∑ i = 1 n ∑ j = 1 m [ g c d ( i , j ) = k ] ,我们来推下式子。
∑ni=1∑mj=1[gcd(i,j)=k]=∑nki=1∑mkj=1[gcd(i,j)=1]=∑nki=1∑mkj=1∑x|gcd(i,j)μ(x) ∑ i = 1 n ∑ j = 1 m [ g c d ( i , j ) = k ] = ∑ i = 1 n k ∑ j = 1 m k [ g c d ( i , j ) = 1 ] = ∑ i = 1 n k ∑ j = 1 m k ∑ x | g c d ( i , j ) μ ( x ) (根据性质1)
接下来要做最重要的一步转换,大家要好好理解: ∑nki=1∑mkj=1∑x|gcd(i,j)μ(x)=∑min(n,m)kx=1μ(x)[nkx][mkx] ∑ i = 1 n k ∑ j = 1 m k ∑ x | g c d ( i , j ) μ ( x ) = ∑ x = 1 m i n ( n , m ) k μ ( x ) [ n k x ] [ m k x ]
这步大家可能不太理解,我刚看的时候也是一脸懵逼,简单的说就是将 x x 提到外面,而右边表示的是 x x 的倍数的个数。
我们发现可以用数论分块,所以直接分块枚举x即可。
这道题还有另外一种推法,需要构造。
我们令 f(x)=∑ni=1∑mj=1[gcd(i,j)=x] f ( x ) = ∑ i = 1 n ∑ j = 1 m [ g c d ( i , j ) = x ] , F(x)=∑ni=1∑mj=1[x|gcd(i,j)] F ( x ) = ∑ i = 1 n ∑ j = 1 m [ x | g c d ( i , j ) ]
那么显然 F(x)=∑x|df(d) F ( x ) = ∑ x | d f ( d ) ,那么根据莫比乌斯反演 f(x)=∑x|dμ(dx)F(d) f ( x ) = ∑ x | d μ ( d x ) F ( d )
我们考虑 F(x) F ( x ) 是什么, F(x) F ( x ) 必须要求 i,j i , j 都是 x x 的倍数,所以
所以答案就是 f(k)=∑k|xμ(xk)[nx][mx]=∑min(n,m)kx=1μ(x)[nkx][mkx] f ( k ) = ∑ k | x μ ( x k ) [ n x ] [ m x ] = ∑ x = 1 m i n ( n , m ) k μ ( x ) [ n k x ] [ m k x ]
与上面的做法殊途同归。
#include<bits/stdc++.h>
using namespace std;
int read(){
char c;int x;while(c=getchar(),c<'0'||c>'9');x=c-'0';
while(c=getchar(),c>='0'&&c<='9') x=x*10+c-'0';return x;
}
void print(int x){
if(x/10) print(x/10);
putchar(x%10+'0');
}
int T,a,b,c,d,k,top,ans,mu[50005],sum[50005],pri[50005],vis[50005];
int calc(int n,int m){
if(n>m) swap(n,m);
int l=1,r=0,res=0;n=n/k;m=m/k;
while(l<=n){
r=min(n/(n/l),m/(m/l));r=min(r,n);
res+=(n/l)*(m/l)*(sum[r]-sum[l-1]);
l=r+1;
}
return res;
}
int main()
{
T=read();mu[1]=sum[1]=1;
for(int i=2;i<=50000;i++){
if(!vis[i]) pri[++top]=i,mu[i]=-1;
for(int j=1;j<=top&&pri[j]*i<=50000;j++){
vis[i*pri[j]]=1;
if(i%pri[j]==0) break;
mu[i*pri[j]]=-mu[i];
}
sum[i]=sum[i-1]+mu[i];
}
while(T--){
a=read();b=read();c=read();d=read();k=read();
ans=calc(b,d)-calc(a-1,d)-calc(b,c-1)+calc(a-1,c-1);
print(ans);puts("");
}
return 0;
}