ICode9

精准搜索请尝试: 精确搜索
首页 > 其他分享> 文章详细

莫比乌斯反演

2020-12-16 17:34:15  阅读:278  来源: 互联网

标签:lfloor frac cdot 乌斯 sum rfloor mu 反演 莫比


莫比乌斯函数定义:

设 \(n = p_1 ^ {k_1} \cdot p_2 ^ {k_2} \cdot\cdots\cdot p_m ^ {k_m}\),其中 p 为素数,则定义如下:

\[\mu(n) = \begin{cases} 1 & n = 1 \\ (-1) ^ m & \prod\limits_{i = 1} ^ {m} k_i = 1 \\ 0 & \textrm{otherwise}(k_i \gt 1) \end{cases} \]

​​
可知有平方因子的数的莫比乌斯函数值为 0。

  1. 由于Mobius函数是积性函数,即 \(\mu(a \cdot b)=\mu(a)\cdot \mu(b), (a,b)=1\),所以
    可以在\(O(n)\)的时间内筛出\([1,n]\)的函数值,一个数 \(x= p_1 ^ {k_1} \cdot p_2 ^ {k_2} \cdot\dots\cdot p_m ^ {k_m}\),
    若\(k_i=1 , k_j>1, j>i\) 则x只会被\(\frac x {p_j} \mod p_j=0\) 删除
int mu[MAXN],prime[MAXN],cnt;
bool noprime[MAXN];
void mublus(){
    mu[1]=1;
    for(int i=2;i<MAXN;++i){
        if(!noprime[i]){
            prime[++cnt]=i;
            mu[i]=-1;//质数为-1
        }
        for(int j=1;j<=cnt&&prime[j]*i<MAXN;++j){
            noprime[prime[j]*i]=1;
            if(i%prime[j]==0){
                mu[prime[j]*i]=0;// prime[j]*i有因子prime[j]^2   ,故为0;
                break;
            }
            else mu[prime[j]*i]=-mu[i];//   i 比prime[j]*i多了一个质因子。
        }
    }
}
  1. \(\sum_{d\mid n}\mu(d) = [n = 1]\)

证明:当 n = 1 时有 \(\mu(1) = 1\),显然成立。否则设 \(n = p_1 ^ {k_1} \cdot p_2 ^ {k_2}\cdot\cdots\cdot p_m ^ {k_m}\)​​ ,\(d = p_1 ^ {x_1} \cdot p_2 ^ {x_2} \cdot\cdots\cdot p_m ^ {x_m}\)

​​ 根据\(\mu(n)\) 的定义,只需考虑 \(x_i = 0\) 或 \(x_i = 1\)的情况。我们设 d 中存在 r 个 \(x_i\)为 1,那么有:

\(\sum_{d\mid n}\mu(d) = \sum_{r = 0}^m\binom m r(-1)^r(n \neq 1)\)

根据二项式定理:

\((x+y)^{n}=\sum _{k=0}^{n}\binom n kx^{n-k}y^{k} (x+y)\)
​​
我们令 x = 1, y = -1,即得证:

\(\sum_{r = 0}^m\binom m r(-1)^r = (1 - 1)^m = 0(n \neq 1)\)

莫比乌斯反演

假设对于数论函数 \(f(n)\) 和\(F(n)\),满足以下关系式:\(F(n)=\sum _{{d|n}}f(d)\)

则将其默比乌斯反演公式定义为:

\[f(n)=\sum _{d|n}\mu (d) F({\frac {n}{d}}) \]

或者满足:\(F(n)=\sum_{n|d}f(d)\)

\[f(n)=\sum_{n|d}\mu(\frac d n)F(d) \]

1.求区间[1,a], [1,b]中gcd=k的个数

\[\sum_{i=1}^{a}\sum_{j=1}^{b}[gcd(i,j)=k] \]

即求

\[\sum_{i=1}^{\lfloor\frac a k\rfloor}\sum_{j=1}^{\lfloor\frac b k\rfloor}[gcd(i,j)=1] \]

令\(f(k)\)为满足gcd=k的个数,\(F(k)\)为满足\(gcd=x \cdot k,x>0\) 的个数.

则有 \(F(n)=\sum_{n|d}f(d)\) (d是n的倍数,\(\sum f(d)\)显然就是F的定义)

所以有反演公式:\(f(n)=\sum_{n|d}\mu(\frac d n)F(d)\)

显然,\(F(n)=\lfloor\frac a n\rfloor\cdot\lfloor\frac b n\rfloor\)

(a/n是第一个区间里n的倍数的个数,两者相乘则为 i,j 都是n的倍数的个数,所以gcd(i,j)是n的倍数,显然就是F的定义)

所以\(f(n)=\sum_{n|d}\mu(\frac d n)F(d)=\sum_{n|d}\mu(\frac d n)\cdot\lfloor\frac a d\rfloor\cdot\lfloor\frac b d\rfloor\)

由于k=1

有\(f(1)=\sum_{1|d}\mu(d)F(d)=\sum_{1|d}\mu(d)\cdot\lfloor\frac a d\rfloor\cdot\lfloor\frac b d\rfloor\)

HDU-1695

求[a,b],[c,d]里gcd=k的个数,a=c=1,gcd(1,2)和gcd(2,1)算作一个。

\(ans=f(1)=\sum_{1|i}^{\min(b,d)}\mu(i)F(i)=\sum_{1|i}\mu(i)\cdot\lfloor\frac b i\rfloor\cdot\lfloor\frac d i\rfloor\)
为何sum的上界是min(b,d),因为i大于min时,F(i)=0,且F(i)表示两个区间里各取一个数,满足gcd=i的倍数,若i大于其中一个区间的上界,则不能满足gcd=i的倍数。
最后,ans是有重复的,例如[1,5]和[1,10]两个区间,[1,3]和[3,1]算作两次,但是[1,6],[6,1]中[6,1]不可能倍计算,所以只算了一次。所以,在[1,min(b,d)]里的答案都计算了两次。

#include<bits/stdc++.h>
using namespace std;
#define Init(arr,val) memset(arr,val,sizeof(arr))
const int inf=0x3f3f3f3f,mod=1000000007,MAXN=1e5+8;
typedef long long ll;
int mu[MAXN],prime[MAXN],cnt;
bool noprime[MAXN];
void mubius(){
    mu[1]=1;
    for(int i=2;i<MAXN;++i){
        if(!noprime[i]){
            prime[++cnt]=i;
            mu[i]=-1;
        }
        for(int j=1;j<=cnt&&prime[j]*i<MAXN;++j){
            noprime[prime[j]*i]=1;
            if(i%prime[j]==0){
                mu[prime[j]*i]=0;
                break;
            }
            else mu[prime[j]*i]=-mu[i];
        }
    }
}
int main() {
    std::ios::sync_with_stdio(0);
    std::cin.tie(0);
    std::cout.tie(0);
    mubius();
    int T,a,b,c,d,k;
    cin>>T;
    for(int Case=1;Case<=T;++Case){
        cin>>a>>b>>c>>d>>k;
        if(k==0){//没有gcd=0
            cout<<"Case "<<Case<<": 0\n";
            continue;
        }
        b/=k,d/=k;//将gcd=k化成gcd=1.
        ll ans1=0,ans2=0;
        if(b>d)swap(b,d);//b,d里b是最小值。
        for(int i=1;i<=b;++i){
        	ans1+=1ll*mu[i]*(b/i)*(d/i);
        	ans2+=1ll*mu[i]*(b/i)*(b/i);
        }
        cout<<"Case "<<Case<<": "<<ans1-(ans2>>1)<<endl;
    }
    return 0;
}

这样做是 O(n)。
例如,b=10000,d=11000,i>6000时,\(\lfloor\frac b i\rfloor=\lfloor\frac d i\rfloor=1\),但是我们还是计算了4000次。

所以,不如先求\(\mu\) 的前缀和,再求出对于b,d两数,\(\lfloor\frac b i\rfloor, \lfloor\frac d i\rfloor\)均不变的范围,在该范围内,

\(ans=sum(\mu)\cdot\lfloor\frac b i\rfloor\cdot\lfloor\frac d i\rfloor\)

对于一个在区间[1 ,x]内的位置 i,设 \(\lfloor\frac x i\rfloor=k\) ,对于向下取整都为k的最后一个位置last,

\(\lfloor\frac x {last}\rfloor=k\),且 \(\lfloor\frac {x} {last+1}\rfloor=k-1,x\mod ({last+1})=0\)

得:\(last=\frac{x}{\lfloor\frac x i\rfloor -1}-1\),因为x和\(\lfloor\frac x i\rfloor -1\)是整除关系,把-1去掉,分子相当于+1,所以\(\frac{x}{\lfloor\frac x i\rfloor }-1\)不是整除,

比原来少了1,再加上1就和原来一样了,所以为了方便,减1全部去掉。

最后再和另一个区间里的取最小值。

for(int i=1;i<=b;++i){
        	ans1+=1ll*mu[i]*(b/i)*(d/i);
        	ans2+=1ll*mu[i]*(b/i)*(b/i);
	}
//由上面变成了下面,只需处理一下前缀和:
//for(int i=1;i<MAXN;++i)sumu[i]=sumu[i-1]+mu[i];
for(int last,i=1;i<=b;i=last+1){
            last=min(b/(b/i),d/(d/i));
            ans1+=1ll*(sumu[last]-sumu[i-1])*(b/i)*(d/i);
            ans2+=1ll*(sumu[last]-sumu[i-1])*(b/i)*(b/i);
        }

复杂度:\(O(\sqrt n + \sqrt m)\)

P2257 YY的GCD

题意:给定N, M,求1<=x<=N, 1<=y<=M且gcd(x, y)为质数的(x, y)有多少对。

还是一样,设\(f(d)=\sum_{i=1}^N \sum_{j=1}^M[gcd(i,j)=d]\) 假设\(N\le M\)

\(F(n)=\sum_{i=1}^N \sum_{j=1}^M[n|gcd(i,j)]=\lfloor\frac N n\rfloor\lfloor\frac M n\rfloor\)

由反演公式得到:

\(f(n)=\sum_{n|d}\mu(\frac d n)F(d)=\sum_{n|d}\mu(\frac d n)\lfloor\frac N d\rfloor\lfloor\frac M d\rfloor\)

\(ans=\sum_{p\in prime}f(p)=\sum_{p\in prime}\sum_{p|d}\mu(\frac d p)\lfloor\frac N d\rfloor\lfloor\frac M d\rfloor\)

\(=\sum_{p\in prime}\sum_{k=1}\mu(k)\lfloor\frac N {dp}\rfloor\lfloor\frac M {kp}\rfloor, d:=kp\)

令\(T=kp\),有:

\(ans=\sum_{p\in prime}\sum_{k=1}\mu(k)\lfloor\frac N {T}\rfloor\lfloor\frac M {T}\rfloor=\sum_{p\in prime}\sum_{p|T}\mu(\frac T p)\lfloor\frac N {T}\rfloor\lfloor\frac M {T}\rfloor\)

重点来了,如何交换两个求和符号的顺序?画图大法好。

\(p=p_1\)时,\(T_1=p_1,2p_1,3p_1,\dots,n_1p_1\)
\(p=p_2\)时,\(T_2=p_2,2p_2,3p_2,\dots,n_2p_2\)
\(\dots\)
\(p=p_n\)时,\(T_n=p_n,2p_n,3p_n,\dots,np_n\)
其中\(T_i\le N\)
所以,T从最小的质数\(p_1\)开始枚举,由于任意合数 \(x\),都能分解为\(x=k_ip_i=k_jp_j\),所以T遍历\(p_1\)开始到x的

所有数,\(\frac T p\)枚举的是\(k_i \dots\),所以有:

\(ans=\sum_{p\in prime}\sum_{p|T}\mu(\frac T p)\lfloor\frac N {T}\rfloor\lfloor\frac M {T}\rfloor\)

\(=\sum_{T=1}^{x}\sum_{p|T, p\in prime}\mu(\frac T p)\lfloor\frac N {T}\rfloor\lfloor\frac M {T}\rfloor\)

\(=\sum_{T=1}^{x}\lfloor\frac N {T}\rfloor\lfloor\frac M {T}\rfloor\sum_{p|T, p\in prime}\mu(\frac T p)\)(T从2开始也是一样的,T=1时后面的求和等于0)

右边的求和可以通过预处理获得前缀和,对于\(\sum_{p|T, p\in prime}\mu(\frac T p)\),我们考虑每个质数p对T=kp的贡献。

#include<bits/stdc++.h>
using namespace std;
#define Init(arr,val) memset(arr,val,sizeof(arr))
const int inf=0x3f3f3f3f,mod=10007,MAXN=1e7+8;
typedef long long ll;
int mu[MAXN],prime[MAXN],cnt;
bool noprime[MAXN];
int sumu[MAXN];
void Mobius() {
   mu[1]=1;
   for(int i=2; i<MAXN; ++i) {
       if(!noprime[i]) {
   			prime[++cnt]=i;
   			mu[i]=-1;
       }
       for(int j=1; j<=cnt&&prime[j]*i<MAXN; ++j) {
           noprime[prime[j]*i]=1;
           if(i%prime[j]==0) {
               mu[i*prime[j]]=0;
               break;
           } else
               mu[i*prime[j]]=-mu[i];
       }
   }
   for(int i=1; i<=cnt; ++i)
       for(int j=1; j*prime[i]<MAXN; ++j)
           sumu[j*prime[i]]+=mu[j];
   for(int i=1; i<MAXN; ++i)
       sumu[i]+=sumu[i-1];
}
int main() {
   std::ios::sync_with_stdio(0);
   std::cin.tie(0);
   std::cout.tie(0);
   Mobius();
   int t;
   cin>>t;
   while(t--){
       int x,y,last;
       cin>>x>>y;
       ll ans=0;
       if(x>y)swap(x,y);
       for(int i=1;i<=x;i=last+1){
           last=min(x/(x/i),y/(y/i));
           ans+=1ll*(x/last)*(y/last)*(sumu[last]-sumu[i-1]);
       }
       cout<<ans<<endl;
   }
   return 0;
}

标签:lfloor,frac,cdot,乌斯,sum,rfloor,mu,反演,莫比
来源: https://www.cnblogs.com/foursmonth/p/14145116.html

本站声明: 1. iCode9 技术分享网(下文简称本站)提供的所有内容,仅供技术学习、探讨和分享;
2. 关于本站的所有留言、评论、转载及引用,纯属内容发起人的个人观点,与本站观点和立场无关;
3. 关于本站的所有言论和文字,纯属内容发起人的个人观点,与本站观点和立场无关;
4. 本站文章均是网友提供,不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属;如您发现该文章侵犯了您的权益,可联系我们第一时间进行删除;
5. 本站为非盈利性的个人网站,所有内容不会用来进行牟利,也不会利用任何形式的广告来间接获益,纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。

专注分享技术,共同学习,共同进步。侵权联系[81616952@qq.com]

Copyright (C)ICode9.com, All Rights Reserved.

ICode9版权所有