算法笔记--莫比乌斯反演
这个B站上的视频讲的不错：http://www.bilibili.com/video/av14325327/
线性筛求莫比乌斯反演函数代码：
void Mobius() { mem(vis,false); mu[1]=1; cnt=0; for(int i=2;i<N;i++) { if(!vis[i]) { prime[cnt++]=i; mu[i]=-1; } for(int j=0;j<cnt&&i*prime[j]<N;j++) { vis[i*prime[j]]=true; if(i%prime[j])mu[i*prime[j]]=-mu[i]; else { mu[i*prime[j]]=0; break; } } } }
例题1：HDU 1695 GCD
代码：
#include<bits/stdc++.h> using namespace std; #define ll long long #define pb push_back #define mem(a,b) memset(a,b,sizeof(a)) const int N=1e5+5; int prime[N]; bool not_prime[N]={false}; int mu[N]; void Mobius() { mu[1]=1; int k=0; for(int i=2;i<N;i++) { if(!not_prime[i]) { prime[k++]=i; mu[i]=-1; } for(int j=0;i*prime[j]<N;j++) { not_prime[i*prime[j]]=true; if(i%prime[j]==0) { mu[i*prime[j]]=0; break; } else mu[i*prime[j]]=-mu[i]; } } } int main() { ios::sync_with_stdio(false); cin.tie(0); int T; int a,b,c,d,k; Mobius(); cin>>T; for(int Case=1;Case<=T;Case++) { cin>>a>>b>>c>>d>>k; if(k==0){ cout<<"Case "<<Case<<": 0"<<endl; continue; } ll ans1=0,ans2=0; b/=k; d/=k; for(int i=1;i<=min(b,d);i++)ans1+=(ll)mu[i]*(b/i)*(d/i); for(int i=1;i<=min(b,d);i++)ans2+=(ll)mu[i]*(min(b,d)/i)*(min(b,d)/i); ans1-=ans2/2; cout<<"Case "<<Case<<": "<<ans1<<endl; } return 0; }
代码：
#include<bits/stdc++.h> using namespace std; #define ll long long #define pb push_back #define mem(a,b) memset(a,b,sizeof(a)) const int N=1e5+5; const int INF=0x3f3f3f3f; const int MOD=1e9+7; int a[N]; int prime[N]; bool not_prime[N]={false}; int mu[N]; int cnt[2*N]; void Mobius() { mu[1]=1; int k=0; for(int i=2;i<N;i++) { if(!not_prime[i]) { prime[k++]=i; mu[i]=-1; } for(int j=0;i*prime[j]<N;j++) { not_prime[i*prime[j]]=true; if(i%prime[j]==0) { mu[i*prime[j]]=0; break; } else mu[i*prime[j]]=-mu[i]; } } } ll q_pow(ll n,ll k) { ll ans=1; while(k) { if(k&1)ans=(ans*n)%MOD; n=(n*n)%MOD; k>>=1; } return ans; } int main() { /*ios::sync_with_stdio(false); cin.tie(0);*/ int T; int n; Mobius(); scanf("%d",&T); for(int Case=1;Case<=T;Case++) { int mn=INF; int mx=0; mem(cnt,0); scanf("%d",&n); for(int i=1;i<=n;i++)scanf("%d",&a[i]),mn=min(mn,a[i]),mx=max(mx,a[i]); for(int i=1;i<=n;i++)cnt[a[i]]++; for(int i=1;i<2*N;i++)cnt[i]+=cnt[i-1]; ll ans=0; for(int i=2;i<=mn;i++) { ll t=1; for(int j=1;j*i<=mx;j++) { t=(t*q_pow(j,cnt[(j+1)*i-1]-cnt[j*i-1]))%MOD; } ans=(ans-mu[i]*t%MOD+MOD)%MOD; } printf("Case #%d: %lld\n",Case,ans); } return 0; }
例题3：Codeforces 547C - Mike and Foam
提示：C(n,2)=n*(n-1)/2=1+2+...+n-1
代码：
#include<bits/stdc++.h> using namespace std; #define ll long long #define pb push_back #define mem(a,b) memset(a,b,sizeof(a)) const int N=5e5+5; const int INF=0x3f3f3f3f; const int MOD=1e9+7; int a[N]; int prime[N]; bool not_prime[N]={false}; bool vis[N]={false}; int mu[N]; int cnt[N]; void Mobius() { mu[1]=1; int k=0; for(int i=2;i<N;i++) { if(!not_prime[i]) { prime[k++]=i; mu[i]=-1; } for(int j=0;i*prime[j]<N;j++) { not_prime[i*prime[j]]=true; if(i%prime[j]==0) { mu[i*prime[j]]=0; break; } else mu[i*prime[j]]=-mu[i]; } } } int main() { ios::sync_with_stdio(false); cin.tie(0); int n,q,t; Mobius(); cin>>n>>q; for(int i=1;i<=n;i++)cin>>a[i]; ll ans=0; while(q--) { cin>>t; if(!vis[t]) { vis[t]=true; for(int i=1;i*i<=a[t];i++) { if(a[t]%i==0) { ans+=(ll)mu[i]*cnt[i]; cnt[i]++; if(i*i!=a[t])ans+=mu[a[t]/i]*cnt[a[t]/i],cnt[a[t]/i]++; } } } else { vis[t]=false; for(int i=1;i*i<=a[t];i++) { if(a[t]%i==0) { cnt[i]--; ans-=(ll)mu[i]*cnt[i]; if(i*i!=a[t])cnt[a[t]/i]--,ans-=mu[a[t]/i]*cnt[a[t]/i]; } } } cout<<ans<<endl; } return 0; }
例题4：大白书P301石头染色方案计数
代码：
#include<iostream> #include<cstdio> #include<cstring> #include<vector> #include<map> using namespace std; #define fi first #define se second #define pi acos(-1.0) #define LL long long //#define mp make_pair #define pb push_back #define ls rt<<1, l, m #define rs rt<<1|1, m+1, r #define ULL unsigned LL #define pll pair<LL, LL> #define pli pair<LL, int> #define pii pair<int, int> #define piii pair<pii, int> #define mem(a, b) memset(a, b, sizeof(a)) #define fio ios::sync_with_stdio(false);cin.tie(0);cout.tie(0); #define fopen freopen("in.txt", "r", stdin);freopen("out.txt", "w", stout); //head const int mod = 1e9 + 7; LL q_pow(LL n, LL k) { LL ans = 1; while(k) { if(k&1) ans = (ans * n) % mod; n = (n * n) % mod; k >>= 1; } return ans; } map<int, int> get_mu(int n) { map<int, int> res; vector<int> p; for (int i = 2; i * i <= n; i++) { if(n % i == 0) { p.pb(i); while(n % i == 0) n /= i; } } if(n > 1) p.pb(n); int m = (int)p.size(); for (int i = 0; i < (1<<m); i++) { int t = 1, mu = 1; for (int j = 0; j < m; j++) { if(i & (1<<j)) { t *= p[j]; mu *= -1; } } res[t] = mu; } return res; } int main() { int n, m; while(~scanf("%d %d", &n, &m)) { LL ans = 0; for(int i = 1; i*i <= n; i++) { if(n % i == 0) { if(i * i == n) { map<int, int> mu = get_mu(i); LL t = 0; for (map<int, int>::iterator it = mu.begin(); it != mu.end(); it++) { t = (t + (it->se) * q_pow(m, i/(it->fi))) % mod; } t = (t * q_pow(i, mod-2)) % mod; ans = (ans + t) % mod; } else { int t1 = i, t2 = n/i; map<int, int> mu = get_mu(t1); LL t = 0; for (map<int, int>::iterator it = mu.begin(); it != mu.end(); it++) { t = (t + (it->se) * q_pow(m, t1/(it->fi))) % mod; } t = (t * q_pow(t1, mod-2)) % mod; ans = (ans + t) % mod; mu = get_mu(t2); t = 0; for (map<int, int>::iterator it = mu.begin(); it != mu.end(); it++) { t = (t + (it->se) * q_pow(m, t2/(it->fi))) % mod; } t = (t * q_pow(t2, mod-2)) % mod; ans = (ans + t) % mod; } } } printf("%lld\n", ans); } return 0; }
例题5：P4318 完全平方数
思路：二分+容斥
代码：
#pragma GCC optimize(2) #pragma GCC optimize(3) #pragma GCC optimize(4) #include<bits/stdc++.h> using namespace std; #define y1 y11 #define fi first #define se second #define pi acos(-1.0) #define LL long long //#define mp make_pair #define pb emplace_back #define ls rt<<1, l, m #define rs rt<<1|1, m+1, r #define ULL unsigned LL #define pll pair<LL, LL> #define pli pair<LL, int> #define pii pair<int, int> #define piii pair<pii, int> #define pdd pair<double, double> #define mem(a, b) memset(a, b, sizeof(a)) #define debug(x) cerr << #x << " = " << x << "\n"; #define fio ios::sync_with_stdio(false);cin.tie(0);cout.tie(0); //head const int N = 5e4 + 10; int prime[N], cnt, mu[N], T, k; bool not_p[N]; void seive() { mu[1] = 1; for (int i = 2; i < N; ++i) { if(!not_p[i]) prime[++cnt] = i, mu[i] = -1; for (int j = 1; i*prime[j] < N && j <= cnt; ++j) { not_p[i*prime[j]] = true; if(i%prime[j] == 0) { mu[i*prime[j]] = 0; break; } else mu[i*prime[j]] = -mu[i]; } } } bool ck(LL n) { LL tot = 0; for (LL i = 1; i*i <= n; ++i) { tot += mu[i]*(n/(i*i)); } return tot >= k; } int main() { seive(); scanf("%d", &T); while(T--) { scanf("%d", &k); LL l = 1, r = 2e9 + 100, m = l+r >> 1; while(l < r) { if(ck(m)) r = m; else l = m+1; m = l+r >> 1; } printf("%lld\n", m); } return 0; }
