<学习笔记> 莫比乌斯反演

f(n)g(n) 是定义在整数集上的两个函数,满足一下关系:

f(n)=dng(n)

但是手摸发现不仅可以从 gf,还可以 fg,发现还满足这个关系:

g(n)=dnf(d)μ(nd)

μ(x)={1x=10(1)kk

性质一: 对于任意的正整数 n,有 dnμ(d)={1n=10n>1

证明

n=1 时,dnμ(d)=1,成立。

n>1 时,设 n=P1x1...pkxk ,设 d=P1y1...Pkyk 。因为 yi>2 全为 0,我们只考虑 yi0/1 的情况,假设 d 中有 rn 的互异质因子,且质因子质数位 1 ,那么对应 μ(d)=(1)r,这样的 dCkr 个,所以 dnμ(d)=r=0kCkr(1)r=(11)k=0,其实就是二项式定理。

性质二: μ(n) 是积性函数。

证明 易证。

莫比乌斯反演的证明#

形式一#

f(n)=dng(d)g(n)=dnf(nd)μ(d)

证明一

g(n)=dnf(nd)μ(d)=dnμ(d)indg(i)

将上式换成以 g(i) 为主体的形式,那么有:

ing(i)id|nμ(d)=ing(i)dniμ(d)

根据性质一,只有 i=n 时才有值,否则全为 0

综上 =g(n) 得证。

证明二

I(n)=1id(n)=nϵ(n)=[n=1]

单位元 ϵ 满足 fϵ=f

根据狄利克雷卷积得到几个性质:

  1. μI=ϵ

  2. φI=id

  3. μid=φ

给出的条件等价于 g=fI

根据性质一 gμ=fIμ=fϵ=f,得证。

变形#

f(i)=d=1nig(di)g(i)=d=1nif(di)μ(d)

证明

g(i)=d=1nif(di)μ(d)=d=1niμ(d)d1=1ndig(d1di)

T=dd1,则有上式

=T=1nig(Ti)dTμ(d)

根据性质一有,只有当 T=1 时,原式为 g(i) , 得证。

莫比乌斯反演的性质#

g(n) 是积性函数的充分必要条件是 f(n) 是积性函数。

莫比乌斯反演的应用#

g(i) 很难求,但是 f(i) 很好求,可以根据反演求 g(i)

f(i)=d=1nig(di)g(i)=d=1nif(di)μ(d)

P4449 于神之怒加强版 #

f(p) 表示 gcd(i,j)=p 的个数,那么 Ans=pf(p)pk

发现 f(p) 不好求,那么我们设 g(p) 表示 pgcd(i,j) 的个数,那么 g(p)=np×mp

g(p)=d=1npf(dp)

f(p)=d=1npμ(d)g(dp)

=d=1npμ(d)ndpmdp

那么

Ans=p=1npkd=1npμ(d)ndpmdp

T=dp 则有 Ans=T=1nnTmTpTμ(Tp)pk

那么我们只需要线性求出 pTμ(Tp)pk 就可以了。

g(x)=xk,设 G(T)=pTμ(Tp)g(p)

首先可以知道 G(x) 为积性函数,所以 G(ab)=G(a)G(b)gcd(a,b)=1

那么我们尝试证明 G(ab)=G(a)bkgcd(a,b)1

我们设 a=p1k1···pnkn,另 b=pttn

因为有 G(a)=i=1nG(piki)

那么就有

G(ab)G(a)=G(ptkt+1)G(ptkt)

G(ptkt)=i=0ctμ(pti)g(ptcti)

只有 i=0/1 时才有值,那么有

G(ptkt)=g(ptct)g(ptct1)

那么

G(ab)G(a)=g(ptct+1)g(ptct)g(ptct)g(ptct1)=ptk

那么我们就可以在线性筛的时候求出,然后计算答案直接整数分块,复杂度 O(Tn)

code
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=5*1e6+10;
const int mod=1e9+7;
int prime[N],nprime[N],mul[N],f[N],sumf[N];
int T,k;
int qpow(int x,int p){
    int ans=1;
    while(p){
        if(p&1) ans=(ans*x)%mod;
        x=(x*x)%mod;
        p>>=1;
    }
    return ans;
}
void get_prime(){
    mul[1]=1;
    f[1]=1;
    for(int i=2;i<=N-10;i++){
        if(!nprime[i]){
            f[i]=(qpow(i,k)-1+mod)%mod;       
            prime[++prime[0]]=i;
            mul[i]=-1;
        }
        for(int j=1;j<=prime[0] && i*prime[j]<=N-10;j++){
            int tmp=i*prime[j];
            nprime[tmp]=1;
            mul[tmp]=mul[i]*mul[prime[j]];
            if(i%prime[j]==0){
                f[tmp]=f[i]*qpow(prime[j],k)%mod;
                mul[tmp]=0;
                break;
            }
            f[tmp]=(f[i]*f[prime[j]]%mod)%mod;
        }
    }
    for(int i=1;i<=N-10;i++) sumf[i]=(sumf[i-1]+f[i])%mod;    
}
signed main(){
    scanf("%lld%lld",&T,&k);
    get_prime();
    while(T--){
		int n,m;
		scanf("%lld%lld",&n,&m);
		if(n>m) swap(n,m);
		int l=1,r=0;
		int ans=0;
		while(l<=n){
			int ta=n/l,tb=m/l;
			if(!ta||!tb) break;
			r=min(n/ta,m/tb);
			ans=(ans+ta*tb%mod*(sumf[r]-sumf[l-1]+mod)%mod)%mod;
			l=r+1;
		}
		printf("%lld\n",ans);
	}
}

DZY Loves Math #

首先它让求 i=1nj=1mF(gcd(i,j))

然后套路的转化为

x=1nF(x)d=1nxμ(d)nxdmxd

T=xd ,则有

T=1nnTmTxTF(x)μ(Tx)

现在只需要求出 G(T)=xTF(x)μ(Tx) 就可以整除分块了。

我们设 T=p1z1...pkzkx=p1y1...pkyk

我们设 a 为最大的 ziq 表示 zi 为最大值的个数,k 为质因子个数。

我们考虑 μ(Td)0,所以 yi=zi/zi1 , 那么 F(x)=a/a1

q=k

  • 若最大值为 F(d)=a,那么只要满足一个为 zi=yi,所以有 (补一位进行二项式定理)

F(d)=aF(d)μ(Td)=aF(d)=aμ(Td)=ai=1k(1)kiCki=a(i=0k(1)kiCki(1)k)=a(1)k

  • 若最大值为 F(d)=a1,就有

F(d)=a1f(d)μ(Td)=(a1)(1)k

那么这部分答案为 g(T)=a(1)k+(a1)(1)k=(1)k+1

q<k

这时不论 F(d)=a 还是 F(d)=a1,都需要考虑另外 kq 个元素也就是在上式乘上个 j=0kq(1)jCkqj 发现这部分为零,所以这部分的答案为 g(T)=0

可以用线性筛来求,发现每个数会被最小的质数筛掉,那么我们维护一个 mnsum(x) 表示 x 的最小质因子个数, excp(x) 表示除了最小质因子之外的数,然后求得时候维护。

code
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=1e7+10;
int g[N],mnsum[N],excp[N];
int prime[N];
bool nprime[N];
void get_prime(){
    for(int i=2;i<=N-5;i++){
        if(!nprime[i]){
            prime[++prime[0]]=i;
            g[i]=mnsum[i]=excp[i]=1;
        }
        for(int j=1;j<=prime[0] && i*prime[j]<=N-5;j++){
            int k=prime[j]*i;
            nprime[k]=1;
            if(i%prime[j]==0){
                excp[k]=excp[i];
                mnsum[k]=mnsum[i]+1;
                if(excp[k]==1) g[k]=1;
                else g[k]=(mnsum[excp[k]]==mnsum[k]?-g[excp[k]]:0);
                break;
            }
            excp[k]=i,mnsum[k]=1;
            g[k]=(mnsum[i]==1?-g[i]:0);
        }
    }
    for(int i=1;i<=N-5;i++) g[i]+=g[i-1];
}
signed main(){
    get_prime();
    int T;
    scanf("%lld",&T);
    while(T--){
        int a,b;
        scanf("%lld%lld",&a,&b);
        if(a>b) swap(a,b);
        int l=1,r=0;
        int ans=0;
        while(l<=a){
            int ta=a/l,tb=b/l;
            if(!ta || !tb) break;
            int r=min(a/ta,b/tb);
            ans+=ta*tb*(g[r]-g[l-1]);
            l=r+1;
        }
        printf("%lld\n",ans);
    }
}

[SDOI2015]约数个数和#

首先给出结论 d(nm)=injm[gcd(i,j)=1]

证明

nm=p1a1...pkak,那么有 d(nm)=i=1kai+1

n=p1b1...pkbkm=p1a1b1...pkakbk

i=p1x1...pkxkj=p1y1...pkyk

要想 gcd(i,j)=1,必须满足 xi=0yi=0

  • x1=0 时,那么 y1 可以取 [0,a1b1]

  • y1=0 时,那么 x1 可以取 [0,b1]

那么对于第一个质因子,符合条件的组合就有 a1b1+1+b1+11=a1+1

以此类推,那么所有符合条件的组合就有 i=1i=kai+1

因此 injm[gcd(i,j)=1]=i=1i=kai+1=d(nm)

回归正题那么原式就为

i=1nj=1mxiyj[gcd(x,y)=1]

x=1ny=1mnxmy[gcd(x,y)=1]

f(x)=x=1ny=1mnxmy[gcd(x,y)=x]

g(x)=x=1ny=1mnxmy[x|gcd(x,y)]

因为 ndx=ndxx ,设 s(n)=i=1nni,那么 g(d)=s(nd)s(md)

根据莫比乌斯反演,就有

f(x)=d=1ng(xd)μ(d)

因为我们只需要求 $f(1),那么答案就是 f(1)=d=1ng(d)μ(d)

整除分块就可以了。

code
//  P3327 [SDOI2015] 约数个数和 
#include<bits/stdc++.h>
#define int long long
using namespace std;
const int N=50005;
int mul[N],suml[N];
int prime[N];
bool nprime[N];
void get_prime(){
    mul[1]=1;
    for(int i=2;i<=N-3;i++){
        if(!nprime[i]){
            prime[++prime[0]]=i;
            mul[i]=-1;
        }
        for(int j=1;j<=prime[0] && i*prime[j]<=N-3;j++){
            int k=prime[j]*i;
            nprime[k]=1;
            if(i%prime[j]==0){
                mul[k]=0;
                break;
            }
            mul[k]=mul[prime[j]]*mul[i];
        }
    }
    for(int i=1;i<=N-3;i++){
        suml[i]=suml[i-1]+mul[i];
    }
}
int s[N];
signed main(){
    get_prime();
    int T;
    scanf("%lld",&T);
    for(int n=1;n<=N-3;n++){
        int l=1,r=0;
        while(l<=n){
            int t=n/l;
            if(!t) break;
            int r=n/t;
            s[n]+=(r-l+1)*t;
            l=r+1;
        }
    }
    while(T--){
        int n,m;
        scanf("%lld%lld",&n,&m);
        if(n>m)  swap(n,m);
        int ans=0;
        int l=1,r=0;
        while(l<=n){
            int ta=n/l,tb=m/l;
            if(!ta || !tb) break;
            int r=min(n/ta,m/tb);
            ans+=s[ta]*s[tb]*(suml[r]-suml[l-1]);
            l=r+1;
        }
        printf("%lld\n",ans);
    }
}

作者:bloss

出处:https://www.cnblogs.com/jinjiaqioi/p/17627668.html

版权:本作品采用「署名-非商业性使用-相同方式共享 4.0 国际」许可协议进行许可。

posted @   _bloss  阅读(41)  评论(0编辑  收藏  举报
相关博文:
阅读排行:
· 单线程的Redis速度为什么快?
· 展开说说关于C#中ORM框架的用法!
· Pantheons:用 TypeScript 打造主流大模型对话的一站式集成库
· SQL Server 2025 AI相关能力初探
· 为什么 退出登录 或 修改密码 无法使 token 失效
点击右上角即可分享
微信分享提示
more_horiz
keyboard_arrow_up dark_mode palette
选择主题
menu