算法笔记--莫比乌斯反演

这个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;
}

例题2：HDU 6053 TrickGCD

代码：

#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
",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
", 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 << "
";
#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
", m);
    }
    return 0;
}