更新日志

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 符号定义

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

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 参考资料

• 小蒟蒻yyb的博客