题意:这题求[1,b],[1,d]gcd为k的对数。
思路:转化之后就是[1,b/k],[1,d/k]之间互质的数的个数。
设f(k)为gcd(x,y)=k的数对(x,y)的对数,我们要求的是f(1)
设g(k)为gcd(x,y)为k的倍数的数对(x,y)的对数,可以想到g(k)=floor(b/k)*floor(d/k),
由莫比乌斯反演得:
令lim=min(b/k,d/k)
f(1)=mu[1]*g(1) + mu[2]*g[2] + ... + mu[lim]*g(lim)
因为(n1,n2)和(n2,n1)算为同一种情况,所以最后结果还要减掉重复的情况。
ps:这道题还可以用容斥做。
代码1:朴素求法
#include<iostream> #include<stdio.h> #include<string.h> using namespace std; const int MAXN=1000000; bool check[MAXN+10]; int prime[MAXN+10]; int mu[MAXN+10]; void Mobius(){ memset(check,false,sizeof(check)); mu[1]=1; int tot=0; for(int i=2;i<=MAXN;i++){ if(!check[i]){ prime[tot++]=i; mu[i]=-1; } for(int j=0;j<tot;j++){ if(i*prime[j]>MAXN)break; check[i*prime[j]]=true; if(i%prime[j]==0){ mu[i*prime[j]]=0; break; } else{ mu[i*prime[j]]=-mu[i]; } } } } int main(){ Mobius(); int t,i; int a,b,c,d,k; scanf("%d",&t); for(i=1;i<=t;i++){ scanf("%d%d%d%d%d",&a,&b,&c,&d,&k); if(k==0){ printf("Case %d: 0 ",i); continue; } b/=k; d/=k; if(b>d)swap(b,d); long long ans1=0; for(int j=1;j<=b;j++) ans1+=(long long)mu[j]*(b/j)*(d/j); long long ans2=0; for(int j=1;j<=b;j++) ans2+=(long long)mu[j]*(b/j)*(b/j); ans1-=ans2/2; printf("Case %d: %I64d ",i,ans1); } return 0; }
代码2:分块优化
#include<iostream> #include<stdio.h> #include<string.h> using namespace std; const int MAXN=100000; bool check[MAXN+10]; int prime[MAXN+10]; int mu[MAXN+10]; void Mobius(){ memset(check,false,sizeof(check)); mu[1]=1; int tot=0; for(int i=2;i<=MAXN;i++){ if(!check[i]){ prime[tot++]=i; mu[i]=-1; } for(int j=0;j<tot;j++){ if(i*prime[j]>MAXN)break; check[i*prime[j]]=true; if(i%prime[j]==0){ mu[i*prime[j]]=0; break; } else{ mu[i*prime[j]]=-mu[i]; } } } } int sum[MAXN+10]; //找[1,n],[1,m]内互质的数的对数 long long solve(int n,int m){ long long ans=0; for(int i=1,la=0;i<=n;i=la+1){ la=min(n/(n/i),m/(m/i)); ans+=(long long)(sum[la]-sum[i-1])*(n/i)*(m/i); } return ans; } int main(){ Mobius(); sum[0]=0; for(int i=1;i<=MAXN;i++) sum[i]=sum[i-1]+mu[i]; int t,i; int a,b,c,d,k; scanf("%d",&t); for(i=1;i<=t;i++){ scanf("%d%d%d%d%d",&a,&b,&c,&d,&k); if(k==0){ printf("Case %d: 0 ",i); continue; } b/=k; d/=k; if(b>d)swap(b,d); long long ans1=solve(b,d); long long ans2=solve(b,b); ans1-=ans2/2; printf("Case %d: %I64d ",i,ans1); } return 0; }