莫比烏斯反演入門 HDOJ 1695:GCD 、BZOJ 2301: [HAOI2011]Problem b
阿新 • • 發佈:2019-01-06
下面我所說的都基於上面這篇部落格的內容。
莫比烏斯反演有兩種形式(mu表示莫比烏斯函式):
HDOJ1695 GCD
求1<=x<=n,1<=y<=m中gcd(x,y)==k的(x,y)組數,注意(a,b)和(b,a)視為同一情況。
相當於計算1<=x<=n/k,1<=y<=m/k的數中有多少對互質(即gcd(x,y)==1)。
進一步轉化出所求即為
f(1)=
令n/=k,m/=k
算出1<=x<=n,1<=y<=m中的互質情況數,減去1<=x<=min(n,m)、1<=y<=min(n,m)中互質情況的一半就是答案
為什麼是min(n,m)呢?因為(a,b)、(b,a)這種重複情況,a、b肯定屬於同一範圍,而且這一範圍是較小的那一個範圍。
假設n<m,x在1~n中取,y在1~m中取,如果y進一步取在n-1~m中,由於和x範圍不交叉,所以肯定沒有重複情況。
所以有重複情況的(x,y),一定在公共(交叉)區間內。
因此上界取min(n,m),省去了一定計算量。
程式碼體現在:
ans1+=(LL)mu[i]*(b/i)*(b/i);
完整程式碼
#include <iostream> #include <cstdio> #include <algorithm> #include <cmath> #include <cstring> #define mst(a,b) memset(a,b,sizeof(a)) typedef long long LL; using namespace std; const int N = 100005; bool check[N+10]; int prime[N+10],mu[N+10]; //莫比烏斯函式模板 void Mobius() { memset(check,false,sizeof(check)); mu[1] = 1; int tot = 0; for(int i=2;i<=N;i++){ if(!check[i]){ prime[tot ++] = i; mu[i] = -1; } for(int j=0;j<tot;j++){ if(i * prime[j] > N) 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() { int a,b,c,d,k,T,cas=1; Mobius(); cin>>T; while(T--){ scanf("%d%d%d%d%d",&a,&b,&c,&d,&k); if(k==0){ printf("Case %d: 0\n",cas++); continue; } b/=k; d/=k; if(b>d) swap(b,d); LL ans=0,ans1=0; for(int i=1;i<=b;i++){ ans+=(LL)mu[i]*(b/i)*(d/i); ans1+=(LL)mu[i]*(b/i)*(b/i); } ans-=ans1/2; printf("Case %d: %lld\n",cas++,ans); } return 0; }
BZOJ 2301: [HAOI2011]Problem b
每次求有多少個數對(x,y),滿足a≤x≤b,c≤y≤d,且gcd(x,y) = k,gcd(x,y)函式為x和y的最大公約數。
#include <iostream> #include <cstdio> #include <algorithm> #include <cmath> #include <cstring> #define mst(a,b) memset(a,b,sizeof(a)) typedef long long LL; using namespace std; const int N = 100000; bool check[N+5]; int prime[N+5],mu[N+5],sum[N+5]; //莫比烏斯函式模板 void Mobius() { memset(check,false,sizeof(check)); mu[1] = 1; int tot = 0; for(int i=2;i<=N;i++){ if(!check[i]){ prime[tot ++] = i; mu[i] = -1; } for(int j=0;j<tot;j++){ if(i * prime[j] > N) break; check[i * prime[j]] = true; if(i % prime[j] == 0){ mu[i * prime[j]] = 0; break; } else mu[i * prime[j]] = -mu[i]; } } } LL cal(int n,int m,int k){ LL ans=0; n/=k; m/=k; if(n>m) swap(n,m); for(int i=1,last;i<=n;i=last+1){ //i是不變塊的左邊界,last是不變塊的右邊界 last = min(n/(n/i),m/(m/i)); //根據左邊界算出右邊界 ans+=(LL)(sum[last]-sum[i-1])*(n/i)*(m/i); } return ans; } int main() { int a,b,c,d,k,T; Mobius(); sum[0]=0; for(int i=1;i<=N;i++) sum[i]=sum[i-1]+mu[i]; cin>>T; while(T--){ scanf("%d%d%d%d%d",&a,&b,&c,&d,&k); if(k==0){ printf("0\n"); continue; } LL ans1=cal(b,d,k); LL ans2=cal(a-1,d,k); LL ans3=cal(b,c-1,k); LL ans4=cal(a-1,c-1,k); LL ans=ans1-ans2-ans3+ans4; //容斥原理 printf("%lld\n",ans); } return 0; }