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

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

莫比乌斯反演入门：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平均最小公倍数