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

最新推荐文章于 2024-09-01 16:58:11 发布

weixin_34378045

最新推荐文章于 2024-09-01 16:58:11 发布

阅读量72

点赞数

CC 4.0 BY-SA版权

原文链接：http://www.cnblogs.com/widsom/p/7701129.html

这个B站上的视频讲的不错：http://www.bilibili.com/video/av14325327/

定理：和是定义在非负整数集合上的两个函数，并且满足条件，那么我们得到结论

在上面的公式中有一个函数，它的定义如下：

（1）若，那么

（2）若，均为互异素数，那么

（3）其它情况下

对于函数，它有如下的常见性质：

（1）对任意正整数有

（2）对任意正整数有

线性筛选求莫比乌斯反演函数代码：

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

View Code

例题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\n",Case,ans);
    } 
    return 0;
}

View Code

例题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;
}

View Code

例题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;
}