浅谈 min25 筛
analysis
·
·
个人记录
更新日志
原版写的太冗杂了，重写了一下。
有一题有版权问题，删掉了。
原版
目录
-
引出问题之前
-
解决 min25 筛应对的问题
- 引出问题
- 讨论质数部分
- 讨论合数部分
- 解决问题
- 代码实现
- 简单总结
-
例题与解析
引出问题之前
符号表
| 符号 |
意义 |
| p_i |
所有以 p 命名的变量默认为质数 |
| P |
质数集合，对应 n 的质数集合我们认为是所有小于等于 \sqrt{n} 的质数 |
| lpf_i |
i 的最小质因子 |
| [x] |
艾弗森括号，当命题 x 成立时为 1，否则为 0。 |
解决的问题
一类特殊的积性函数前缀和。
（积性函数在自变量为质数或其若干次方时可以使用关于质数的多项式表示）
以及一些有关质数的和式计算。
主要是利用了质数的性质快速处理。
解决 min25 筛应对的问题
以下所有式子汇总在简单总结小节。
引出问题
以洛谷模版题为例：
定义积性函数 f(x)，且 f(p^k)=p^k(p^k-1)（p 是一个质数），求
\sum_{i=1}^n f(i)
对 10^9+7 取模。
对于 100\% 的数据，保证 1\le n\le 10^{10}
首先，它是积性函数。
它把在自变量为质数或其若干次方时关于质数的多项式表示给出来了。
所以，这是道 min25 筛解决的问题。
可以分别讨论质数与合数。
\sum_{i=1}^{n}f(i) = \sum_{i=1}^{n}[i \in P]f(i) + \sum_{i=1}^{n}[i \not\in P]f(i)
讨论质数部分
\sum_{i=1}^{n}[i \in P]f(i)
为了解决这个问题，我们要引入一个神奇的东西：
g_k(n,j) = \sum_{i=1}^{n}[i \in P 或 lpf_i > j]F_k(i),F_k(i) = i^k
最重要的是，根据前面提出的“**积性函数**在**自变量为质数或其若干次方**时可以使用**关于质数的多项式**表示”这条条件中指出了 $f$ 可以表达为多项式。
也就是可以使用 $F_k$ 组合而成。
如本题：
$$
f(p^k) = p^k(p^k-1) = (p^k)^2 - p^k = F_2(p^k) - F_1(p^k)
$$
这样一来，我们可以分别计算 $F_2$ 与 $F_1$ 的前缀和，再相减得出答案。
根据上面的讨论，我们可以得出：
$$
\sum_{i=1}^{n}[i \in P]f(i) = g_2(n,|P|) - g_1(n,|p|)
$$
（显然 $n$ 内所有合数的 $lpf$ 不会大于 $|P|$）
然后我们回到 $g$。
$$
g_k(n,j) = \sum_{i=1}^{n}[i \in P 或 lpf_i > j]F_k(i),F_k(i) = i^k
$$
我有个无脑推式子的方法（详见原版），但现在再看觉得没那个必要了。
显然有：
$$
g_k(n,j) = g_k(n,j-1) - \sum_{i=1}^n[i \not\in P 且 lpf_i = j]i^k
$$
记 $sp_i = \sum_{i=1}^{p_i}[i \in P]i^k$。
我们考虑表示 $\sum_{i=1}^n[i \not\in P 且 lpf_i = p_j]i^k$。
去掉质数剩下部分可以提出 $p_j$，因为 $lpf=p_k$。
$F_k$ 是完全积性函数，所以我们强制提出 $p_j$，即：
$$
g_k(n,j)= g_k(n,j-1) - F_k(p_j)(g_k(\lfloor\frac{n}{p_j}\rfloor,j-1) - sp_{j-1})
$$
$lpf = j$，所以去掉 $sp_{j-1}$。
然后发现，当 $p_j^2 > n$ 时，也即没有 $lpf = j$ 的合数时，更新没有意义，所以（这可以在实现时剪枝，以免 TLE）：
$$
g_k(n,j) = \begin{Bmatrix}
g_k(n,j-1) &&&&& p_j^2 > n\\
g_k(n,j-1) - p_j^k(g_k(\lfloor\frac{n}{p_j}\rfloor,j-1) - sp_{j-1}) &&&&&p_j^2\leq n
\end{Bmatrix}
$$
## 讨论合数部分
这里不讨论 $1$。
看到式子：
$$
\sum_{i=1}^n f(i)
$$
式子里是 $1$ 到 $n$ 的所有值的前缀和，但我们只有自变量为质数或其若干次方时关于质数的多项式，没有合数的。
但是，我们可以质因数分解一下。
对于 $i = p_1^{a_1} p_2^{a_2} \cdots p_k^{a_k}$，有：
$$
f(i) = \prod_{j=1}^kf(p_j^{a_j})
$$
上式可以再改写一下，变成递归的形式，这是我们需要的形式。
我们枚举 $i$ 的最小质因子及次数 $c$：
$$
f(i) = f(lpf_i^c)f(\frac{i}{lpf_i^c})
$$
于是我们就可以愉快地表示合数了。
## 解决问题
然后原式就可以写成（枚举合数的最小质因子及次数）：
$$
\sum_{i=1}^{n}[i \in P]f(i) + \sum_{p_i^a \leq n,i \leq |p|}f(p_i^a)(\sum_{i=2}^{\lfloor\frac{n}{p_i^a}\rfloor}[lpf_i > p_i]f(i) + [a \not= 1]f(1)) + f(1)
$$
（如果 $a=1$ 且 $i=1$ 就变成质数了）
你会发现，若令
$$
S(n,j) = \sum_{i=2}^{n}[lpf_i>p_j]f(i)
$$
上式又可以递归：
$$
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))
$$
至此，我们就得到了答案计算的递归式了。
## 代码实现
这也是一个重点，光看推导是实现不出代码的。
先处理质数集与 $Sp$，用线性筛可以搞定。
```cpp
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(\lfloor\frac{n}{x}\rfloor,j)
$$
所以我们只需要存这些点的值即可。
根据我们学习整除分块时的经验，可以知道不超过 $O(\sqrt{n})$ 个点值。
> 简单证明
>
> 当 $x \leq \sqrt{n}$，$x$ 有 $\sqrt{n}$ 个取值，$\lfloor\frac{n}{x}\rfloor$ 不多于 $\sqrt{n}$ 个取值。
>
> 当 $x > \sqrt{n}$，$\lfloor\frac{n}{x}\rfloor < \sqrt{n}$，$\lfloor\frac{n}{x}\rfloor$ 不多于 $\sqrt{n}$ 个取值。
>
> 证毕。
根据上面的证明，我们想到离散化储存这些点值的方法。
当 $\lfloor\frac{n}{x}\rfloor \leq \sqrt{n}$ 时，以 $\lfloor\frac{n}{x}\rfloor$ 为关键字存在 `ind1`。
当 $\lfloor\frac{n}{x}\rfloor > \sqrt{n}$ 时，以 $x$ 为关键字存在 `ind2`。
```cpp
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];
```
```cpp
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(n,0)$。（上面的代码滚动掉了第二维）
然后根据公式求解 $g$。
```cpp
int getid(int x) // x to id
{
if(x <= sqr)
return ind1[x];
else
return ind2[n/x];
}
```
```cpp
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$。
```cpp
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;
}
```
完整代码：
```cpp
#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;
}
```
## 简单总结
### 定义
$$
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}
$$
# 例题与解析
## 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
$$
于是我们很快就打出代码：
```cpp
#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 && 1LL * pri[j] * pri[j] <= w[i]; i++) {
g[i] -= (g[getid(w[i] / pri[j])] - (j-1));
}
}
printf("%lld",g[getid(n)]);
}
```