技术教程 编程教程 莫比乌斯反演/线性筛/积性函数/杜教筛/min25筛 学习笔记 正文

莫比乌斯反演/线性筛/积性函数/杜教筛/min25筛 学习笔记

2023-06-20,,

最近重新系统地学了下这几个知识点，以前没发现他们的联系，这次总结一下。

莫比乌斯反演入门：https://blog.csdn.net/litble/article/details/72804050

线性筛筛常见积性函数及其代码：https://blog.masterliu.net/algorithm/sieve/

积性函数与线性筛(包括普通线性函数)：https://blog.csdn.net/weixin_42562050/article/details/87997582

bzoj2154/bzoj2693/洛谷P1829 Crash的数字表格

这道题求  最终化到的式子是这样的 

那么主要矛盾就是怎么求后面的那坨，要深刻理解线性筛，然后分析i%pime[j]==0的时候prime[j]的新贡献和i%prime[j]!=0的时候prime[j]的贡献。

#include<bits/stdc++.h>

using namespace std;

typedef long long LL;

const int N=1e7+;

const int MOD=1e8+;

int n,m;

bool vis[N];

int tot=,pri[N],f[N],s[N];

void prework(int n) {

    f[]=;

    for (int i=;i<=n;i++) {

        if (!vis[i]) pri[++tot]=i,f[i]=(-i+MOD)%MOD;

        for (int j=;j<=tot&&i*pri[j]<=n;j++) {

            int k=i*pri[j]; vis[k]=;

            if (i%pri[j]) f[k]=(LL)f[i]*f[pri[j]]%MOD;

            else { f[k]=f[i]; break; }

        }

    }

    for (int i=;i<=n;i++) f[i]=(f[i-]+(LL)i*f[i]%MOD)%MOD;

    for (int i=;i<=n;i++) s[i]=((LL)i*(i+)/)%MOD;

}

int main()

{

    int T; cin>>T;

    prework();

    while (T--) {

        scanf("%d%d",&n,&m);

        if (n>m) swap(n,m);

        LL ans=;

        for (int i=,j;i<=n;i=j+) {

            j=min(n/(n/i),m/(m/i));

            LL tmp=(LL)s[n/i]*s[m/i]%MOD*(f[j]-f[i-]+MOD)%MOD;

            ans=(ans+tmp)%MOD;

        }

        printf("%lld\n",ans);

    }

    return ;

}

bzoj3994/洛谷P3327 约数个数和

要求 必须得先直到一个结论才能往下化简

化简结果是：  其中

那么我们发现g(x)其实就是d(x)(约数x的个数)的前缀和，g(n)=sigma(d(x)) (x=1~n) 。

那么线性筛之后求前缀和，直接分块即可。

#include<bits/stdc++.h>

using namespace std;

typedef long long LL;

const int N=5e4+;

int n,m;

bool vis[N];

LL tot=,pri[N],minp[N],d[N],u[N];

void prework(int n) {  //线性筛求莫比乌斯函数和约数个数

    d[]=; u[]=;

    for (int i=;i<=n;i++) {

        if (!vis[i]) pri[++tot]=i,minp[i]=,d[i]=,u[i]=-;

        for (int j=;j<=tot&&i*pri[j]<=n;j++) {

            int k=i*pri[j]; vis[k]=;

            if (i%pri[j]==) {

                u[k]=;

                minp[k]=minp[i]+;

                d[k]=d[i]/(minp[i]+)*(minp[i]+);

                break;

            }

            u[k]=-u[i];

            minp[k]=;

            d[k]=d[i]*;

        }

    }

    for (int i=;i<=n;i++) d[i]+=d[i-],u[i]+=u[i-];

}

int main()

{

    prework();

    int T; cin>>T;

    while (T--) {

        scanf("%d%d",&n,&m);

        if (n>m) swap(n,m);

        LL ans=;

        for (int i=,j;i<=n;i=j+) {

            j=min(n/(n/i),m/(m/i));

            ans+=(LL)d[n/i]*d[m/i]*(u[j]-u[i-]);

        }

        printf("%lld\n",ans);

    }

    return ;

}

bzoj3529/洛谷P3312 数表

要求 (其中σ(x)表示x的约数和) ，化简式子得到

注意到后面的一坨东西是迪利克雷卷积可以线性筛，但是这题有多组询问(n,m,a)，这个a的作用是限制σ(x)>a的不作贡献。这样的话就不能无脑筛前缀和然后直接分块计算了。

那么我们考虑离线，对询问的a从小到大排序。那么对于该次询问的a，我们逐一加入σ(x)<=a的σ(x)对答案的贡献，那么我们此时分块计算的就是正确的贡献了。那么就要需要一种能插入值又能区间和查询的数据结构了，BIT就是这样一个数据结构并且常数更小。每次加入一个σ(x)我们考虑对所有x的倍数更新σ(x)*μ(T/d)。  更新的时间复杂度是n/1+n/2+...n/n是无穷级数=nlogn不会超时。

#include<bits/stdc++.h>

using namespace std;

typedef long long LL;

const int N=1e5+;

int n,m,ans[N];

struct dat{

    int n,m,a,id;

    bool operator < (const dat &rhs) const {

        return a<rhs.a;

    }

}Q[N];

struct sumdat{

    int id,sum;

    bool operator < (const sumdat &rhs) const {

        return sum<rhs.sum;

    }

}d[N];

bool vis[N];

int tot=,pri[N],u[N],minp[N],sump[N];

void prework(int n) {

    u[]=; minp[]=; sump[]=;

    for (int i=;i<=n;i++) {

        if (!vis[i]) pri[++tot]=i,minp[i]=i,sump[i]=i+,u[i]=-;

        for (int j=;j<=tot&&i*pri[j]<=n;j++) {

            int k=i*pri[j]; vis[k]=;

            if (i%pri[j]==) {

                u[k]=;

                minp[k]=minp[i]*pri[j];

                LL tmp=(LL)sump[i]*((LL)minp[k]*pri[j]-);

                sump[k]=(int)(tmp/(minp[k]-));

                break;

            }

            u[k]=-u[i];

            minp[k]=pri[j];

            sump[k]=sump[i]*(+pri[j]);

        }

    }

    for (int i=;i<=n;i++) d[i].id=i,d[i].sum=sump[i];

}

int sum[N];

void update(int x,int v) {

    for (;x<=;x+=x&-x) sum[x]+=v;

}

int query(int x) {

    int ret=;

    for (;x;x-=x&-x) ret+=sum[x];

    return ret;

}

int solve(int n,int m) {

    if (n>m) swap(n,m);

    int ret=;

    for (int i=,j;i<=n;i=j+) {

        j=min(n/(n/i),m/(m/i));

        ret+=(n/i)*(m/i)*(query(j)-query(i-));

    }

    return ret;

}

int main()

{

    prework();

    int T; cin>>T;

    for (int i=;i<=T;i++) {

        scanf("%d%d%d",&Q[i].n,&Q[i].m,&Q[i].a);

        Q[i].id=i;

    }

    sort(Q+,Q+T+);

    sort(d+,d++);

    int now=;

    for (int i=;i<=T;i++) {

        while (now<= && d[now].sum<=Q[i].a) {

            for (int i=d[now].id;i<=;i+=d[now].id)

                update(i,d[now].sum*u[i/d[now].id]);

            now++;

        }

        ans[Q[i].id]=solve(Q[i].n,Q[i].m);

    }

    for (int i=;i<=T;i++) {

        if (ans[i]<) ans[i]+=,ans[i]++;

        printf("%d\n",ans[i]);

    }

    return ;

}

洛谷P3768 简单的数学题

要求     常规套路推式子得到

后面的可以分块得到，那么主要矛盾变成怎么求T^2*Φ(T)的前缀和。这里就要用到杜教筛了。

令g(i)=i^1，那么他两的迪利克雷卷积就是

带入杜教筛的那条式子得到 积累两条公式 ①sigma(i^3)(i=1~n) = ( sigma(i)(i=1~n) )^2  ②sigma(i^2) = (i)(i+1)(2i+1)/6 。

那么预处理前2e6项，后面的记忆化搜索就可以了。

#include<bits/stdc++.h>

using namespace std;

typedef long long LL;

const int N=2e6+;

const int Pr=2e6;

int MOD; LL inv2,inv6;

map<LL,LL> mp;

LL power(LL x,LL p) {

    LL ret=;

    for (;p;p>>=) {

        if (p&) ret=ret*x%MOD;

        x=x*x%MOD;

    }

    return ret;

}

bool vis[N];

int tot=,pri[N]; LL mu[N],phi[N];

void prework() {

    vis[]=; mu[]=phi[]=;

    for (int i=;i<=Pr;i++) {

        if (!vis[i]) pri[++tot]=i,mu[i]=-,phi[i]=i-;

        for (int j=;j<=tot&&i*pri[j]<=Pr;j++) {

            int k=i*pri[j]; vis[k]=;

            if (i%pri[j]==) {

                mu[k]=; phi[k]=phi[i]*pri[j]%MOD;

                break;

            }

            mu[k]=-mu[i]; phi[k]=phi[i]*phi[pri[j]]%MOD;

        }

    }

    for (int i=;i<=Pr;i++) phi[i]=phi[i]*i%MOD*i%MOD;

    for (int i=;i<=Pr;i++) phi[i]+=phi[i-],phi[i]%=MOD;

}

LL sq2(LL n) { return n%MOD*((n+)%MOD)%MOD*((*n+)%MOD)%MOD*inv6%MOD; }

LL sq3(LL n) {

    LL ret=(n%MOD)*((n+)%MOD)%MOD*inv2%MOD;

    return ret*ret%MOD;

}

LL solve(LL n) {

    if (n<=Pr) return phi[n];

    if (mp.count(n)) return mp[n];

    LL ret=,i=,j;

    while (i<=n) {

        j=n/(n/i);

        ret+=solve(n/i)*((sq2(j)-sq2(i-)+MOD)%MOD)%MOD;

        ret%=MOD;

        i=j+;

    }

    ret=(sq3(n)-ret)%MOD; ret=(ret+MOD)%MOD;

    return mp[n]=ret;

}

int main()

{

    LL n; cin>>MOD>>n;

    prework();

    inv2=power(,MOD-);

    inv6=power(,MOD-);

    LL i=,j,ans=;

    while (i<=n) {

        j=n/(n/i);

        LL sum=(n/i)%MOD*((n/i+)%MOD)%MOD*inv2%MOD;

        ans+=(solve(j)-solve(i-)+MOD)%MOD*sum%MOD*sum%MOD; ans%=MOD;

        i=j+;

    }

    cout<<ans<<endl;

    return ;

}

洛谷P3704 数字表格

这道题要是蒟蒻自己想想破头也想不出来，题解参考洛谷题解。

要求  枚举gcd(i,j)变成  常规套路化简得到

接下来是关键一步：把上面的sigma拉下来    （原理相当于指数的相加拉下来变成相乘 a^(x+y)=a^x * a^y）

继续莫比乌斯反演套路，改为枚举T=td 得到 发现后面的[N/T] [M/T] 可以用分块，那么问题变成怎么快速计算中间的一坨

想了想感觉还是不能线性筛，但是一看数据范围1e6，我们可以暴力筛时间复杂度是 n/1+n/2+n3/+...n/n=nlogn 可以通过。

那么就是暴力预处理然后对于每个询问分块回答即可。

#include<bits/stdc++.h>

using namespace std;

typedef long long LL;

const int MOD=1e9+;

const int N=1e6+;

LL power(LL x,LL p) {

    LL ret=;

    for (;p;p>>=) {

        if (p&) ret=ret*x%MOD;

        x=x*x%MOD;

    }

    return ret;

}

bool vis[N]; LL tot=,pri[N],mu[N],f[N],inv[N],s[N];

void prework(int n) {

    vis[]=; mu[]=;

    for (int i=;i<=n;i++) {

        if (!vis[i]) pri[++tot]=i,mu[i]=-;

        for (int j=;j<=tot&&i*pri[j]<=n;j++) {

            int k=i*pri[j]; vis[k]=;

            if (i%pri[j]==) {

                mu[k]=; break;

            }

            mu[k]=-mu[i];

        }

    }

    f[]=; f[]=; for (int i=;i<=n;i++) f[i]=(f[i-]+f[i-])%MOD;

    for (int i=;i<=n;i++) s[i]=;

    for (int i=;i<=n;i++) inv[i]=power(f[i],MOD-);

    for (int i=;i<=n;i++)

        for (int j=i;j<=n;j+=i)

            if (mu[j/i]==) s[j]=s[j]*f[i]%MOD;

            else if (mu[j/i]==-) s[j]=s[j]*inv[i]%MOD;

    for (int i=;i<=n;i++) s[i]=(s[i-]*s[i])%MOD;

}

int main()

{

    int T; cin>>T;

    prework();

    while (T--) {

        LL n,m; scanf("%lld%lld",&n,&m);

        if (n>m) swap(n,m);

        LL i=,j,tmp,ans=;

        while (i<=n) {

            j=min(n/(n/i),m/(m/i));

            tmp=s[j]*power(s[i-],MOD-)%MOD;

            ans*=power(tmp,(n/i)*(m/i)%(MOD-)); ans%=MOD;

            i=j+;

        }

        printf("%lld\n",ans);

    }

    return ;

}

杜教筛

比起Min25筛杜教筛还是比较好学的。

当成黑盒算法的话其实就一条式子：

s(n)=sigma(u[i]) (i=1~n) 就是我们要求的积性函数u的前缀和，然后我们找到一个g函数使得g函数和所求u函数的迪利克雷卷积非常好算（这样的话式子的前半部分就好算，后半部分可以通过分块计算），那么我们就可以通过上面这条式子快速计算s(n)。

给出求莫比乌斯函数前缀和和欧拉函数前缀和的模板。(因为使用map，交到OJ上会TLE)

#include<bits/stdc++.h>

using namespace std;

typedef long long LL;

const int N=3e6+;

const int Pr=3e6;

unordered_map<LL,LL> mm,mp;

bool vis[N];

int tot=,pri[N]; LL mu[N],phi[N];

void prework() {

    vis[]=; mu[]=phi[]=;

    for (int i=;i<=Pr;i++) {

        if (!vis[i]) pri[++tot]=i,mu[i]=-,phi[i]=i-;

        for (int j=;j<=tot&&i*pri[j]<=Pr;j++) {

            int k=i*pri[j]; vis[k]=;

            if (i%pri[j]==) {

                mu[k]=; phi[k]=phi[i]*pri[j];

                break;

            }

            mu[k]=-mu[i]; phi[k]=phi[i]*phi[pri[j]];

        }

    }

    for (int i=;i<=Pr;i++) mu[i]+=mu[i-];

    for (int i=;i<=Pr;i++) phi[i]+=phi[i-];

}

LL getmu(LL n) {

    if (n<=Pr) return mu[n];

    if (mm.count(n)) return mm[n];

    LL ret=,i=,j;

    while (i<=n) {

        j=n/(n/i);

        ret+=(j-i+)*getmu(n/i);

        i=j+;

    }

    return mm[n]=-ret;

}

LL getphi(LL n) {

    if (n<=Pr) return phi[n];  //预处理的

    if (mp.count(n)) return mp[n];  //算过的

    LL ret=,i=,j;

    while (i<=n) {  //数论分块求出后面一坨

        j=n/(n/i);

        ret+=(j-i+)*getphi(n/i);

        i=j+;

    }

    return mp[n]=n*(n+)/-ret;

}

int main()

{

    prework();

    int T; cin>>T;

    while (T--) {

        LL n; scanf("%lld",&n);

        printf("%lld %lld\n",getphi(n),getmu(n));

    }

    return ;

}

BZOJ-4916

主要是求phi(i^2)的前缀和，f(i) = phi(i^2)=i*phi(i) 。像办法构造一个g(i)和f(i)做迪利克雷卷积结果方便求。

那么sigma f(i)*g(n/i)=sigma i*phi(i)*g(n/i)，我们也知道sigma( phi(i) ) (i|n) = n 那么g是什么会使式子变得简单呢？纯粹的想法就是把式子的这个i消掉，那么g(i)=i就行了。

#include<bits/stdc++.h>

using namespace std;

typedef long long LL;

const int N=1e7+;

const int Pr=1e7;

const int MOD=1e9+;

const int inv=;  //6的逆元

map<LL,LL> mp;

bool vis[N];

int tot=,pri[N]; LL mu[N],phi[N];

void prework() {

    vis[]=; mu[]=phi[]=;

    for (int i=;i<=Pr;i++) {

        if (!vis[i]) pri[++tot]=i,mu[i]=-,phi[i]=i-;

        for (int j=;j<=tot&&i*pri[j]<=Pr;j++) {

            int k=i*pri[j]; vis[k]=;

            if (i%pri[j]==) {

                mu[k]=; phi[k]=phi[i]*pri[j]%MOD;

                break;

            }

            mu[k]=-mu[i]; phi[k]=phi[i]*phi[pri[j]]%MOD;

        }

    }

    for (int i=;i<=Pr;i++) phi[i]=phi[i]*i%MOD;

    for (int i=;i<=Pr;i++) phi[i]+=phi[i-],phi[i]%=MOD;

}

LL solve(LL n) {

    if (n<=Pr) return phi[n];

    if (mp.count(n)) return mp[n];

    LL ret=,i=,j;

    while (i<=n) {

        j=n/(n/i);

        ret+=(i+j)*(j-i+)/%MOD*solve(n/i)%MOD;

        ret%=MOD;

        i=j+;

    }

    ret=n*(n+)%MOD*(*n+)%MOD*inv%MOD-ret;

    ret=(ret%MOD+MOD)%MOD;

    return mp[n]=ret;

}

int main()

{

    prework();

    LL n; cin>>n;

    printf("%lld\n%lld\n",,solve(n));

    return ;

}

洛谷P3172 选数

这道题用到的莫比乌斯反演部分还是比较简单的。但是注意到n<=1e9，求莫比乌斯函数就不能用线性筛要杜教筛。

求  莫比乌斯套路一下   （这个式子有一点小错误，后面应该是除以x而不是d）

这道题还有一些小细节：①注意因为要求gcd为k，所以一开始就令R/=k,--L/=k后面会好处理很多。②因为这里数论分块要求分到R而不是较小的L，所以当L/i==0的时候使其为INF。

#include<bits/stdc++.h>

using namespace std;

typedef long long LL;

const int P=1e9+;

const int Pr=2e6;

const int N=2e6+;

const int INF=0x3f3f3f3f;

int l,r,n,k;

unordered_map<LL,LL> mm;

bool vis[N];

int tot=,pri[N]; LL mu[N];

void prework() {

    vis[]=; mu[]=;

    for (int i=;i<=Pr;i++) {

        if (!vis[i]) pri[++tot]=i,mu[i]=-;

        for (int j=;j<=tot&&i*pri[j]<=Pr;j++) {

            int k=i*pri[j]; vis[k]=;

            if (i%pri[j]==) {

                mu[k]=; break;

            }

            mu[k]=-mu[i];

        }

    }

    for (int i=;i<=Pr;i++) mu[i]+=mu[i-];

}

LL getmu(LL n) {

    if (n<=Pr) return mu[n];

    if (mm.count(n)) return mm[n];

    LL ret=,i=,j;

    while (i<=n) {

        j=n/(n/i);

        ret+=(j-i+)*getmu(n/i)%P;

        i=j+;

    }

    ret=-ret; ret=(ret%P+P)%P;

    return mm[n]=ret;

}

LL power(LL x,LL p) {

    LL ret=;

    for (;p;p>>=) {

        if (p&) ret=ret*x%P;

        x=x*x%P;

    }

    return ret;

}

int main()

{

    cin>>n>>k>>l>>r;

    r/=k; l=(l-)/k;

    prework();

    LL i=,j,ans=;

    while (i<=r) {

        j=min(l/i?l/(l/i):INF,r/(r/i));

        ans+=(getmu(j)-getmu(i-))*power(r/i-l/i,n)%P;

        ans=(ans%P+P)%P;

        i=j+;

    }

    cout<<ans<<endl;

    return ;

}

51Nod平均最小公倍数

莫比乌斯反演/线性筛/积性函数/杜教筛/min25筛 学习笔记的相关教程结束。

《莫比乌斯反演/线性筛/积性函数/杜教筛/min25筛 学习笔记.doc》

下载本文的Word格式文档，以方便收藏与打印。

Copyright 2023. 昆居客 - IT编程,互联网信息技术文档大全!