Min_25 筛入门

Lijiangjun4 · 2026-04-22 00:09:47 · 算法·理论

前言

校内训练时被抓起来讲解 Min_25 筛，口干舌燥讲解十分钟，忘了三次公式，结果没人听懂。遂写文纪念之。

本文尽量以比较简单易懂的方式讲解这个算法。虽然文章比较长，公式比较多，但只要一点一点看，总能学会的。

算法简介

Min_25 筛是在低于线性时间的复杂度内求出积性函数 f 前缀和的算法。

前置定义和知识

数论函数：定义域在正整数的函数。即对于数论函数 f(x)，要求 x 为正整数；
积性函数：如果对于互质的两个正整数 a，b 和数论函数 f 满足 f(ab)=f(a)f(b)，则称 f 为积性函数；
- 常见积性函数有欧拉函数 \varphi 和莫比乌斯函数 \mu。
完全积性函数：如果对于任意的两个正整数 a，b 和数论函数 f 满足 f(ab)=f(a)f(b)，则称 f 为完全积性函数；
- 常见完全积性函数有单位函数 \varepsilon(n)=[n=1]，常数函数 I(n)=1，单项式 f(x)=x^k。
埃拉托斯特尼筛法（埃筛）：一种在 O(n\log\log n) 的时间复杂度内求出 n 以内的所有质数的算法。其思路为从小到大，把每个质数的的倍数筛掉，剩下的就是质数。

记号约定

下文若无特殊说明，所有变量均为正整数，名为 p 的变量均属于质数集。

适用范围

数论函数 f 可以用 Min_25 筛求出前缀和，需要满足：

比如如果 $x$ 可以被表示为 $p^k$，则 $f(x)=x^2+x+1$，由于其中 $x^2$，$x$ 和 $1$ 都是完全积性函数，所以 $f$ 是符合这个要求的。

算法思路

第一部分

这一步，我们将从积性函数最基本的定义把 f 的前缀和进行简单拆解。

根据积性函数的定义，当我们求 f(n) 的值时，我们可以把它分解质因数，将最小质因子的那一项拆出来，递归求其剩余的部分，即 f(n)=f(\text{minp}(n)^c)f\left(\dfrac{n}{\text{minp}(n)^c}\right)（c 为 n 质因数分解中最小质因子的次数）。

比如 n=60=2^2\times 3\times 5，我们有 f(n)=f(60)=f(2^2)f(3\times 5)=f(4)f(15)。根据“适用范围”中的条件，f(\text{minp}(n)^c) 是容易求的。

根据上文性质，我们可以尝试设计状态。其中一维跟问题规模（即你要求前多少个函数值的前缀和）有关，一维和最小质因子有关。

令 S(n,j) 表示 2\sim n 中最小质因子大于 p_j 的数的 f 值之和，即：

S(n,j)=\sum_{i=2}^{n}[\text{minp}(i)>p_j]f(i)

那么最终的答案为 S(n,0)+f(1)。

读者可以在这里暂停，想一想怎么把 S 写成递归形式。

我们枚举这些 i 的最小质因子以及其次数，有：

\def\bracket#1{\left(#1\right)} \def\floor#1{\left\lfloor#1\right\rfloor} S(n,j)= \sum_{\substack{k>j \\ p_{k}^{2}\le n}} \sum_{\substack{c\ge 1 \\ p_{k}^{c}\le n}} f(p_{k}^{c})\bracket{ S\bracket{ \floor{\frac{n}{p_{k}^{c}}}, k } +[c>1] } +\sum_{\substack{k>j \\ p_k\le n}} f(p_k)

:::info[公式解释]{open}

首先，这个式子被分成了两个部分：以最后一步的加号为界，前半部分计算 \text{minp}(x)>p_j 的合数 x 的 f 值，后半部分计算 \text{minp}(x)>p_j 的质数 x 的 f 值。

对于前半部分，我们枚举合数的最小质因子（\displaystyle\sum_{\substack{k>j \\ p_{k}^{2}\le n}}）和最小质因子的指数（\displaystyle\sum_{\substack{c\ge 1 \\ p_{k}^{c}\le n}}）。由于是合数，我们要求 p_k^2\le n。

对于所有最小质因子项为 p_k^c 的合数 x，我们将 f(x) 拆成 f(p_k^c)f(\frac{x}{p_k^c}) 的形式（因为两者肯定互质），则：
\begin{aligned} &\sum_{x}f(x)\\ =&\sum_{x}f(p_k^c)f\left(\frac{x}{p_k^c}\right)\\ =&f(p_k^c)\sum_{x}f\left(\frac{x}{p_k^c}\right) \end{aligned}
考虑 \displaystyle\sum_{x}f\left(\frac{x}{p_k^c}\right) 等于什么。分类讨论：
- 当 \frac{x}{p_k^c}>1 时，由于 x 的最小质因子为 p_k，所以 \text{minp}\left(\frac{x}{p_k^c}\right)>p_k。由于 \frac{x}{p_k^c} 的最大值为 \lfloor\frac{n}{p_k^c}\rfloor（显然），则有 S(\lfloor\frac{n}{p_k^c}\rfloor,k)。
- 当 c\neq 1 时 \frac{x}{p_k^c} 才可以取到 1，因为 c=1 且 \frac{x}{p_k^c}=1 时 x 为质数。因此还有 [c>1]。
所以，\displaystyle\sum_{x}f\left(\frac{x}{p_k^c}\right)=S\left(\frac{n}{p_k^c}\right)+[c>1]。
对于后半部分，非常简单，枚举大于 p_j 的质数，把它们的 f 加起来。

:::

我们不妨设 q(n) 表示 \sum_{p\le n}f(p)，则改写原式后半部分为：

\def\bracket#1{\left(#1\right)} \def\floor#1{\left\lfloor#1\right\rfloor} S(n,j)=\sum_{\substack{k>j \\ p_{k}^{2}\le n}}\sum_{\substack{c\ge 1 \\ p_{k}^{c}\le n}}f(p_{k}^{c})\bracket{S\bracket{\floor{\frac{n}{p_{k}^{c}}},k}+[c>1]}+q(n)-\sum_{i=1}^{j}f(p_i)

现在，只要能求 q 的值，就能求 S 的值了。

第二部分

求 q 的过程，我们使用埃筛的思想。

由于 n 很大，我们不能直接对质数下手，而是考虑把 [1,n] 中的合数都筛掉，剩下的就是质数的答案。我们发现埃筛的过程，是从小到大枚举每一个质数，筛掉这个质数的所有倍数。我们可以利用这种思路求 q。

我们令 g(n,j) 表示 [2,n] 的数在第 j 轮埃筛后，剩下的数字的 f 值之和。剩下的数字要么是质数，要么是最小质因子大于 p_j 的合数。即：

g(n,j)=\sum_{i=2}^{n}[i\in\mathbb{P}\operatorname{or}\text{minp}(i)>p_j]f(i)

由定义得 q(n)=g(n,|\mathbb{P}\cap[1,\sqrt{n}]|)。

有了第一部分的经验，大家肯定迫不及待把 i 的最小质因子提出来了。但这里仍然有一个问题：f 并不是完全积性函数，当 f(n) 被 \text{minp}(n) 筛到时，f(\text{minp}(n))f\left(\frac{n}{\text{minp}(n)}\right)\neq f(n)。

注意到我们在这一部分的最终目的是求出所有质数的 f 值之和，其在合数处的取值我们是不关心的。由“适用范围”得，f(x) 在 x=p^c 处可以拆成 c 个完全积性函数之和，即：\displaystyle f(x)=\sum_{k=1}^{c}f'_k(x)，f' 为完全积性函数。我们可以分别对这些 f' 进行计算，最后拼回 f 的 q 值。

比如：如果 f(p^c)=p^{2c}+p^{c}，现在求 q(10)。我们可以令 f'(p^c) 先等于 p^{2c}，再等于 p^{c}，分别求出其 q 值为 2^2+3^2+5^2+7^2=87 和 2+3+5+7=17，则 q(n)=87+17=104。

既然 f 只有在 p^c 处的取值是我们关心的，我们可以更改 g(n,j) 的定义：[2,n] 的数第 j 轮埃筛后剩下的数的 f' 值（f' 为组成 f(p^c) 的某个 f'_k，你可以理解为我们要求 c 次 g）。

g(n,j)=\sum_{i=2}^{n}[i\in\mathbb{P}\operatorname{or}\text{minp}(i)>p_j]f'(i)

尝试从第 j-1 次埃筛推出第 j 次埃筛，我们有递推式：

\def\bracket#1{\left(#1\right)} \def\floor#1{\left\lfloor#1\right\rfloor} g(n,j)= \left\{ \begin{aligned} &g(n,j-1)&,p_{j}^{2}>n\\ &g(n,j-1)-f'(p_j)\bracket{g\bracket{\floor{\frac{n}{p_j}},j-1}-\sum_{i=1}^{j-1}f'(p_i)}&,p_{j}^{2}\le n \end{aligned} \right.

初值显然为 g(n,0)=\sum_{i=1}^{n}f'(i)。

:::info[解释]{open}

递推式的思路是从第 j-1 轮埃筛 g(n,j-1)，减掉第 j 轮埃筛涉及到的数（最小质因子为 p_j 的合数）的 f' 值，得到 g(n,j) 的值。

首先，显然地，当 p_j^2>n，没有最小质因子等于 p_j 的合数，g(n,j)=g(n,j-1)。

当存在 \text{minp}(x)=p_j 的合数 x 时，我们要减去这些合数的 f' 值之和。模仿 S 的求解过程，我们将这些合数的 f'(x) 值提取公因数 f'(p_j)，剩下部分 \text{minp}(\frac{x}{p_j}) 应当大于等于 p_j，即公式中的 \displaystyle g\left(\left\lfloor\frac{n}{p_j}\right\rfloor,j-1\right)，但是 \displaystyle g\left(\left\lfloor\frac{n}{p_j}\right\rfloor,j-1\right) 中还包含了一些质数，我们要减掉 \sum_{i=1}^{j-1}f'(p_i) 才行。

:::

观察 S 和 g 的第一维：每次递归往下递归都是除以一个数下取整。根据 \displaystyle\left\lfloor\frac{\left\lfloor\frac{n}{a}\right\rfloor}{b}\right\rfloor=\left\lfloor\frac{n}{ab}\right\rfloor，第一维我们一定可以表现为 \lfloor\frac{n}{t}\rfloor 的形式，这意味着第一维只有 O(\sqrt{n}) 种取值，内存问题就解决了。

复杂度

时间复杂度为 O(\frac{n^{\frac{3}{4}}}{\log n})，证明比较困难，在此不展开，感性理解质数个数不多，值域还是根号，复杂度应该是亚线性。空间复杂度为 O(\sqrt{n}) 显然。

参考代码

代码细节我放在注释了，代码难度相比数据结构都算小清新了。

:::info[P5325 【模板】Min_25 筛]

#include<bits/stdc++.h>
#define N 2000005
#define mod 1000000007
#define inv2 500000004
#define inv6 166666668
using namespace std;
//由于 f(p^c)=p^{2c}-p^c，所以若无特殊解释，下文名字里带数字 1 的数组计算 f'=p^c 时的情况，带数字 2 的数组计算 f'=p^{2c} 的情况
long long n,sqr;
bool p[N];
long long tot,pri[N],s1[N],s2[N];
void sieve()
{
    //筛 sqrt(n) 以内质数
    for(int i=2;i<N;i++)
    {
        if(!p[i]) pri[++tot]=i;
        for(int j=1;j<=tot&&i*pri[j]<N;j++)
        {
            p[i*pri[j]]=true;
            if(i%pri[j]==0) break;
        }
    }

    //这是用于第二部分 g 的递推式和第一部分 f 的递推式中，求 f' 在前若干个质数的取值之和
    for(int i=1;i<=tot;i++) s1[i]=(s1[i-1]+pri[i])%mod;
    for(int i=1;i<=tot;i++) s2[i]=(s2[i-1]+pri[i]*pri[i]%mod)%mod;
    return;
}
long long m,id1[N],id2[N],w[N],g1[N],g2[N];
//这是用于计算 g(n,0) 的初值
inline long long f1(long long x){x%=mod;return x*(x+1)%mod*inv2%mod;}
inline long long f2(long long x){x%=mod;return x*(x+1)%mod*(2*x%mod+1)%mod*inv6%mod;}//sum_{i=1}^{n} i^2 = n(n+1)(2n+1)/6，可用整数裂项证明

//这是用于 g 数组第一维的压空间
//根据上文，g 的第一维只有 n/1 ~ n/n 等取值，不同的只有 sqrt(n) 种
//我们给不同的 n/1 ~ n/n 编号，在整除分块下可做
//但下文代码中显然从 n/1 ~ n/n 映射到编号更频繁
//我们需要以 n/1 ~ n/n 为数组的下标
//但有的 n/t 太大了
//考虑到 min(t,n/t) <= sqrt(n)
//我们记两个数组：id1 和 id2
//如果 n/t 不是很大，直接作为 id1 的下标进行存储
//否则，应存到 id2 的 n/(n/t) 的位置上
//查询同理，getid 根据查询的 x 的大小判断返回哪个数组存储的编号
inline long long getid(long long x)
{
    if(x<=sqr) return id1[x];
    else return id2[n/x];
}

//S 的递归计算函数
long long S(long long n,long long k)
{
    //这显然吧
    if(pri[k]>n) return 0;

    //S 函数递推式的后半部分，记得在使用 g 数组时要 getid
    long long res=((g2[getid(n)]-g1[getid(n)]+mod)%mod-(s2[k]-s1[k]+mod)%mod+mod)%mod;

    //枚举 p_k,c
    for(int i=k+1;i<=tot&&pri[i]*pri[i]<=n;i++)
    {
        for(long long c=1,pw=pri[i];pw<=n;c++,pw*=pri[i])
        {
            //抄公式
            res=(res+pw%mod*(pw%mod-1)%mod*(S(n/pw,i)+(c>1))%mod)%mod;
        }
    }
    return res;
}
int main()
{
    sieve();
    scanf("%lld",&n);
    sqr=(long long)sqrt(n);

    //整除分块
    for(long long l=1,r;l<=n;l=r+1)
    {
        r=n/(n/l);
        w[++m]=n/l;//小巧思：注意到 w 实际上是从大到小的
        //顺手计算 g(n,0) 的初值
        g1[m]=(f1(w[m])-1+mod)%mod;
        g2[m]=(f2(w[m])-1+mod)%mod;
        if(w[m]<=sqr) id1[w[m]]=m;
        else id2[n/w[m]]=m;
    }

    //g 的递推计算
    //注意到 g 的第二维在计算 S 的时候不重要，所以我们滚掉第二维
    //发现 g(n,j) 的计算只依赖 $g(n/p_k,j-1)$，所以我们应从大到小枚举 g 的第一维来保证滚动数组的正确性
    //恰好 w 也是从大到小的，所以顺序遍历 w 数组即可
    for(int j=1;j<=tot;j++)//先滚第二维
    {
        for(int i=1;i<=m&&pri[j]*pri[j]<=w[i];i++)//再顺序遍历 w
        {
            g1[i]=(g1[i]-pri[j]*(g1[getid(w[i]/pri[j])]-s1[j-1]+mod)%mod+mod)%mod;
            g2[i]=(g2[i]-pri[j]*pri[j]%mod*(g2[getid(w[i]/pri[j])]-s2[j-1]+mod)%mod+mod)%mod;
        }
    }

    printf("%lld\n",(S(n,0)+1)%mod);
    return 0;
}