「学习笔记」莫比乌斯反演

wdgm4 · 2023-09-10 18:41:46 · 个人记录

开新坑了。QWQ

更好的阅读体验？

前置芝士：数论分块。（之后再说。QWQ）

积性函数

定义

一个数论函数 f(n) 满足 f(xy)=f(x) \times f(y)（\gcd(x,y)=1），则称 f(n) 是积性函数。

莫比乌斯函数：

杜教筛

杜教筛代码真是好短啊。

一种时间复杂度为 O(n^{\frac{2}{3}}) 的筛，可以求出一些积性函数的前缀和。

具体来说，求 S(n)=\sum\limits_{i=1}^{n} f(i)，f 是积性函数。

构造函数 h,g 满足 h=f*g。

那么：

那么：

（这里就是把 d=1 的给提出来了）

对于 \mu 函数，知道 \mu * I = \epsilon，而 \epsilon 和 \mu 的前缀和都很好求，所以可以用杜教筛。

而对于 \varphi 函数，知道 \varphi * I = id，也可以用杜教筛。

P4213 【模板】杜教筛（Sum）

代码：

#include<bits/stdc++.h>
#define XD 114514
#define MAXN 10000010
#define int long long
using namespace std;
int T,n;
int phi[MAXN],mu[MAXN],sumphi[MAXN],summu[MAXN];
bool vis[MAXN];
vector<int> prm;
map<int,int> mphi,mmu;
void sieve(){
    mu[1]=phi[1]=sumphi[1]=summu[1]=1;
    for(int i=2;i<=1e7;i++){
        if(!vis[i]){
            prm.push_back(i);
            phi[i]=i-1;mu[i]=-1;
        }
        sumphi[i]=sumphi[i-1]+phi[i];
        summu[i]=summu[i-1]+mu[i];
        for(auto j:prm){
            if(i*j>1e7) break;
            vis[i*j]=true;
            if(i%j==0){
                phi[i*j]=phi[i]*j;
                break;
            }
            phi[i*j]=phi[i]*phi[j];
            mu[i*j]=-mu[i];
        }
    }
}
int getmu(int x){
    if(x<=1e7) return summu[x];
    if(mmu[x]) return mmu[x];
    int nem=1,r;
    for(int l=2;l<=x;l=r+1){
        r=x/(x/l);
        nem-=(r-l+1)*getmu(x/l);
    }
    return mmu[x]=nem;
}
int getphi(int x){
    if(x<=1e7) return sumphi[x];
    if(mphi[x]) return mphi[x];
    int nem=x*(x+1)/2,r;
    for(int l=2;l<=x;l=r+1){
        r=x/(x/l);
        nem-=(r-l+1)*getphi(x/l);
    }
    return mphi[x]=nem;
}
signed main(){
    cin>>T;
    sieve();
    while(T--){
        cin>>n;
        cout<<getphi(n)<<" "<<getmu(n)<<"\n";
    }
    return 0;
}

P3455 [POI2007] ZAP-Queries

下面两道题的手推的式子

求：

经典把 \gcd(i,j)=k 转化为 \gcd(i,j)=1。

线性筛一遍莫比乌斯函数再用数论分块即可。

#include<bits/stdc++.h>
#define XD 114514
#define MAXN 50010
using namespace std;
int t,n,m,a,b,k,ans;
int mu[MAXN],num[MAXN];
bool vis[MAXN];
vector<int> prm;
void sieve(){
    mu[1]=1;num[1]=1;
    for(int i=2;i<=5e4;i++){
        if(!vis[i]){
            prm.push_back(i);
            mu[i]=-1;
        }
        num[i]=num[i-1]+mu[i];
        for(int j=0;j<prm.size() and i*prm[j]<=5e4;j++){
            vis[i*prm[j]]=true;
            if(i%prm[j]==0) break;
            mu[i*prm[j]]=-mu[i];
        }
    }
}
void calc(){
    int r;
    for(int l=1;l<=a;l=r+1){
        r=min(a/(a/l),b/(b/l));
        ans+=(num[r]-num[l-1])*(a/l)*(b/l);
    }
}
int main(){
    cin>>t;
    sieve();
    while(t--){
        ans=0;
        cin>>n>>m>>k;
        a=n/k;b=m/k;
        if(a>b) swap(a,b);
        calc();
        cout<<ans<<"\n";
    }
    return 0;
}

P2522 [HAOI2011] Problem b

就是用容斥原理搞一下，之后就和上面的一样了。

#include<bits/stdc++.h>
#define XD 114514
#define MAXN 50010
using namespace std;
int t,n,m,a,b,c,d,k;
int mu[MAXN],num[MAXN];
bool vis[MAXN];
vector<int> prm;
void sieve(){
    mu[1]=1;num[1]=1;
    for(int i=2;i<=5e4;i++){
        if(!vis[i]){
            prm.push_back(i);
            mu[i]=-1;
        }
        num[i]=num[i-1]+mu[i];
        for(int j=0;j<prm.size() and i*prm[j]<=5e4;j++){
            vis[i*prm[j]]=true;
            if(i%prm[j]==0) break;
            mu[i*prm[j]]=-mu[i];
        }
    }
}
int calc(int x,int y){
    n=x/k;m=y/k;
    if(n>m) swap(n,m);
    int r,ans=0;
    for(int l=1;l<=n;l=r+1){
        r=min(n/(n/l),m/(m/l));
        ans+=(num[r]-num[l-1])*(n/l)*(m/l);
    }
    return ans;
}
int main(){
    sieve();
    cin>>t;
    while(t--){
        cin>>a>>b>>c>>d>>k;
        cout<<calc(b,d)-calc(a-1,d)-calc(b,c-1)+calc(a-1,c-1)<<"\n";
    }
    return 0;
}

P2257 YY的GCD 和 PGCD - Primes in GCD Table

说句题外话，我是先写的这道题，再写的上面那两道题。QWQ

用图画手推的式子，很抽象，就当图个乐就行了。

下面正式推式子。

求：

经典把 \gcd(i,j)=k 转化为 \gcd(i,j)=1。

然后用莫比乌斯反演。

这里就把 n 当做 \min(n,m)。

我们设 T=pd，并将式子转换成枚举 T。

可以发现右面的式子可以提前求出。左面的式子数论分块即可。

#include<bits/stdc++.h>
#define XD 114514
#define MAXN 10000010
using namespace std;
int t,n,m;
int mu[MAXN],f[MAXN],sum[MAXN];
bool vis[MAXN];
vector<int> prm;
void sieve(){
    mu[1]=1;
    for(int i=2;i<=1e7;i++){
        if(!vis[i]){
            prm.push_back(i);
            mu[i]=-1;
        } 
        for(int j=0;j<prm.size();j++){
            if(i*prm[j]>1e7) break;
            vis[i*prm[j]]=true;
            if(i%prm[j]==0) break;
            mu[i*prm[j]]=-mu[i];
        }
    }
    for(int i=0;i<prm.size();i++){
        for(int j=1;j<=1e7;j++){
            if(prm[i]*j>1e7) break;
            f[prm[i]*j]+=mu[j];
        }
    }
    for(int i=1;i<=1e7;i++) sum[i]=sum[i-1]+f[i];
}
long long ans;
void solve(){
    int r=0;if(n>m) swap(n,m);
    for(int l=1;l<=n;l=r+1){
        int nem1=n/l,nem2=m/l;
        r=min(n/nem1,m/nem2);
        ans+=1ll*(sum[r]-sum[l-1])*(n/l)*(m/l);
    }
}
int main(){
    ios::sync_with_stdio(false);
    cin.tie(0);cout.tie(0);
    cin>>t;
    sieve();
    while(t--){
        ans=0;
        cin>>n>>m;
        solve();
        cout<<ans<<"\n";
    }
    return 0;
}

P3327 [SDOI2015] 约数个数和

传统异能

求：

首先，要知道 d 函数的一个性质：

所以原式就变成了：

设 p=\left\lfloor\frac{n}{d}\right\rfloor，q=\left\lfloor\frac{m}{d}\right\rfloor，则:

可以发现 \sum\limits_{x=1}^p{\left\lfloor\frac{p}{x}\right\rfloor} 和 \sum\limits_{y=1}^q \left\lfloor\frac{q}{y}\right\rfloor 可以预处理时存在数组里，我们设 f(x)=\sum\limits_{i=1}^x{\left\lfloor\frac{x}{i}\right\rfloor}。

于是查询时数论分块即可。（相当于数论分块套数论分块 QWQ）

#include<bits/stdc++.h>
#define XD 114514
#define MAXN 50010
#define int long long
using namespace std;
int T,n,m,ans;
bool vis[MAXN];
int mu[MAXN],sum[MAXN],num[MAXN];
vector<int> prm;
void sieve(){
    mu[1]=1;sum[1]=1;
    for(int i=2;i<=5e4;i++){
        if(!vis[i]){
            prm.push_back(i);
            mu[i]=-1;
        }
        for(int j=0;j<prm.size() and prm[j]*i<=5e4;j++){
            vis[i*prm[j]]=true;
            if(i%prm[j]==0) break;
            mu[i*prm[j]]=-mu[i];
        }
        sum[i]=sum[i-1]+mu[i];
    }
    for(int i=1;i<=5e4;i++){
        int r,nem=0;
        for(int l=1;l<=i;l=r+1){
            r=i/(i/l);
            nem+=(r-l+1)*(i/l);
        }
        num[i]=nem;
    }
}
void calc(){
    int r;
    for(int l=1;l<=n;l=r+1){
        r=min(n/(n/l),m/(m/l));
        ans+=(sum[r]-sum[l-1])*num[n/l]*num[m/l];
    }
}
signed main(){
    cin>>T;
    sieve();
    while(T--){
        ans=0;
        cin>>n>>m;
        if(n>m) swap(n,m);
        calc();
        cout<<ans<<"\n";
    }
    return 0;
}

P3768 简单的数学题

求：

为方便，设 m=\left\lfloor\frac{n}{d}\right\rfloor。