UOJ Logo analysis的博客

博客

浅谈min25筛

2023-07-26 14:19:06 By analysis

更新日志

2023-05-31

博客第一版

2023-06-23

更新了博客中个别地方的表达。

对“例题导入”部分作了补充。

(想加几道例题,但快生地会考了,不得不先把注意力放在生地上 qwq,只能先咕咕咕了)

2023-07-01

生地考完了,期末考完了,放一天假(今天)然后就回学校竞赛。

根据大佬 RP_INT_MAX 的要求,暂时先加一道例题,随便简单补充了下 $g$ 的存储部分的讲解。

(等到以后轻松一点再考虑要不要再加例题)

添加了目录


目录

0x00 前言

0x01 符号定义

0x02 min25 筛的用处

0x10 第一层:求解 $\sum_{i=1}^{n}f(i)$

0x11 第二层:求解 $\sum_{i=2}^{n}[i \in P]f(i)$

0x12 回到第一层:求解 $\sum_{i=1}^{n}f(i)$

0x20 例题导入

0x21 补充例题

0x30 结论汇总

0x40 参考资料


0x00 前言

由于本蒟蒻的搜索能力较差,很多地方在学习是没有找到明确的推导过程或定义的讲解,所以只好手推或自行定义,如有错误,欢迎指正。

本篇主要是在讲 min25 筛的各条公式如何无脑的用推式子的方式推出,所以会有很多推式子的内容,我会在结论前标注以供公式恐惧症患者快速查看。

0x01 符号定义

下文中意义不明的符号,都可以在这找到。

符号/函数 含义
$P$ 质数集合
$p_j$ 第 $j$ 小的质数,特殊的,$p_0 = 0$。
$lpf_i$ $i$ 的最小质因子
$P'$ 质数集合的子集,包含 $[1,n]$ 内所有质数
$\vert P'\vert$ $P'$ 集合的大小(即 $n$ 内的质数个数)

0x02 min25 筛的用处

min25 筛被用于处理一类积性函数的前缀和。

也即,对于一个函数 $f$,可以使用 min25 筛的先决条件有:

一、函数 $f$ 为积性函数。

二、若 $p$ 为质数,则 $f(p^k)$ 可以写成一个低阶多项式形式(低阶没有明确的定义,但要求你能在代码中敲出来)。

三、若 $p$ 为质数,则 $f(p)$ 可以快速求解(快速求解并没有明确的定义,但通常是 $\Theta(1)$ 的)。

0x10 第一层:求解 $\sum_{i=1}^{n}f(i)$

如果你认真地阅读了 min25 筛的用处,你会发现使用条件几乎都是在限制质数,这是因为在 min25 筛的处理中需要用到大量的质数的性质。

我们可以将上式按照代入的自变量划分成三个部分:质数,合数与 1。

$$ \sum_{i=1}^{n}f(i)=\sum_{i=2}^{n}[i \in P]f(i) + \sum_{i=2}^{n}[i \not\in P]f(i) + f(1) $$

合数想与质数挂上钩,最好的方法就是分解质因数。

我们可以枚举一个合数的最小质因子及其次数。

由于 $f$ 是积性函数,所以可以把他所有的最小质因子提出(显然,这个最小质因子不会超过 $\sqrt{n}$)。

$$ =f(1) + \sum_{i=2}^{n}[i \in P]f(i) + \sum_{p_k,1 \leq p_k^a \leq n}^{\sqrt{n}}f(p_k^a)(\sum_{i=2}^{\lfloor\frac{n}{p_k^a}\rfloor}[lpf_i>p_k]f(i) + [a \not= 1]f(1)) $$

第二个求和符号枚举 $k$ 与 $a$。

上式中,为了不出现 $a=1$ 且 $i=1$ 的情况(这样子乘起来就是质数了,但这里只能是合数),我们令 $i$ 从 2 开始枚举,如果 $a \not= 1$ 再补上 $f(1)$ 的情况。

我们不妨设 $S(n,j) = \sum_{i=2}^{n}[lpf_i>p_j]f(i)$,那么原式(我们要求的前缀和)就等于 $S(n,0) + f(1)$(任何数的最小质因子都大于 0)。

根据上面的式子,我们可以依照相同的思路得到 $S$ 的递推式:

一、分质数与合数(这里定义的 $S$ 是从二开始枚举的,所以不需要 1)。

$$ S(n,j) = \sum_{i=2}^n[i \in P \&\& i > p_j]f(i) + \sum_{i=2}^n[i \not\in P]f(i) $$

二、合数分解质因数

$$ S(n,j) = \sum_{i=2}^{n}[i \in P \&\& i > p_j]f(i) + \sum_{p_k,1 \leq p_k^a \leq n}^{\sqrt{n}}f(p_k^a)(S(\lfloor\frac{n}{p_k^a}\rfloor,k) + [a \not= 1]f(1)) $$

还可以再展开一下(结论一):

$$ S(n,j) = \sum_{i=2}^{n}[i \in P]f(i) - \sum_{i=2}^{p_j}[i \in P]f(i) + \sum_{p_k,1 \leq p_k^a \leq n}^{\sqrt{n}}f(p_k^a)(S(\lfloor\frac{n}{p_k^a}\rfloor,k) + [a \not= 1]f(1)) $$

至此,我们就得到了答案计算的递归式了。

接下来就需要处理前面的质数部分了。

0x11 第二层:求解 $\sum_{i=2}^{n}[i \in P]f(i)$

显然,1 不是质数。

所以原式可以变成:

$$ \sum_{i=1}^{n}[i \in P]f(i) $$

为了解决它,我们引入一个函数:

$$ g(n,j) = \sum_{i=1}^{n}[i \in P || lpf_i > p_j]i^k $$

那么,$\sum_{i=1}^{n}[i \in P]i^k = g(n,|P'|)$。

由于 $\sum_{i=1}^{n}[i \in P]f(i)$ 中的 $i$ 全是质数,根据限制条件二,我们可以将 $f(i)$ 展开,然后用 $i^k$ 的线性组合表示。


可以举个例子(究竟如何线性组合需要依照具体题目进行分析):

$$ f(i) = i^2 - i (i \in P) $$

它由 $i^2$ 与 $i$ 线性组合(减法)而成。

于是,我们对 $i^2$ 与 $i$ 分别求 $g$,再相减。


接下来考虑递推 $g$ 。

显然根据定义有:

$$ g(n,j) = g(n,j-1) - \sum_{i=1}^{n}[lpf_i = p_j]i^k $$

如果 $i$ 为 $lpf = p_j$ 的最小合数,则 $i = p_j^2$。

如果 $i = p_j^2 > n$,显然 $g(n,j) = g(n,j-1)$。

否则继续讨论。

$i^k$ 是完全积性函数,可以直接提出质因子 $p_j$。

$$ g(n,j) = g(n,j-1) - p_j^k\sum_{i=1}^{n}[lpf_i \geq p_j](\frac{i}{p_j})^k $$

$$ = g(n,j-1) - p_j^k\sum_{i=1}^{n}[lpf_i > p_{j-1}](\frac{i}{p_j})^k $$

更改求和指标,使枚举 $i$ 变成枚举 $\frac{i}{p_j}$。

$$ = g(n,j-1) - p_j^k\sum_{i=1}^{\frac{n}{p_j}}[lpf_{ip_j} > p_{j-1}](\frac{ip_j}{p_j})^k $$

只要 $lpf_i > p_{j-1}$,$lpf_{ip_j} > p_{j-1}$ 一定成立。

$$ = g(n,j-1) - p_j^k\sum_{i=1}^{\lfloor\frac{n}{p_j}\rfloor}[lpf_{i} > p_{j-1}]i^k $$

将后面向 $g$ 的定义转化:

$$ = g(n,j-1) - p_j^k(\sum_{i=1}^{\lfloor\frac{n}{p_j}\rfloor}[i \in P||lpf_{i} > p_{j-1}]i^k - \sum_{i=1}^{\lfloor\frac{n}{p_j}\rfloor}[i \in P]i^k) $$

显然上式是有错的,对于大于 $p_{j-1}$ 的质数,会满足条件 $lpf_{i} > p_{j-1}$,不应重复筛去,更正为:

$$ = g(n,j-1) - p_j^k(\sum_{i=1}^{\lfloor\frac{n}{p_j}\rfloor}[i \in P||lpf_{i} > p_{j-1}]i^k - \sum_{i=1}^{p_{j-1}}[i \in P]i^k) $$

我们为了方便,再加入一个定义:

$$ sp_i = \sum_{i=1}^{p_i}[i \in P]i^k $$

他经过线性组合,得到 $Sp_i$(下文会结合例题详细讲解)。

于是:

$$ = g(n,j-1) - p_j^k(g(\lfloor\frac{n}{p_j}\rfloor,j-1) - sp_{j-1}) $$

总结一下(结论二):

$$ g(n,j) = \begin{Bmatrix} g(n,j-1) &&&&& p_j^2 > n\\ g(n,j-1) - p_j^k(g(\lfloor\frac{n}{p_j}\rfloor,j-1) - sp_{j-1}) &&&&&p_j^2\leq n \end{Bmatrix} $$

通过线性组合 $g$,我们可以得到 $G$(下文会结合例题详细讲解)。

0x12 回到第一层:求解 $\sum_{i=1}^{n}f(i)$

我们此前已经推出的结论:

$$ S(n,j) = \sum_{i=2}^{n}[i \in P]f(i) - \sum_{i=2}^{p_j}[i \in P]f(i) + \sum_{p_k,1 \leq p_k^a \leq n}^{\sqrt{n}}f(p_k^a)(S(\lfloor\frac{n}{p_k^a}\rfloor,k) + [a \not= 1]f(1)) $$

带入上面求出的结果(最终结论):

$$ S(n,j) = G(n,|P'|) - Sp_j + \sum_{p_k,1 \leq p_k^a \leq n}^{\sqrt{n}}f(p_k^a)(S(\lfloor\frac{n}{p_k^a}\rfloor,k) + [a \not= 1]f(1)) $$

自此,我们已经完成了大部分的理论推导,接下来我会结合例题进行讲解。

0x20 例题导入

例题解析

请先读题:传送门

首先,他肯定适用 min25 筛。

我们可以简单的改动一下题中的式子(展开为多项式):

$$ f(p^k) = p^k(p^k-1) = (p^k)^2 - p^k $$

显然他可以由 $i^k$ 线性组合得到(设 $F_i(x) = x^i$):

$$ f(p^k) = F_2(p^k) - F_1(p^k) $$

于是,我们可以分别求出 $F_2$ 和 $F_1$ 的结果,相减(即这道题的线性组合策略)得到 $f$ 的结果。

分析完题目,我们开始做:

第一步是要处理出 $Sp_i$。

$Sp_i$ 为 $sp_i$ 线性组合得到。

根据定义:

$$ sp_i = \sum_{i=1}^{p_i}[i \in P]F_k(i) $$

于是:

$$ Sp_i = \sum_{i=1}^{p_i}[i \in P]f(i) = \sum_{i=1}^{p_i}[i \in P]F_2(i)+F_1(i) $$

能这样拆是因为 $i$ 必定为质数。

容易发现,这里 $Sp_i$ 的 $p_i$ 必定在 $\sqrt{n}$ 范围内。

我们可以使用线性筛处理:

int pri[N],cnt;
int vis[N];
int sp1[N],sp2[N];//sp1存F_1的sp,sp2存F_2的sp
void init(int maxn)
{
    for(int i=2;i<=maxn;i++)
    {
        if(!vis[i])
        {
            pri[++cnt] = i;
            sp1[cnt] = (sp1[cnt-1] + i)%mod;
            sp2[cnt] = (sp2[cnt-1] + i * i)%mod;
        }
        for(int j=1;j<=cnt;j++)
        {
            if(i * pri[j] > maxn)break;
            vis[i * pri[j]] = 1;
            if(i % pri[j] == 0)break;
        }
    }
}

接下来处理 $G$。

由于我们使用 $G$,必然是这个形式的:$G(x,|P'|)$

所以第二个参数可以滚动掉。

但是由于 $n \leq 10^{10}$,不可以直接存起来。

我们发现,$x$ 只有可能是这样的 $\lfloor\frac{n}{a}\rfloor$。

只有 $\sqrt{n}$ 个。

但是如果直接存,下标还是到 $n$。

所以考虑这样的储存策略:

使用离散化的方法,使用序号 $tot$ 表示一个数。

使用 $ind1$ 与 $ind2$ 建立数字到序号的映射,使用 $w$ 建立序号到数字的映射。

当数字 $w[tot] = \frac{n}{x}$ 小于(等于) $\sqrt{n}$ 时,$x$ 必然大于(等于) $\sqrt{n}$。

所以令 $ind1[w[tot]] = tot$。

因为存在 $ind1$ 中的 $w[tot] \leq \sqrt{n}$,所以 $ind1$ 下标小于 $\sqrt{n}$。

反之则 $x$ 必然小于 $\sqrt{n}$,故令 $ind2[x] = tot$。

因为存在 $ind2$ 中的 $x \leq \sqrt{n}$,所以 $ind2$ 下标小于 $\sqrt{n}$。

代码如下:

其中 $f1$ 为计算 $g_1(x,0) = \sum_{i=1}^{x}[i \in P || lpf_i > 0]f_1(i) = \sum_{i=1}^{x}i = \frac{x(x+1)}{2}$

$f_2$ 计算 $g_2(x,0) = \sum_{i=1}^{x}i^2 = \frac{x(x+1)(2x+1)}{6}$

int f1(int x)
{
    x%=mod;
    return x * (x + 1)/2 % mod;
}
int f2(int x)
{
    x%=mod;
    return x * (x+1) % mod * ((2 * x  + 1)% mod) % mod * inv6 % mod;
}
int w[N],ind1[N],ind2[N],tot,g1[N],g2[N];
for(int l = 1,r;l<=n;l = r + 1) //枚举x可能出现的值(整除分块)
{
    r = n / (n / l);
    //w数组的含义:id(第几个可能的x值)转实际值
    w[++tot] = n/l; //表示第tot个值是n/l
    //计算g(x,0)
    //实际处理中,我们一般不含f(1)(这里f(1)=1),只要最后加回即可
    g1[tot] = f1(w[tot])-1;
    g2[tot] = f2(w[tot])-1;
    //ind1[x]为x的id
    //ind2[x]为n/x的id
    if(w[tot] <= sqr) // 当x<=sqrt(n)时,存于ind1,否则存于ind2
        ind1[w[tot]] = tot;
    else
        ind2[n/w[tot]] = tot;
}

然后,根据公式递推 $g$:

int getid(int x) // x to id
{
    if(x <= sqr)
        return ind1[x];
    else
        return ind2[n/x];
}
for(int i=1;i<=cnt;i++)
{
    for(int j=1;j<=tot && pri[i] * pri[i] <= w[j];j++)
    {
        g1[j] =(g1[j] - pri[i] * (g1[getid(w[j]/pri[i])] - sp1[i-1]) % mod + mod)%mod;
        g2[j] =(g2[j] - pri[i] * pri[i] % mod * (g2[getid(w[j]/pri[i])] - sp2[i-1]) % mod+mod)%mod;
    }
}

最后计算 $S$:

int s(int x,int k)
{
    if(pri[k] > x)
        return 0;
    int ans = ((g2[getid(x)] - g1[getid(x)] + mod)%mod - (sp2[k] - sp1[k] + mod)%mod + mod)%mod;
    for(int i=k+1;i<=cnt && pri[i] * pri[i] <= x;i++)
    {
        for(int a = 1,p = pri[i];p <= x;p *= pri[i],a++)
        {
            ans =(ans + p % mod * ((p% mod-1) % mod * (s(x/p,i) + (a!=1)) % mod))%mod;
        }
    }
    return ans;
}

汇总代码:

#include<bits/stdc++.h>
#define int long long
const int N = (int)2e5+7,mod = (int)1e9+7;
const int inv2 = 500000004,inv6 = 166666668;
using namespace std;
int pri[N],cnt;
int vis[N];
int sp1[N],sp2[N];
void init(int maxn)
{
    for(int i=2;i<=maxn;i++)
    {
        if(!vis[i])
        {
            pri[++cnt] = i;
            sp1[cnt] = (sp1[cnt-1] + i)%mod;
            sp2[cnt] = (sp2[cnt-1] + i * i)%mod;
        }
        for(int j=1;j<=cnt;j++)
        {
            if(i * pri[j] > maxn)break;
            vis[i * pri[j]] = 1;
            if(i % pri[j] == 0)break;
        }
    }
}
int n,sqr;
int w[N],ind1[N],ind2[N],tot,g1[N],g2[N];
int getid(int x)
{
    if(x <= sqr)
        return ind1[x];
    else
        return ind2[n/x];
}
int f1(int x)
{
    x%=mod;
    return x * (x + 1)/2 % mod;
}
int f2(int x)
{
    x%=mod;
    return x * (x+1) % mod * ((2 * x  + 1)% mod) % mod * inv6 % mod;
}
int s(int x,int k)
{
    if(pri[k] > x)
        return 0;
    int ans = ((g2[getid(x)] - g1[getid(x)] + mod)%mod - (sp2[k] - sp1[k] + mod)%mod + mod)%mod;
    for(int i=k+1;i<=cnt && pri[i] * pri[i] <= x;i++)
    {
        for(int a = 1,p = pri[i];p <= x;p *= pri[i],a++)
        {
            ans =(ans + p % mod * ((p% mod-1) % mod * (s(x/p,i) + (a!=1)) % mod))%mod;
        }
    }
    return ans;
}
signed main()
{
    cin>>n;
    sqr = sqrt(n);
    init(sqr);
    for(int l = 1,r;l<=n;l = r + 1)
    {
        r = n / (n / l);
        w[++tot] = n/l;
        g1[tot] = f1(w[tot])-1;
        g2[tot] = f2(w[tot])-1;
        if(w[tot] <= sqr)
            ind1[w[tot]] = tot;
        else
            ind2[n/w[tot]] = tot;
    }
    for(int i=1;i<=cnt;i++)
    {
        for(int j=1;j<=tot && pri[i] * pri[i] <= w[j];j++)
        {
            g1[j] =(g1[j] - pri[i] * (g1[getid(w[j]/pri[i])] - sp1[i-1]) % mod + mod)%mod;
            g2[j] =(g2[j] - pri[i] * pri[i] % mod * (g2[getid(w[j]/pri[i])] - sp2[i-1]) % mod+mod)%mod;
        }
    }
    //f(1)=1
    cout<<(s(n,0)+1)%mod;
    return 0;
}

于是,你就拿到了 $100pts$。

对本例题的汇总

min25 筛的过程中,需要求解的有:

  • $Sp$ 函数

  • $G$ 函数

  • $S$ 函数

对于 $Sp$ 函数,我们可以利用线性筛预处理出来。

查看整个 min25 筛的过程,使用到的 $Sp$ 不超过 $\vert P' \vert$。

所以预处理出 $Sp_1$ ~ $Sp_{\sqrt{n}}$ 的值。

对于 $G$ 函数,我们只需要使用 $G(x,\vert P' \vert)$,所以第二维可以滚动掉。

对于 $x$ 的 $\sqrt{n}$ 个可能值,使用类似离散化的手法。

$W_x$ 为 $x$ 这个 $id$ 的真实值,$ind1/2$ 为值对应的 $id$。

然后根据推导出的公式递推。

对于 $S$ 函数,递归求解即可。

0x21 补充例题

Loj6235 区间素数个数

首先,我们按照需求,构造出需要求得的函数。

$$ ans = \sum_{i=1}^{n}[i \in P] $$

如果你对 min25 筛的过程足够熟悉,你就会发现

$$ \begin{Bmatrix} g(n,j) = \sum_{i=1}^{n}[i \in P || lpf_i > p_j]i^k &&&&& g的定义 \\ ans = \sum_{i=1}^{n}[i \in P] &&&&& 目标式子 \end{Bmatrix} $$

的相似之处。

不妨变换一下:

$$ ans = \sum_{i=1}^{n}[i \in P] = \sum_{i=1}^{n}[i \in P]i^0 = g(n,|P'|) , k=0 $$

于是我们很快就打出代码:

#include <bits/stdc++.h>
#define ll long long
using namespace std;
//变量名的意义和之前时一样的
ll n, sqr,w[1000005],g[1000005];
int pri[1000005], cnt;
int vis[1000005];
int ind1[1000005], ind2[1000005], tot;

void init(int maxn) {
    //k值为0,sp[x]即为1~p_x的质数个数
    for (int i = 2; i <= maxn; i++) {
        if (!vis[i]) {
            pri[++cnt] = i;
            //sp[cnt] = cnt;
            //没有必要多用一个数组储存sp
        }

        for (int j = 1; j <= cnt; j++) {
            if (i * pri[j] > maxn)
                break;
            vis[i * pri[j]] = 1;
            if (i % pri[j] == 0)
                break;
        }
    }
}

int getid(ll x) {
    if (x <= sqr)
        return ind1[x];
    else
        return ind2[n / x];
}
signed main() {
    scanf("%lld",&n);
    sqr = sqrt(n);
    init(sqr);
    for (ll l = 1, r; l <= n; l = r + 1) {
        r = n / (n / l);
        w[++tot] = n / l;
        g[tot] = w[tot] - 1;
        if (w[tot] <= sqr) {
            ind1[w[tot]] = tot;
        } else {
            ind2[n / w[tot]] = tot;
        }
    }
    for (int j = 1; j <= cnt; j++) {
        for (int i = 1; i <= tot; i++) {
            //k=0
            if ((ll)pri[j] * pri[j] > w[i])break;
            g[i] -= (g[getid(w[i] / pri[j])] - (j-1));
        }
    }
    printf("%lld",g[getid(n)]);
}

100pts

但请记住,在这里:

for (int j = 1; j <= cnt; j++) {
    for (int i = 1; i <= tot; i++) {
        //k=0
        if ((ll)pri[j] * pri[j] > w[i])break;
        g[i] -= (g[getid(w[i] / pri[j])] - (j-1));
    }
}

请记住别写成这样(60pts TLE 的痛):

for (int j = 1; j <= cnt; j++) {
    for (int i = 1; i <= tot; i++) {
        //k=0
        if ((ll)pri[j] * pri[j] <= w[i])
        g[i] -= (g[getid(w[i] / pri[j])] - (j-1));
    }
}

0x30 结论汇总

定义

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

$$g(n,j) = \sum_{i=1}^{n}[i \in P || lpf_i > p_j]i^k$$

$$sp_i = \sum_{i=1}^{p_i}[i \in P]i^k$$

公式

$$ S(n,j) = G(n,|P'|) - Sp_j + \sum_{p_k,1 \leq p_k^a \leq n}^{\sqrt{n}}f(p_k^a)(S(\lfloor\frac{n}{p_k^a}\rfloor,k) + [a \not= 1]f(1)) $$

$$ g(n,j) = \begin{Bmatrix} g(n,j-1) &&&&& p_j^2 > n\\ g(n,j-1) - p_j^k(g(\lfloor\frac{n}{p_j}\rfloor,j-1) - sp_{j-1}) &&&&&p_j^2\leq n \end{Bmatrix} $$

0x40 参考资料

analysis Avatar