莫比乌斯反演与杜教筛
积性函数
定义
对于一个数论函数 \(f\),若满足 \(\forall (a,b)=1\) 都有 \(f(ab)=f(a)f(b)\),那么称 \(f\) 为积性函数。
例子
常见的积性函数有很多,比如以下:
- \(id(x)=x\)。
- \(\epsilon(x)=[x=1]\)。
- \(d(x)=\sum_{i=1}[i|x]\)。
狄利克雷卷积
对于两个数论函数 \(f,g\) 定义他们的狄利克雷卷积
\[f*g(x)=\sum_{ij=x} f_ig_j \]这个玩意有个漂亮的性质,我们后面再谈。
【例题】
求 \(1*\varphi\)。
解答:原式相当于 \(\sum_{i|x} \varphi(i)\),下面介绍两种推法:
法1:展开法 考虑 \(x\) 的标准形式 \(x=p_1^{\alpha _1}p_2^{\alpha_2}\dots p_k^{\alpha_k}\),那么观察 \(\varphi(x)\) 的计算式,发现它 \(p_i^{\alpha_i}\) 的系数都在一下的集合中挑选:
\[\{1,p-1,(p-1)p,\dots,(p-1)p^{\alpha-1}\} \]不难发现上述集合所有元素的和恰好为 \(p^{\alpha}\),全部乘起来结果恰为 \(n\)。
法2:一一对应法 考虑所有 \(\leq x\) 的正整数 \(n\)。那么 \(\dfrac{n}{\gcd(n,x)}\) 必然和 \(\dfrac{x}{\gcd(n,x)}\) 互质,这就形成了一个一一对应,所以 \(1*\varphi=id\)。
性质
【性质1】 只需要知道所有 \(x\) 为质数的方幂时的值就可以确定整个函数,换句话说假设 \(x\) 的标准分解为 \(p_1^{\alpha _1}p_2^{\alpha_2}\dots p_k^{\alpha_k}\) 就会有:
\[f(x)=\prod f(p_i^{\alpha_i}) \]【性质2】 积性函数的狄利克雷卷积为积性函数。
证明:考虑所有的 \(x\) 为质数方幂 \(p^a\) 的值,\(f*g(p^a)=\sum_{i+j=a} f(p^i)g(p^j)\)。不难发现其他的 \(x\) 也可以可以写成标准分解的形式,证毕。
莫比乌斯函数 \(\mu\)
我们需要寻找一个函数 \(\mu\) 使得 \(1*\mu=\epsilon\)。
- 因为 \(1\) 为积性函数,所以 \(\mu\) 为积性函数。
- \(1*\mu(1)=1\),所以 \(\mu(1)=1\)。
- 假设 \(p\) 为质数,那么 \(1*\mu(p)=\mu(1)+\mu(p)=0\),所以 \(\mu(p)=-1\)。
- 假设 \(p\) 为质数,\(k\geq 2\)。那么 \(1*\mu(p)=\mu(1)+\mu(p)+\dots +\mu(p^k)\),不难得到 \(\mu(p^k)=0\)。
- 有了前面的条件可以知道:
给一个线性筛的 \(\mu\) 求法:
mu[1]=1;
for(int i=2;i<MAXN;i++){
if(!vis[i]) mu[i]=-1,p[++cnt]=i;
for(int j=1;j<=cnt&&p[j]*i<MAXN;j++){
vis[i*p[j]]=1;
if(i%p[j]==0) break;
mu[i*p[j]]=-mu[i];
}
}
莫比乌斯反演
简介
假设我们需要函数 \(f\) 不好求,但是 \(1*f=g\) 是好求的函数,我们就可以有 \(\mu*1*f=g*\mu\) 所以 \(\epsilon*f=g*\mu\),也就是 \(f(1)=g*\mu\)。\(f(2,3,\dots)\) 是同理的,直接暴力复杂度是 \(O(n\ln n)\) 的。
【习题】补全代码
- 给定数论函数 \(g\) 求 \(f\) 使得 \(1*f=g\)。
- \(n\leq 10^6\)。
应用
【例题1】P3455 [POI2007]ZAP-Queries 改编
- 求 \(\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)=1]\)。
- 多测,\(T,n,m\leq 5\times 10^4\)。
对 \(\gcd=1\) 问题一般使用莫反处理
考虑如下变换:
\[\begin{aligned}\sum_{i=1}^a\sum_{j=1}^b[\gcd(i,j)=1]&=\sum_{i=1}^a\sum_{j=1}^b\sum_{k|\gcd(i,j)}\mu(k)\\&=\sum_{k=1}\sum_{k|i}^a\sum_{k|j}^b\mu(k)\\&=\sum_{k=1}\ \mu(k)\times [\dfrac ak]\times[\dfrac bk]\end{aligned} \]这里第一步 \(\epsilon=1*\mu\) 的性质。对最后的式子数论分块可以做到 \(O(T\sqrt n)\)。
ll solve(){
ll ans=0;
for(int l=1,r;l<=min(a,b);l=r+1){
r=min(a/(a/l),b/(b/l));
ans+=(a/l)*(b/l)*(mu[r]-mu[l-1]);
}
return ans;
}
【例题2】P2257 YY的GCD
- 求 \(\sum_{i=1}^n\sum_{j=1}^m[\gcd(i,j)\in \textbf P]\)。
- 多测,\(T\leq 10^4,n,m\leq10^7\)。
我们还是像上一题一样向枚举 \(\gcd\) 得到这个奇怪式子:
\[\sum_{p\in\textbf P}\sum_{k=1}\ \mu(k)\times [\dfrac a{pk}]\times[\dfrac b{pk}] \]这个时间复杂度是 \(T\sqrt n\ln n\) 的,不可接受,考虑简化。
我们可以尝试枚举 \(pk\),这样子就变成了:
\[\begin{aligned} \sum_{T}\sum_{pk=T} \mu(k)\times [\dfrac aT]\times [\dfrac bT]&=\sum_T\sum_{k|T} \mu(k)\times [\dfrac aT]\times [\dfrac bT] \end{aligned} \]不难发现我们可以把 \(\sum_{k|T} \mu(k)\) 提前预处理出来,最后直接乘起来就好了,时间复杂度是 \(O(T\sqrt n)\) 可以接受。
#include<iostream>
using namespace std;
#define ll long long
const int MAXN=1e7+5;
ll mu[MAXN],pr[MAXN],cnt=0;bool vis[MAXN];
ll f[MAXN];
void init(){
mu[1]=1;
for(int i=2;i<MAXN;i++){
if(!vis[i]) pr[++cnt]=i,mu[i]=-1;
for(int j=1;j<=cnt&&i*pr[j]<MAXN;j++){
vis[i*pr[j]]=1;
if(i%pr[j]==0) break;
mu[i*pr[j]]=-mu[i];
}
}
for(int i=1;i<MAXN;i++)
for(int j=1;j<=cnt&&i*pr[j]<MAXN;j++)
f[i*pr[j]]+=mu[i];
for(int i=1;i<MAXN;i++)
f[i]+=f[i-1];
}
void solve(){
ll a,b;cin>>a>>b;
ll ans=0;
for(int l=1,r;l<=min(a,b);l=r+1){
r=min(a/(a/l),b/(b/l));
ans+=(a/l)*(b/l)*(f[r]-f[l-1]);
}
cout<<ans<<endl;
return;
}
int main(){
init();
int _;cin>>_;
while(_--) solve();
}
【习题】P2522 [HAOI2011]Problem b
【习题】P1829 [国家集训队]Crash的数字表格 / JZPTAB
【例题3】P3327 [SDOI2015]约数个数和
- 求 \(\sum_{i=1}^n\sum_{j=1}^md(ij)\)。
- 多测,\(T,n,m\leq5\times 10^4\)。
为了解决本题需要给一个有力的结论把 \(d\) 函数转化成 \(\gcd\) 函数:
\[d(xy)=\sum_{i|x}\sum_{j|y} [\gcd(i,j)=1] \]证明留给读者,一般来说解决数论等式的证明问题有两种手法,一种是直接展开,一种是构造一一映射。
后面的推导留给读者。
【习题】AT_agc038_c [AGC038C] LCMs
【习题】P3312 [SDOI2014]数表
【习题】P3704 [SDOI2017]数字表格
杜教筛
我们需要解决数论函数快速求前缀和问题,我们需要得到某一项的前缀和,伟大的杜老师发明了 \(O(n^{0.66})\) 筛出来的方法。
简单来说,对于给定函数 \(f\),我们需要找一个函数 \(g\) 使得 \(g,f*g\) 的前缀和都能 \(O(1)\)。那么考虑如下变换:
\[\sum_{i=1}^n f*g(i)=\sum_{i=1}^n\sum_{xy=i} f(x)g(y)=\sum_{i=1}^n (\sum_{j=1}^if(i))g([\dfrac ni]) \]假设 \(f\) 的前缀函数叫 \(S\),那么不难
\[\sum_{i=1}^n f*g(i)=g(1)S(n)+\sum_{i=2}^n S([\dfrac ni])g(i) \]分治,等式左边我们知道,最右边那一坨可以分治求,解决。时间复杂度:\(O(n^{0.66})\)。具体来说,我们要设置一个阈值 \(X\) 使得 \(\leq X\) 的前缀和暴力求,\(>X\) 的用 unordered_map
记录维护即可。
【例题】P4213 【模板】杜教筛(Sum)
- 求 \(\varphi,\mu\) 的前缀和。
- 多测,\(T\leq 10,n\leq 2^{31}\)。
不难发现:
- \(\varphi*1=id\)(上面例题)
- \(\mu *1=\epsilon\)(\(\mu\) 定义)
套公式即可。
const int MAXN=1e7+5;
ll phi[MAXN],mu[MAXN];bool vis[MAXN];
int p[MAXN],cnt=0;
int tmp=0;
unordered_map<int,ll> M1,M2;
void init(){
phi[1]=mu[1]=1;
for(int i=2;i<MAXN;i++){
if(!vis[i]) p[++cnt]=i,phi[i]=i-1,mu[i]=-1;
for(int j=1;j<=cnt&&i*p[j]<MAXN;j++){
vis[i*p[j]]=1;
if(i%p[j]==0){
phi[i*p[j]]=phi[i]*p[j];
break;
}
phi[i*p[j]]=phi[i]*(p[j]-1);
mu[i*p[j]]=-mu[i];
}
}
for(int i=1;i<MAXN;i++)
phi[i]+=phi[i-1],mu[i]+=mu[i-1];
}
ll pre_phi(ll x){
if(x<MAXN) return phi[x];
if(M1[x]) return M1[x];
ll ans=x*(x+1)/2;
for(ll l=2,r;l<=x;l=r+1){
r=x/(x/l);
ans-=(r-l+1)*pre_phi(x/l);
}
return M1[x]=ans;
}
ll pre_mu(ll x){
if(x<MAXN) return mu[x];
if(M2[x]) return M2[x];
ll ans=1;
for(ll l=2,r;l<=x;l=r+1){
r=x/(x/l);
ans-=(r-l+1)*pre_mu(x/l);
}
return M2[x]=ans;
}
int main(){
init();
int _;cin>>_;
while(_--){
ll n;cin>>n;
cout<<pre_phi(n)<<" "<<pre_mu(n)<<endl;
}
}
【习题】P3768 简单的数学题
标签:函数,dfrac,sum,杜教,mu,反演,莫比,times,gcd From: https://www.cnblogs.com/KawaiiOU/p/17121316.html