【学习笔记】狄利克雷卷积与高级筛法

合集 - 学习笔记(42)

1.【学习笔记】Kruskal 重构树2022-11-252.【学习笔记】网络流2022-12-153.【学习笔记】高级数据结构2022-12-214.【学习笔记】线性基2022-12-275.【学习笔记】Link Cut Tree2022-12-286.【学习笔记】字符串后缀算法2022-12-307.【学习笔记】字符串回文算法2023-01-128.【学习笔记】组合数学2023-01-279.【学习笔记】多项式 1：基础操作2023-01-3010.【学习笔记】多项式 2：集合幂级数2023-02-0311.【学习笔记】多项式 3：多项式运算2023-02-1112.【学习笔记】Prufer 序列2023-02-1513.【学习笔记】多项式 4：生成函数2023-02-1514.【学习笔记】DP 套 DP2023-02-1915.【学习笔记】图的连通性2023-02-2216.【学习笔记】差分约束2023-03-2017.【学习笔记】长链剖分2023-04-2218.【学习笔记】2-SAT2023-04-2319.【学习笔记】根号算法2023-06-0620.【学习笔记】Primal-Dual 原始对偶算法2023-06-1521.【学习笔记】Bostan-Mori 算法2023-07-02

22.【学习笔记】狄利克雷卷积与高级筛法2023-07-03

23.【学习笔记】DP 优化 1：基础优化2023-07-0324.【学习笔记】DP 优化 2：动态 DP2023-07-0425.【学习笔记】李超线段树2023-07-0926.【学习笔记】优化建图2023-07-1127.【学习笔记】Segment Tree Beats2023-07-1128.【学习笔记】插头 DP2023-07-1329.【学习笔记】任意模数多项式乘法2023-07-1630.【学习笔记】SG 函数与 SG 定理2023-08-0531.【学习笔记】类欧几里得算法2023-08-0632.【学习笔记】狄利克雷前/后缀和/差分2023-08-1133.【学习笔记】DSU on Tree2023-08-2234.【学习笔记】DP 优化 3：闵可夫斯基和优化 DP2023-09-0135.【学习笔记】笛卡尔树2023-09-0536.【学习笔记】Miller-Rabin 算法2023-10-2437.【学习笔记】DP 优化 4：决策单调性2023-11-1538.【学习笔记】DP 优化 5：wqs 二分优化 DP2023-12-0439.【学习笔记】边分治2023-12-2840.【学习笔记】KMP 相关算法2024-01-1041.【学习笔记】概率生成函数2024-01-1142.【学习笔记】离散对数和剩余2024-01-25

收起

狄利克雷卷积#

概念#

对于数论函数 $f,g$，定义其狄利克雷卷积 $h=f*g$，满足：

$$h(n)=(f*g)(n)=\sum_{d\mid n} f(d)g\left(\dfrac{n}{d}\right)$$

运算律：

• 满足交换律，显然具有对称性。

• 满足结合律，等价于三个 $d_i$ 贡献到 $n$。

• 满足加法的分配率。

常见数论函数：

• $\mathrm{I}$，常函数，值恒为 $1$。

• $\mathrm{Id}_k$，幂函数，值为 $n^k$。

• $\epsilon$，狄利克雷卷积的单位元，除 $\epsilon(1)=1$ 外其余均为 $0$。

• $\mu$，莫比乌斯函数，狄利克雷卷积的逆元。

• $\varphi$，欧拉函数。

• $\sigma_k$，除数函数，定义 $\sigma_k(n)=\sum_{d\mid n}d^k$。

常见数论函数卷积变换#

• $\mu *\mathrm{I}=\epsilon$

• $\varphi *\mathrm{I}=\mathrm{Id}$

• $\mathrm{Id}_k *\mathrm{I}=\sigma_k$

• $\mu* \mathrm{Id}=\varphi$

• $(\mu \cdot \mathrm{Id}_k)*\mathrm{Id}_k=\epsilon$

• $(\varphi\cdot \mathrm{Id}_k)*\mathrm{Id}_k=\mathrm{Id}_{k+1}$

杜教筛#

原理#

求数论函数前缀和。

设 $h=f*g,S=\sum f$，容易推得：

$$\begin{aligned} \sum_{i=1}^n h(i)&=\sum_{i=1}^n \sum_{j\mid i} f(j)g\left(\dfrac{i}{j}\right)\\ &=\sum_{j=1}^n g(j)\sum_{i=1}^{\left\lfloor n/i\right\rfloor} f(i)\\ &=\sum_{i=1}^n g(i)S_f\left(\left\lfloor \dfrac{n}{i}\right\rfloor\right) \end{aligned}$$

整理一下可以写成：

$$S_f(n)=\sum_{i=1}^n h(i)-\sum_{i=2}^{n} g(i)S_f\left(\left\lfloor \dfrac{n}{i}\right\rfloor\right)$$

如果构造出可以快速计算前缀和的函数 $g,h$，就可以递归求解了。

时间复杂度#

首先需要证明的是：

$$\left\lfloor\dfrac{\left\lfloor\tfrac{n}{a}\right\rfloor}{b}\right\rfloor=\left\lfloor\dfrac{n}{ab}\right\rfloor$$

设 $\mathrm{RHS}=k$，则有：

$$k\times ab\le n<(k+1)\times ab$$

于是：

$$k\times b\le \dfrac{n}{a}<(k+1)\times b$$

从而：

$$k\times b\le \left\lfloor\dfrac{n}{a}\right\rfloor<(k+1)\times b$$

$$k\le \dfrac{\left\lfloor\tfrac{n}{a}\right\rfloor}{b}<k+1$$

即：

$$\left\lfloor\dfrac{\left\lfloor\tfrac{n}{a}\right\rfloor}{b}\right\rfloor=\left\lfloor\dfrac{n}{ab}\right\rfloor$$

---

这说明当我们数论分块向下递归时，递归到 $n'=\left\lfloor \dfrac{n}{d}\right\rfloor$，枚举 $d'$ 得到 $\left\lfloor \dfrac{n'}{d'}\right\rfloor=\left\lfloor \dfrac{\left\lfloor \tfrac{n}{d}\right\rfloor}{d'}\right\rfloor$，一定可以表示为 $\left\lfloor \dfrac{n}{k}\right\rfloor$ 的形式，而这一定在第一层递归时处理到（不考虑先后顺序），于是在记忆化的加持下，计算复杂度只需要考虑第一层递归即可。

设当前块长 $i$，对复杂度的贡献就是 $O\left(\sqrt{\dfrac{n}{i}}\right)$。

当 $i> \sqrt{n}$ 时，第二层最大 $O(n^{1/4})$，所以总复杂度 $O(n^{3/4})$。

当 $i\le \sqrt{n}$ 时，就写成：

$$\begin{aligned} T(n)&=O\left(\sum_{i=1}^{\sqrt{n}} \sqrt{\dfrac{n}{i}}\right)\\ &=O\left(\sqrt{n}\int_{1}^n \dfrac{1}{\sqrt{x}} \mathrm{d}x\right)\\ &=O(n^{3/4}) \end{aligned}$$

但较小的前缀和是可以预处理出的，假设预处理出了 $[1,n^k]$ 的前缀和，那么也就是只有 $\dfrac{n}{i}>n^k$ 的部分才需要递归，也就变成：

$$\begin{aligned} T(n)&=O\left(n^k+\sum_{i=1}^{n^{1-k}} \sqrt{\dfrac{n}{i}}\right)\\ &=O\left(n^k+\sqrt{n}\int_{1}^{n^{1-k}} \dfrac{1}{\sqrt{x}}\mathrm{d}x\right)\\ &=O(n^k+n^{1-k/2}) \end{aligned}$$

取 $k=\dfrac{2}{3}$ 得到理论最优复杂度 $O(n^{2/3})$。

数论分块时（无论嵌套几层）使用杜教筛求数论函数前缀和，由于所有询问的区间右端点为 $\left\lfloor \dfrac{n}{\left\lfloor \tfrac{n}{d}\right\rfloor}\right\rfloor$，可以表示为上述记忆化的形式，所以并不会增加复杂度。

点击查看代码

namespace Phi{
    int pr[maxn],mn[maxn],mnpw[maxn],phi[maxn];
    bool vis[maxn];
    ull s1[maxn],s2[maxn];
    bool mark[maxn];
    inline void linear_sieve(){
        mn[1]=1,phi[1]=1;
        for(int i=2;i<=lim;++i){
            if(!vis[i]) pr[++pr[0]]=i,mn[i]=i,mnpw[i]=i,phi[i]=i-1;
            for(int j=1;j<=pr[0]&&i*pr[j]<=lim;++j){
                vis[i*pr[j]]=1,mn[i*pr[j]]=pr[j];
                if(mn[i]==pr[j]) mnpw[i*pr[j]]=mnpw[i]*pr[j];
                else mnpw[i*pr[j]]=pr[j];
                if(i*pr[j]==mnpw[i*pr[j]]) phi[i*pr[j]]=phi[i]*pr[j];
                else phi[i*pr[j]]=phi[i*pr[j]/mnpw[i*pr[j]]]*phi[mnpw[i*pr[j]]];
                if(i%pr[j]==0) break;
            }
        }
        for(int i=1;i<=lim;++i) s1[i]=s1[i-1]+phi[i];
    }
    ull get_sum(int N){
        if(N<=lim) return s1[N];
        if(mark[n/N]) return s2[n/N];
        mark[n/N]=1;
        ull sum=0;
        for(unsigned int l=2,r;l<=N;l=r+1){
            r=N/(N/l);
            sum+=1ull*(r-l+1)*get_sum(N/l);
        }
        if(N&1) sum=1ull*(N+1)/2*N-sum;
        else sum=1ull*N/2*(N+1)-sum;
        return s2[n/N]=sum;
    }
    inline ull solve(){
        memset(mark,0,sizeof(mark));
        return get_sum(n);
    }
}

namespace Mu{
    int pr[maxn],mu[maxn];
    bool vis[maxn];
    int s1[maxn],s2[maxn];
    bool mark[maxn];
    inline void linear_sieve(){
        mu[1]=1;
        for(int i=2;i<=lim;++i){
            if(!vis[i]) pr[++pr[0]]=i,mu[i]=-1;
            for(int j=1;j<=pr[0]&&i*pr[j]<=lim;++j){
                vis[i*pr[j]]=1,mu[i*pr[j]]=-mu[i];
                if(i%pr[j]==0){
                    mu[i*pr[j]]=0;
                    break;
                }
            }
        }
        for(int i=1;i<=lim;++i) s1[i]=s1[i-1]+mu[i];
    }
    int get_sum(int N){
        if(N<=lim) return s1[N];
        if(mark[n/N]) return s2[n/N];
        mark[n/N]=1;
        ull sum=0;
        for(unsigned int l=2,r;l<=N;l=r+1){
            r=N/(N/l);
            sum+=(r-l+1)*get_sum(N/l);
        }
        sum=1-sum;
        return s2[n/N]=sum;
    }
    inline int solve(){
        memset(mark,0,sizeof(mark));
        return get_sum(n);
    }
}

常用构造#

最常用的是 $\mu *\mathrm{I}=\epsilon$ 以及 $\varphi *\mathrm{I}=\mathrm{Id}$。

在此基础上有：$(\mu\cdot \mathrm{Id}_k)*\mathrm{Id}_k=\epsilon$，$(\varphi\cdot \mathrm{Id}_k) * \mathrm{Id}_{k}=\mathrm{Id}_{k+1}$。

另外根据 $\sigma_k$ 的定义有 $\mathrm{Id}_k*\mathrm{I}=\sigma_k$。

例题#

Luogu-P5218 无聊的水题 II#

容易发现是要求选出一个集合满足 $\gcd=1$。

设 $F(n)$ 为值域 $[1,n]$ 的方案数，那么 $\gcd=k$ 就是 $F\left(\left\lfloor \dfrac{n}{k}\right\rfloor\right)$，可以单步容斥。

$$F(n)=2^n-1-\sum_{k=2}^n F\left(\left\lfloor \dfrac{n}{k}\right\rfloor\right)$$

式子长得像杜教筛，复杂度 $O(n^{3/4})$ 不太好过。

考虑设 $f(d)$ 为 $\gcd=d$ 的方案数，$g(d)$ 为 $d \mid \gcd$ 的方案数，则 $g(d)=2^{\left\lfloor n/d\right\rfloor}-1$，且有：

$$g(d)=\sum_{d\mid n} f(n)$$

根据莫比乌斯反演有：

$$f(d)=\sum_{d\mid n} \mu\left(\dfrac{n}{d}\right) g(n)$$

所以：

$$f(1)=\sum_{i=1}^n \mu(i)g(i)=\sum_{i=1}^n \mu(i)(2^{\left\lfloor n/i\right\rfloor}-1)$$

杜教筛预处理+数论分块即可，使用光速幂复杂度 $O(n^{2/3})$。

LibreOJ-6229 这是一道简单的数学题#

把 $j$ 的上下界改成 $[1,n]$，处理一下 $i=j$ 的情况即可。

简单反演得到：

$$\sum_{T=1}^n\sum_{d\mid T} \mu(d)d^2 S_1\left(\left\lfloor \dfrac{n}{d}\right\rfloor\right)^2$$

$S_1(x)$ 是 $1$ 次自然数幂和。

于是要求 $\sum_{d\mid n} \mu(d)d^2$ 的前缀和，不难发现这个函数是 $(\mu \cdot \mathrm{Id}_2)* \mathrm{I}$，根据狄利克雷卷积的交换律以及结合律，可以得到 $(\mu \cdot \mathrm{Id}_2)*\mathrm{I}*\mathrm{Id}_2=(\mu \cdot \mathrm{Id}_2)*\mathrm{Id}_2*\mathrm{I}=\epsilon*\mathrm{I}=\mathrm{I}$。

于是有：

$$S_f(n)=n-\sum_{i=2}^{n} i^2S_f\left(\left\lfloor \dfrac{n}{i}\right\rfloor\right)$$

杜教筛处理即可。

LibreOJ-6491 XXOI2018 简单的最大公约数#

和上面的题差不多，先设 $f(m)$ 表示值域 $[1,m]$，$n$ 个数 $\gcd=1$ 的方案数，答案是 $\sum_{k=1}^m k\times f\left(\left\lfloor \dfrac{m}{k}\right\rfloor\right)$，最终求解用数论分块。

$$f(m)=m^n-\sum_{i=2}^m f\left(\left\lfloor \dfrac{m}{i}\right\rfloor\right)$$

直接算还是 $O(n^{3/4})$。

希望预处理一些 $f$。

依旧是莫比乌斯反演，把 $f(m,k)$ 扩充到值域 $[1,m]$，$\gcd=k$ 的方案数，$g(m,k)$ 定义为值域 $[1,m]$，$k\mid \gcd$ 的方案数，于是反演得到：

$$g(m,k)=\left\lfloor \dfrac{m}{k}\right\rfloor^{n}$$

$$f(m,k)=\sum_{k\mid m'} \mu\left(\dfrac{m'}{k}\right) \left\lfloor \dfrac{m}{m'}\right\rfloor^{n}$$

只关心 $k=1$，所以有：

$$f(m)=\sum_{k=1}^m \mu(k) \left\lfloor \dfrac{m}{k}\right\rfloor^{n}$$

这个直接求一行比较困难，考虑差分：

$$\begin{aligned} \Delta f(m)&=f(m)-f(m-1)\\ &=\sum_{k=1}^m \mu(k) \left(\left\lfloor \dfrac{m}{k}\right\rfloor^{n}-\left\lfloor \dfrac{m-1}{k}\right\rfloor^{n}\right) \end{aligned}$$

分析一下什么时候 $\left\lfloor \dfrac{m}{k}\right\rfloor\neq \left\lfloor \dfrac{m-1}{k}\right\rfloor$，显然一定是 $k\mid m$，即 $m$ 恰好到了下一个块里。

于是这个预处理就是调和级数 $O(n\log n)$ 的，由于 $\mu$ 的性质，实际产生贡献的会更少。

这样大致处理 $2\times 10^6$ 左右，再套类似杜教筛的东西就可以通过了。

---

另一个做法是使用欧拉反演，即利用 $\varphi *\mathrm{I}=\mathrm{Id}$ 的性质，可以得到：

$$\begin{aligned} &\sum_{i_1=1}^m\sum_{i_2=1}^m\cdots\sum_{i_n=1}^m \gcd(i_1,i_2,\cdots,i_n)\\ =&\sum_{i_1=1}^m\sum_{i_2=1}^m\cdots\sum_{i_n=1}^m \sum_{d\mid \gcd(i_1,i_2,\cdots,i_n)} \varphi(d)\\ =&\sum_{d=1}^m \varphi(d) \left\lfloor\dfrac{m}{d}\right\rfloor^n \end{aligned}$$

杜教筛朴素计算即可。

PN 筛#

Powerful Number#

定义一个数 $n$ 是 Powerful Number(PN) 当且仅当 $n=\prod_{p\in \mathbf{P}} p_i^{c_i},c_i>1$。

每个 PN 一定可以被唯一表示为 $a^2b^3$ 的形式，次数为偶数的放在 $a$ 部分，次数为奇数的将 $3$ 次放在 $b$ 部分，剩余放在 $a$ 部分。

于是可以通过枚举 $a$ 的值来计算 PN 的个数，大致是：

$$O\left(\sum_{i=1}^{\sqrt{n}}\sqrt[3]{\dfrac{n}{i^2}}\right)=O\left(\int_1^{\sqrt{n}} \sqrt[3]{\dfrac{n}{i^2}}\right)=O(\sqrt{n})$$

所以可以通过 DFS 来得出所有 PN。

求积性函数前缀和#

以 模板题 为例。

给定积性函数 $f$，求其前缀和。

考虑构造积性函数 $g$，满足 $g(p)=f(p)$，这里 $g(p)=f(p)=p(p-1)$，可以构造 $g=\varphi\cdot\mathrm{Id}$。

再找到积性函数 $h$，满足 $f=g*h$，由于 $f(p)=g(1)h(p)+g(p)h(1)$，积性函数 $1$ 处点值一定为 $1$，所以 $f(p)=h(p)+g(p)$，即 $h(p)=0$，这样所有可以被质因子筛到的点值处都是 $0$，也就是所有非 PN 点值都是 $0$。

这样化简 $f$ 前缀和的式子：

$$\begin{aligned} &\sum_{i=1}^n f(i)\\ =&\sum_{i=1}^n \sum_{j\mid i}g(i)h\left(\dfrac{i}{j}\right)\\ =&\sum_{i=1}^n h(i)S_g\left(\left\lfloor\dfrac{n}{i}\right\rfloor\right)\\ =&\sum_{\substack{i=1\\ i\text{ is PN}}}^n h(i)S_g\left(\left\lfloor\dfrac{n}{i}\right\rfloor\right) \end{aligned}$$

$g$ 的前缀和可以用杜教筛求出（有更优秀的也可以），计算答案在 DFS 求出所有 PN 时进行，其中求 $h$ 可以求出所有 $h(p^c)$ 再累乘。

求 $h(p^c)$ 最常规的办法是利用 $f=g*h$，移项得到 $h(p^c)=f(p^c)-\sum_{i=0}^{c-1} h(p^i)g(p^{c-i})$。

点击查看代码

ll n;
int pr[maptrn],phi[maptrn];
bool vis[maptrn];
int s1[maptrn],s2[maptrn];
bool mark[maptrn];
int pw[maptrn/10][40],h[maptrn/10][40];
inline void linear_sieve(){
    phi[1]=1;
    for(int i=2;i<=lim;++i){
        if(!vis[i]) pr[++pr[0]]=i,phi[i]=i-1;
        for(int j=1;j<=pr[0]&&i*pr[j]<=lim;++j){
            vis[i*pr[j]]=1,phi[i*pr[j]]=phi[i]*phi[pr[j]];
            if(i%pr[j]==0){
                phi[i*pr[j]]=phi[i]*pr[j];
                break;
            }        
        }
    }
    for(int i=1;i<=lim;++i) s1[i]=(s1[i-1]+1ll*i*phi[i]%mod)%mod;
    for(int i=1;i<=pr[0];++i){
        pw[i][0]=1,h[i][0]=1;
        for(ll j=pr[i],e=1;;j*=pr[i],++e){
            pw[i][e]=j%mod;
            int now=1ll*pw[i][e]*(pw[i][e]-1)%mod;
            for(int k=0;k<e;++k){
                now=(now-1ll*h[i][k]*(pr[i]-1)%mod*pw[i][e-k-1]%mod*pw[i][e-k]%mod+mod)%mod;
            }
            h[i][e]=now;
            if(j>n/pr[i]) break;
        }
    }
}
int get_sumg(ll N){
    if(N<=lim) return s1[N];
    if(mark[n/N]) return s2[n/N]; 
    mark[n/N]=1;
    int tmpn=N%mod;
    int res=1ll*tmpn*(tmpn+1)%mod*(2ll*tmpn%mod+1)%mod*inv6%mod;
    for(ll l=2,r;l<=N;l=r+1){
        r=N/(N/l);
        int tmpl=(l-1)%mod,tmpr=r%mod;
        int now=(1ll*tmpr*(tmpr+1)/2-1ll*tmpl*(tmpl+1)/2)%mod;
        res=(res-1ll*now*get_sumg(N/l)%mod+mod)%mod;
    }
    return s2[n/N]=res;
}
int ans;
void dfs(ll N,int ptr,int res){
    ans=(ans+1ll*res*get_sumg(n/N)%mod)%mod;
    if(ptr>pr[0]) return;
    for(int i=ptr;i<=pr[0];++i){
        if(N>n/(1ll*pr[i]*pr[i])) break;
        for(ll j=N*pr[i]*pr[i],e=2;;j*=pr[i],++e){
            dfs(j,i+1,1ll*res*h[i][e]%mod);
            if(j>n/pr[i]) break;
        }
    }
}

int main(){
    scanf("%lld",&n);
    linear_sieve();
    dfs(1,1,1);
    printf("%d\n",ans);
    return 0;
}

复杂度分析#

求 $g$ 前缀和复杂度认为是杜教筛复杂度 $O(n^{2/3})$。

只预处理 $[1,\sqrt{n}]$ 的质数，每个质数枚举 $p^c$ 并计算 $h$ 的复杂度是 $O(\log^2 n)$，素数密度 $\pi(n)=O\left(\dfrac{n}{\log n}\right)$，所以总复杂度是 $O(\sqrt{n}\log n)$，实际这个上界比较松。

于是得到大致是 $O(n^{2/3})$ 的复杂度。

算法流程总结#

1. 用积性函数拟合出合适的 $g$。

2. 构造快速求 $g$ 前缀和的方法。

3. 预处理出 $h(p^c)$。

4. DFS 出所有 PN 并计算答案。

例题#

LibreOJ-6053 简单的函数#

除了 $2$ 以外的 $p$ 处点值是 $p-1$，先不考虑特例的话拟合成 $\varphi$ 是非常好的，在此基础上 $\varphi(2)=1$，所以将 $2$ 的倍数处拟合成 $3\varphi$。（肯定不是 $\varphi+2$ 吧）

$$g(n)=\begin{cases}3\varphi(n)&2\mid n\\\varphi(n)&2\nmid n\end{cases}$$

把 $\varphi$ 提出一个会好一点，顺便把系数也带出来：

$$l(n)=\begin{cases}\varphi(n)&2\mid n\\0&2\nmid n\end{cases}$$

于是 $g(n)=\varphi(n)+2l(n)$，也即 $S_g(n)=S_\varphi(n)+2S_l(n)$。

注意到 $S_l(n)$ 有意义的 $n$ 都是偶数，可以用一个 $S(n)=\sum_{i=1}^n \varphi(2i)$ 来代替，即 $S_l(n)=S\left(\left\lfloor\dfrac{n}{2}\right\rfloor\right)$。

把 $\varphi(2n)$ 详写成关于 $n$ 的：

$$\varphi(2n)=\begin{cases}2\varphi(n)&2\mid n\\\varphi(n)&2\nmid n\end{cases}$$

参照上面的化简方法，得到：

$$\varphi(2n)-\varphi(n)=\begin{cases}\varphi(n)&2\mid n\\0&2\nmid n\end{cases}$$

这个和 $l$ 一模一样，这里当然再把 $S_l$ 换成 $S$。

$$S(n)=S_\varphi(n)+S\left(\left\lfloor\dfrac{n}{2}\right\rfloor\right)$$

放回最初的式子 $S_g(n)=S_\varphi(n)+2S_h(n)=S_\varphi(n)+2S\left(\left\lfloor\dfrac{n}{2}\right\rfloor\right)$，这样一直把 $S$ 展开，就是一个倍增的形式，单次计算是 $O(\log n)$，所求点值都在杜教筛射程范围内。

构造的另一函数 $h$ 就是朴素的 $O(\sqrt{n}\log n)$ 预处理。

根据 PN 的数量级和预处理的复杂度，总体仍是 $O(n^{2/3})$。

Min_25 筛#

Min_25 筛是一种基于埃筛的筛法，可以对数论函数 $f(n)$ 求前缀和，要求 $f(p)$ 是关于 $p$ 的低次多项式，且 $f(p^k)$ 可快速计算。

构造质数点值前缀和函数#

在求解的开始，先用线性筛筛出小于等于 $\sqrt{n}$ 的质数，定义 $\mathrm{minp}(i)$ 表示 $i$ 的最小质因子。

设 $g(n)=\sum_{p\in \mathbf{P},p\le n} f(p)$，这个式子不好直接求，先设 $g_k(n)=\sum_{p\in \mathbf{P},p\le n} p^k$，即 $g(n)=\sum_{k} a_kg_k(n)$，其中 $a_k$ 是 $f(p)$ 这个关于 $p$ 的低次多项式 $k$ 次项系数。

设 $g_k(n,j)=\sum_{i=1}^n i^k[i\in\mathbf{P}\lor\mathrm{minp}(i)>p_j]$，这个定义相当于埃筛筛掉了含小于等于 $p_j$ 质因子的合数，容易发现找到最小的质数 $p_x$ 满足 $p_x\ge \sqrt{n}$，那么 $g_k(n)=g_k(n,x)$，因为此时不再有最小值因子大于 $p_x$ 的合数。

考虑递推，有：

$$g_k(n,j)= \begin{cases} g_k(n,j-1)&p_j^2\ge n\\ g_k(n,j-1)-p_j^k\left(g\left(\left\lfloor\dfrac{n}{p_j}\right\rfloor,j-1\right)-g(p_{j-1},j-1)\right)&p_j^2<n \end{cases}$$

第二个递推式实际上就是把 $\mathrm{minp}(i)=p_j$ 的数拆成 $p_j$ 和另一部分的乘积，这一部分就表现在 $g\left(\left\lfloor\dfrac{n}{p_j}\right\rfloor,j-1\right)$ 中，但这里面还包含了质数，所以要减去。

于是只关心 $O(\sqrt{n})$ 个块筛位置的点值，而 $p_{j-1}\le \sqrt{n}$，所以一定可以被块筛筛到，这样计算复杂度也是经典方法了，具体是：

$$\begin{aligned} &O\left(\sum_{i=1}^{\sqrt{n}} \pi(\sqrt{i})+\sum_{i=1}^{\sqrt{n}}\pi\left(\sqrt{\dfrac{n}{i}}\right)\right)\\ =&O\left(\sum_{i=1}^{\sqrt{n}}\pi\left(\sqrt{\dfrac{n}{i}}\right)\right)\\ =&O\left(\int_1^{\sqrt{n}} \dfrac{\sqrt{\tfrac{n}{x}}}{\log \sqrt{\tfrac{n}{x}}} \mathrm{d}x\right)\\ =&O\left(\dfrac{\sqrt{n}\displaystyle\int_1^{\sqrt{n}} \dfrac{1}{\sqrt{x}}\mathrm{d}x}{\log n}\right)\\ =&O\left(\dfrac{n^{3/4}}{\log n}\right) \end{aligned}$$

构造前缀和函数#

方法一：递归#

设 $S(n)=\sum_{i\le n}f(i)$，仿照上面，扩充 $S$ 的定义为 $S(n,j)=\sum_{i\le n}f(i)[\mathrm{minp}(i)>p_j]$，这样答案就是 $S(n,0)+1$。

直接递归，式子是：

$$S(n,j)=g(n)-g(p_j)+\sum_{\substack{p_j<p_k\le \sqrt{n}\\p_k^2\le p_k^{e+1}\le n}} f(p_k^e)S\left(\left\lfloor\dfrac{n}{p_k^e}\right\rfloor,k\right)+f(p_k^{e+1})$$

就是分质数合数来计算，后半部分相当于枚举了一个最小值因子及幂次，然后递归求解，$f(p_k^{e+1})$ 的加入源于 $S$ 并不计算 $1$ 的点值，因此要补充计算，而 $f(p_k)$ 已经在质数处算过了。

复杂度证明见论文，结论是复杂度为 $O(n^{1-\epsilon})$，但 $n\le 10^{13}$ 时为 $O\left(\dfrac{n^{3/4}}{\log n}\right)$。

实现常数很小，但只能求单点。

点击查看代码

inline int S1(ll N){
    N%=mod;
    return N*(N+1)%mod*inv2%mod;
}
inline int S2(ll N){
    N%=mod;
    return N*(N+1)%mod*(2*N+1)%mod*inv6%mod;
}
inline int f(ll N){
    N%=mod;
    return N*(N-1+mod)%mod;
}

ll n;
int lim;
int pr[maxn],cntpr;
bool vis[maxn];
ll id[maxn<<1];
int id1[maxn],id2[maxn];
int g1[maxn<<1],g2[maxn<<1];

inline int get_id(ll N){
    if(N<=lim) return id1[N];
    else return id2[n/N];
}

inline void linear_sieve(){
    for(int i=2;i<=lim;++i){
        if(!vis[i]) pr[++cntpr]=i;
        for(int j=1;j<=cntpr&&i*pr[j]<=lim;++j){
            vis[i*pr[j]]=1;
            if(i%pr[j]==0) break;
        }
    }
}

int S(ll N,int j){
    if(pr[j]>=N) return 0;
    int res=0;
    res=(g1[get_id(N)]-g1[get_id(pr[j])]+mod)%mod;
    for(int k=j+1;k<=cntpr&&1ll*pr[k]*pr[k]<=N;++k){
        ll now=pr[k];
        for(int e=1;now<=N;++e,now*=pr[k]){
            res=(res+1ll*f(now)*(S(N/now,k)+(e>1))%mod)%mod;
        }
    }
    return res;
}

int main(){
    scanf("%lld",&n);
    lim=sqrt(n)+1;
    linear_sieve();
    for(ll l=1,r;l<=n;l=r+1){
        r=n/(n/l);
        id[++id[0]]=n/l;
        if(n/l<=lim) id1[n/l]=id[0];
        else id2[n/(n/l)]=id[0];
        g1[id[0]]=(S1(n/l)-1+mod)%mod,g2[id[0]]=(S2(n/l)-1+mod)%mod;
    }
    for(int j=1;j<=cntpr;++j){
        for(int i=1;i<=id[0]&&1ll*pr[j]*pr[j]<=id[i];++i){
            g1[i]=(g1[i]-1ll*pr[j]*(g1[get_id(id[i]/pr[j])]-g1[get_id(pr[j-1])]+mod)%mod+mod)%mod;
            g2[i]=(g2[i]-1ll*pr[j]*pr[j]%mod*(g2[get_id(id[i]/pr[j])]-g2[get_id(pr[j-1])]+mod)%mod+mod)%mod;
        }
    }
    for(int i=1;i<=id[0];++i) g1[i]=(g2[i]-g1[i]+mod)%mod;
    printf("%d\n",(S(n,0)+1)%mod);
    return 0;
}

方法二：递推#

考虑把 $S$ 的状态设计和 $g$ 类比为 $S(n,j)=\sum_{i\le n} f(i)[i\in\mathbf{P}\lor\mathrm{minp}(i)>p_j]$。

这样目标是得到 $S(n,0)+1$，初始 $S(n,+\infty)=g(n)$（实际上就是只有质数部分，第二维取第一个大于 $\sqrt{n}$ 的质数编号即可）。

递推过程：

$$S(n,j)= \begin{cases} S(n,j+1)&p_{j+1}^2> n\\ S(n,j+1)+\displaystyle\sum_{p_{j+1}^2\le p_{j+1}^{e+1}\le n}f(p_{j+1}^e)\left(S\left(\left\lfloor\dfrac{n}{p_{j+1}^e}\right\rfloor,j+1\right)-S(p_{j+1},j+1)\right)+f(p_{j+1}^{e+1})&p_{j+1}^2\le n \end{cases}$$

这里对 $p_{j+1}^{e+1}$ 做限制的原因是保证了 $\left\lfloor\frac{n}{p_{j+1}^e}\right\rfloor\ge p_{j+1}$，从而消去 $S(p_{j+1},j+1)$ 是正确的。

复杂度同上，但常数较大，优越性在于递推算法求出了块筛。

点击查看代码

inline int S1(ll N){
    N%=mod;
    return N*(N+1)%mod*inv2%mod;
}
inline int S2(ll N){
    N%=mod;
    return N*(N+1)%mod*(2*N+1)%mod*inv6%mod;
}
inline int f(ll N){
    N%=mod;
    return N*(N-1+mod)%mod;
}

ll n;
int lim;
int pr[maxn],cntpr;
bool vis[maxn];
ll id[maxn<<1];
int id1[maxn],id2[maxn];
int g1[maxn<<1],g2[maxn<<1];
int S[maxn<<1];

inline int get_id(ll N){
    if(N<=lim) return id1[N];
    else return id2[n/N];
}

inline void linear_sieve(){
    for(int i=2;i<=lim;++i){
        if(!vis[i]) pr[++cntpr]=i;
        for(int j=1;j<=cntpr&&i*pr[j]<=lim;++j){
            vis[i*pr[j]]=1;
            if(i%pr[j]==0) break;
        }
    }
}

int main(){
    scanf("%lld",&n);
    lim=sqrt(n)+1;
    linear_sieve();
    for(ll l=1,r;l<=n;l=r+1){
        r=n/(n/l);
        id[++id[0]]=n/l;
        if(n/l<=lim) id1[n/l]=id[0];
        else id2[n/(n/l)]=id[0];
        g1[id[0]]=(S1(n/l)-1+mod)%mod,g2[id[0]]=(S2(n/l)-1+mod)%mod;
    }
    for(int j=1;j<=cntpr;++j){
        for(int i=1;i<=id[0]&&1ll*pr[j]*pr[j]<=id[i];++i){
            g1[i]=(g1[i]-1ll*pr[j]*(g1[get_id(id[i]/pr[j])]-g1[get_id(pr[j-1])]+mod)%mod+mod)%mod;
            g2[i]=(g2[i]-1ll*pr[j]*pr[j]%mod*(g2[get_id(id[i]/pr[j])]-g2[get_id(pr[j-1])]+mod)%mod+mod)%mod;
        }
    }
    for(int i=1;i<=id[0];++i) S[i]=g1[i]=(g2[i]-g1[i]+mod)%mod;
    for(int j=cntpr;j>=1;--j){
        for(int i=1;i<=id[0]&&1ll*pr[j]*pr[j]<=id[i];++i){
            ll now=pr[j];
            for(int e=1;now*pr[j]<=id[i];++e,now*=pr[j]){
                S[i]=(S[i]+1ll*f(now)*(S[get_id(id[i]/now)]-S[get_id(pr[j])]+mod)%mod+f(now*pr[j]))%mod;
            }
        }
    }
    printf("%d\n",(S[1]+1)%mod);
    return 0;
}

例题#

LibreOJ-6235 区间素数个数#

只需要做 Min_25 的第一步，求 $g$ 即可。

LibreOJ-6027 from CommonAnt 质数计数 I#

就是求模 $4$ 余 $1$ 的质数个数，注意到求 $g$ 部分是把 $\mathrm{minp}(i)=p_j$ 分成了 $p_j$ 和其他两部分，而 $g\left(\left\lfloor\dfrac{n}{p_j}\right\rfloor,j-1\right)$ 中每个数的余数乘上 $p_j$ 的余数才是得到 $\mathrm{minp}(i)=p_j$ 合数的余数。因此不能只求余数为 $1$ 的，计算时按余数分组即可。

LibreOJ-6028 from CommonAnt 质数计数 II#

与上一题类似

LibreOJ-6053 简单的函数#

先按 $f(p)=p-1$ 求 $g$，第二步算之前给前缀和再加上 $2$ 即可。

注意第一步只需要保证完全积性，第二步只需要保证积性，所以这样处理是正确的。

LibreOJ-6785 简单的函数 10^13 版#

要求块筛，用递推方法，需要大量卡常。

LibreOJ-6181 某个套路求和题#

$f(n)$ 不是积性函数，但是很简洁，容易发现质数位置是 $-1$，含平方因子位置是 $0$。

不含平方因子的合数位置，设共 $k$ 个质因子，那么只关心 $\mu(d)=-1$ 的位置，也就是 $k$ 中选出奇数个，方案数 $2^{k-1}$，因此不含平方因子的合数位置是 $1$。

只看合数位置可以拟合成 $\mu^2(n)$，这样质数位置从 $-1$ 变成了 $1$，答案就应该是：

$$\sum_{i\le n}\mu^2(i)-2\sum_{i\le n}[i\in \mathrm{P}]$$

前一个式子直接 Min_25 前缀和即可，容易发现后一个式子等价于第一步的 $g$ 函数。

参考资料#

狄利克雷卷积#

• 狄利克雷卷积 - command_block

杜教筛#

• OI Wiki

• 杜教筛 - command_block

PN 筛#

• OI Wiki

• PN 筛 - fjzzq2002

Min_25 筛#

• 洛谷日报 #406 Min_25 筛