初涉莫比乌斯反演

今天我们来讨论莫比乌斯反演。我承认,反演这个东西对于数学不好的人来说确实很痛苦(比如我)。但是真正学透了,还是会发现这个东西非常巧妙。


预备知识

数论分块:

关于数论分块,我写过一篇博客,也介绍了一些例题,这里再做一个简介。
比如我们要求式子 nd=1nd ∑ d = 1 n ⌊ n d ⌋ ,那么其实我们会发现对于一段区间 i[l,r] i ∈ [ l , r ] ni ⌊ n i ⌋ 的值相同。所以我们枚举这个 l,r l , r ,就可以在 n n 的复杂度内求出这段区间值了。


莫比乌斯函数:

莫比乌斯函数 μ(x) μ ( x ) 的定义是:
1.当 d=1 d = 1 时, μ(d)=1 μ ( d ) = 1
2.当 d=Πki=1pi d = Π i = 1 k p i 其中 pi p i 为互不相同的质数,那么 μ(d)=1k μ ( d ) = − 1 k (也就是 d d 分解质因数后全是一次项,那么μ(d)为-1的项数次)
3.当 d d 含有任何质因子的次数大于2,则μ(d)=0

并且这里要介绍一些关于莫比乌斯函数的重要性质:
n|dμ(n)=[d==1] ∑ n | d μ ( n ) = [ d == 1 ] [v] [ v ] (v是一个语句)为当该语句成立时, [v]=1 [ v ] = 1 ,否则 [v]=0 [ v ] = 0
这个性质可以证明,首先 d=1 d = 1 ,那么只有一个 μ(1) μ ( 1 ) 所以值是1。否则我们将 n n 分解质因数,假设其有k个质因数,那么其实上式为 C0kC1k+C2k+Ckk C k 0 − C k 1 + C k 2 − … … + C k k (k为奇数则最后一项为负)。那么根据二项式定理,原式为0。(原式= (1+(1))k ( 1 + ( − 1 ) ) k )

n|dμ(n)n=φ(d)d ∑ n | d μ ( n ) n = φ ( d ) d 这个定理十分奇妙,然而遗憾的是我并不是特别会证。

接下去我们要讲怎么求出莫比乌斯函数,由它是积性函数的性质我们可以得到:莫比乌斯函数可以由线性筛得到,代码如下:

for(i=2;i<=MAXN;i++){
        if(!vis[i]) pri[++top]=i,mu[i]=-1;
        for(j=1;j<=top&&pri[j]*i<=MAXN;j++){
            vis[pri[j]*i]=1;
            if(i%pri[j]==0) break;
            mu[i*pri[j]]=-mu[i];
        }
    }

莫比乌斯反演

解决完莫比乌斯函数后,我们就可以开始反演了。首先,我应该不合时宜的直接指出莫比乌斯反演的内容:
定理:对于两个函数 F(n) F ( n ) f(n) f ( n ) 是定义在正整数集上的两个函数,并且满足条件:
F(n)=d|nf(d) F ( n ) = ∑ d | n f ( d ) 那么 f(n)=d|nμ(d)F(nd) f ( n ) = ∑ d | n μ ( d ) ∗ F ( n d )
同样也有若 F(n)=n|df(d) F ( n ) = ∑ n | d f ( d ) 那么 f(n)=n|dμ(dn)F(d) f ( n ) = ∑ n | d μ ( d n ) ∗ F ( d )
接下来我们来证明一下,这里就证明第一条,因为第二条同理。

d|nμ(d)F(nd)=d|nμ(d)d|ndf(d) ∑ d | n μ ( d ) F ( n d ) = ∑ d | n μ ( d ) ∑ d ′ | n d f ( d ′ )
我们设 nd=kd n d = k d ′ ,则 d=nkd d = n k d ′ ,所以 d|nd d | n d ′ ,因此我们可以将上面的式子改写:
d|nμ(d)d|ndf(d)=d|nf(d)d|ndμ(d) ∑ d | n μ ( d ) ∑ d ′ | n d f ( d ′ ) = ∑ d ′ | n f ( d ′ ) ∑ d | n d ′ μ ( d )
那么我们根据性质1,仅当 nd==1 n d ′ == 1 时, d|ndμ(d)=1 ∑ d | n d ′ μ ( d ) = 1 ,其它时候其都等于0。于是原式就等同于 f(n) f ( n )
所以 f(n)=d|nμ(d)F(nd) f ( n ) = ∑ d | n μ ( d ) F ( n d ) ,证明完毕。


例题

今天实在是写不动了,例题明天再补上。

好吧一天等于7月16号到8月7号,我承认我太懒了,本来打算17号写的,结果17号颓了一天,7月18号到8月2号去美国,一直也没写,于是就一直拖到今天。

总共四道例题

BZOJ2301

BZOJ2301[HAOI2011]problem b                       比较模板的题,也是我做的第一道,只要容斥一下即可。

#include<bits/stdc++.h>
using namespace std;
int read(){
    char c;int x;while(c=getchar(),c<'0'||c>'9');x=c-'0';
    while(c=getchar(),c>='0'&&c<='9') x=x*10+c-'0';return x;
}
void print(int x){
    if(x/10) print(x/10);
    putchar(x%10+'0');
}
int T,a,b,c,d,k,top,ans,mu[50005],sum[50005],pri[50005],vis[50005];
int calc(int n,int m){
    if(n>m) swap(n,m);
    int l=1,r=0,res=0;n=n/k;m=m/k;
    while(l<=n){
        r=min(n/(n/l),m/(m/l));r=min(r,n);
        res+=(n/l)*(m/l)*(sum[r]-sum[l-1]);
        l=r+1;
    }
    return res;
}
int main()
{
    T=read();mu[1]=sum[1]=1;
    for(int i=2;i<=50000;i++){
        if(!vis[i]) pri[++top]=i,mu[i]=-1;
        for(int j=1;j<=top&&pri[j]*i<=50000;j++){
            vis[i*pri[j]]=1;
            if(i%pri[j]==0) break;
            mu[i*pri[j]]=-mu[i];
        }
        sum[i]=sum[i-1]+mu[i];
    }
    while(T--){
        a=read();b=read();c=read();d=read();k=read();
        ans=calc(b,d)-calc(a-1,d)-calc(b,c-1)+calc(a-1,c-1);
        print(ans);puts("");
    }
    return 0;
}

BZOJ1101

BZOJ1101[POI2007]Zap                       比上一道还简单。

#include<bits/stdc++.h>
#define MAXN 50000
using namespace std;
int read(){
    char c;int x;while(c=getchar(),c<'0'||c>'9');x=c-'0';
    while(c=getchar(),c>='0'&&c<='9') x=x*10+c-'0';return x;
}
void print(int x){
    if(x/10) print(x/10);
    putchar(x%10+'0');
}
int top,T,a,b,k,mu[MAXN+5],sum[MAXN+5],pri[MAXN+5],vis[MAXN+5];
int calc(){
    int l=1,r=0,res=0;
    a=a/k;b=b/k;
    while(l<=a){
        r=min(a/(a/l),b/(b/l));r=min(r,a);
        res+=(sum[r]-sum[l-1])*(a/l)*(b/l);
        l=r+1;
    }
    return res;
}
int main()
{
    T=read();mu[1]=sum[1]=1;
    for(int i=2;i<=MAXN;i++){
        if(!vis[i]) pri[++top]=i,mu[i]=-1;
        for(int j=1;j<=top&&pri[j]*i<=MAXN;j++){
            vis[pri[j]*i]=1;
            if(i%pri[j]==0) break;
            mu[i*pri[j]]=-mu[i];
        }
        sum[i]=sum[i-1]+mu[i];
    }
    while(T--){
        a=read();b=read();k=read();
        if(a>b) swap(a,b);
        print(calc());puts("");
    }
    return 0;
}

LUOGU2257

Luogu P2257 YY的GCD                       这道题求的是质数,所以在原来的基础上还要再推一步,并要预处理。

#include<bits/stdc++.h>
#define MAXN 10000000
#define ll long long
using namespace std;
ll read(){
    char c;ll x;while(c=getchar(),c<'0'||c>'9');x=c-'0';
    while(c=getchar(),c>='0'&&c<='9') x=x*10+c-'0';return x;
}
void print(ll x){
    if(x/10) print(x/10);
    putchar(x%10+'0');
}
int T,n,m,top,l,r,mu[MAXN+5],pri[MAXN+5],vis[MAXN+5];
ll ans,sum[MAXN+5],g[MAXN+5];
int main()
{
    T=read();mu[1]=1;
    register int i,j;
    for(i=2;i<=MAXN;i++){
        if(!vis[i]) pri[++top]=i,mu[i]=-1;
        for(j=1;j<=top&&pri[j]*i<=MAXN;j++){
            vis[pri[j]*i]=1;
            if(i%pri[j]==0) break;
            mu[i*pri[j]]=-mu[i];
        }
    }
    for(j=1;j<=top;j++)
     for(i=1;i*pri[j]<=MAXN;i++) g[i*pri[j]]+=(ll)mu[i];
    for(i=1;i<=MAXN;i++) sum[i]=sum[i-1]+g[i];
    while(T--){
        n=read();m=read();ans=0;
        if(n>m) swap(n,m);l=1;
        while(l<=n){
            r=min(n/(n/l),m/(m/l));r=min(r,n);
            ans+=1ll*(n/l)*(m/l)*(sum[r]-sum[l-1]);
            l=r+1;
        }
        print(ans);puts("");
    }
    return 0;
}
BZOJ2005

BZOJ2005 [Noi2010]能量采集                       这道题由于是求值,所以要做两个数论分块,并相互套起来

#include<bits/stdc++.h>
#define MAXN 100005
#define ll long long
using namespace std;
int read(){
    char c;int x;while(c=getchar(),c<'0'||c>'9');x=c-'0';
    while(c=getchar(),c>='0'&&c<='9') x=x*10+c-'0';return x;
}
ll n,m,l,r,prime[MAXN],mu[MAXN],vis[MAXN],top,gcd,sum[MAXN];
ll calc(int x){
    ll l=1,ans=0,N=n/x,M=m/x;
    while(l<=N){
        ll r=min(N/(N/l),M/(M/l));r=min(r,N);
        ans+=(sum[r]-sum[l-1])*(N/l)*(M/l);
        l=r+1;
    }
    return ans;
}
int main()
{
    n=read();m=read();
    if(n>m) swap(n,m);mu[1]=sum[1]=1;
    for(int i=2;i<=m;i++){
        if(!vis[i]) prime[++top]=i,mu[i]=-1;
        for(int j=1;j<=top&&i*prime[j]<=m;j++){
            vis[i*prime[j]]=1;
            if(i%prime[j]==0) break;
            mu[i*prime[j]]=-mu[i];
        }
        sum[i]=sum[i-1]+mu[i];
    }
    l=1;
    while(l<=n){
        r=min(n/(n/l),m/(m/l));r=min(r,n);
        gcd+=calc(l)*(r-l+1)*(l+r)/2;
        l=r+1;
    }
    printf("%lld",2*gcd-n*m);
    return 0;
}

坑终于补完了

评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包
实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值