本来没有学习这种较难的算法的想法的,因为比赛也做不到这种难度的题, 但是最近打牛客多校02,有一题要求 [ 1 , n ] [1,n] [1,n] 中素数的个数,我以为是像莫反一样容斥,但是后面感觉不行。赛后知道是用 min_25 筛来求,赛时过了一车,因此我也不得不学习这个算法了。
我打算拿洛谷里面的模板来举例。
就是给你积性函数
f ( x ) = { k x = 1 g ( x ) = ∑ i = 0 k a i x i x = p c , p ∈ p r i m e , c ≥ 1 f ( x 1 ) f ( x 2 ) ( x 1 , x 2 ) = 1 f(x)=\begin{cases} k & x=1 \\ g(x)=\sum_{i=0}^ka_ix^i & x=p^c,p\in prime,c\ge 1 \\ f(x_1)f(x_2) & (x_1,x_2)=1 \end{cases} f(x)=⎩
⎨
⎧kg(x)=∑i=0kaixif(x1)f(x2)x=1x=pc,p∈prime,c≥1(x1,x2)=1
然后让你求
∑ i = 1 n f ( i ) \sum_{i=1}^nf(i) i=1∑nf(i)
其中 n ≤ 1 0 13 n\le10^{13} n≤1013
在牛客多校的那个题中,也有 n ≤ 1 0 9 n\le 10^9 n≤109
一般来讲,想解决这种问题都要从 i i i 最小的质因子入手,因为一个数要么是质数,要么最小的质因子 ≤ n 0.5 \le n^{0.5} ≤n0.5,而 n 0.5 n^{0.5} n0.5 一般都是 1 0 6 10^6 106~ 1 0 7 10^7 107 数量级,可以直接用线性筛得到,也能全部枚举一次。因此总的思路大概就是先求出 ∑ i = 1 n g ( i ) \sum_{i=1}^ng(i) ∑i=1ng(i) 然后通过枚举最小的质因子来修正有超过一个质因子的数的贡献,这样能够很好的解决大质数无法直接得到的问题。
如果没有学过要怎么求 ∑ i = 1 n i k \sum_{i=1}^ni^k ∑i=1nik,可以看我的博客
基于这个大的思路,我自己感觉 min_25 筛的思路是非常自然,很容易就能记下来。
对于式子
∑ i = 1 n f ( i ) \sum_{i=1}^nf(i) i=1∑nf(i)
先将 1 1 1 去掉,因为其不包含任何素数,计算的时候容易出现奇怪的错误,虽然按照积性函数的定义, f ( 1 ) = 1 f(1)=1 f(1)=1。
于是变成
f ( 1 ) + ∑ i = 2 n f ( i ) f(1)+\sum_{i=2}^nf(i) f(1)+i=2∑nf(i)
自己感觉 min_25 筛很关键的一点是其能够快速求出
g ( b ) = ∑ p ∈ p r i m e , p ≤ b p k g(b)=\sum_{p\in prime,p\le b}p^k g(b)=p∈prime,p≤b∑pk其中 b = ⌊ n a ⌋ b=\left\lfloor \frac{n}{a} \right\rfloor b=⌊an⌋,这也意味着 b b b 的个数只有 O ( n 0.5 ) O(n^{0.5}) O(n0.5)
虽然只能求出若干次方的和,但只要把每一项都求一次就可以得到原函数了。
然后剩下的和数就能通过枚举小最小的质因子来考虑
因此 min_25 筛有两步,分别是求出质数出的函数和,求出所有的函数和
第一步:求出质数处的函数和
这一步的思路就是像我刚才说的,先求出所有的函数和,然后通过枚举最小的质因子把和数的贡献剪掉即可。
我们从小到大枚举素数,并且同时维护所有的 g ( b ) g(b) g(b),假设当前的素数为 p p p,假设其为第 i i i 个素数。那么显然有 p ≤ n 0.5 p\le n^{0.5} p≤n0.5。 我们记第 i i i 个素数为 p i p_i pi
对于 g ( b ) g(b) g(b),我们应该减去 1 1 1 ~ b b b 中最小质因子为 p p p 的和数,如果我们直接固定 p p p,也就是只需要考虑 1 1 1 ~ ⌊ b p ⌋ \left\lfloor \frac{b}{p} \right\rfloor ⌊pb⌋ 中最小质因子 ≥ p \ge p ≥p 的数,然后你就会发现这个数刚好就是 g ( ⌊ b p ⌋ ) − ∑ j = 1 i − 1 p j k g(\left\lfloor \frac{b}{p} \right\rfloor)-\sum_{j=1}^{i-1}p_j^k g(⌊pb⌋)−∑j=1i−1pjk,(因为 g g g 函数中包含 < p <p <p 的素数)后者只需要线性筛的时候预处理即可。
然后容斥完以后我们就可以得到所有的 g ( b ) g(b) g(b) 了。
第二步:求出原函数
这一步的思路是也是拆分成质数贡献与和数贡献,质数贡献可以直接通过前面的 g g g 函数得到,和数贡献可以通过枚举最小质因子不断递归求解。
我们另 S ( n , x ) S(n,x) S(n,x) 表示 1 1 1 ~ n n n 中最小质因子 > p x >p_x >px 的数的贡献之和,我们认为质数本身就是其最小质因子。
然后有
S ( n , x ) = g ( n ) − ∑ i = 1 x f ( p i ) + ∑ p i c ≤ n , i > x f ( p i c ) ( S ( ⌊ n p i c ⌋ , i + 1 ) + [ c > 1 ] ) S(n,x)=g(n)-\sum_{i=1}^xf(p_i)+\sum_{p_i^c\le n,i>x}f(p_i^c)(S(\left\lfloor\frac{n}{p_i^c}\right\rfloor,i+1)+[c>1]) S(n,x)=g(n)−i=1∑xf(pi)+pic≤n,i>x∑f(pic)(S(⌊picn⌋,i+1)+[c>1])
因为 c > 1 c>1 c>1 的时候 p i c p_i^c pic 本身也能够产生贡献。
以下是我洛谷模板的代码
#include <bits/stdc++.h>
#define ll long long
using namespace std;
const int mod = 1e9 + 7;
int power(int a, int b) {
int ret = 1;
while (b) {
if (b & 1) ret = 1ll * ret * a % mod;
a = 1ll * a * a % mod; b >>= 1;
}
return ret;
}
const int inv2 = power(2, mod - 2), inv3 = power(3, mod - 2);
const int N = 1e6 + 10;
int cnt, p[N]; bool v[N];
int sp1[N], sp2[N];
ll n, Sqr, w[N];
int tot, g1[N], g2[N], ind1[N], ind2[N];
void pre(int n) {
v[1] = 1;
for (int i = 2; i <= n; i++) {
if (!v[i]) {
p[++cnt] = i;
sp1[cnt] = (sp1[cnt - 1] + i) % mod;
sp2[cnt] = (sp2[cnt - 1] + 1ll * i * i) % mod;
}
for (int j = 1; j <= cnt && i * p[j] <= n; j++) {
v[i * p[j]] = 1;
if (i % p[j] == 0) break;
}
}
}
int S(ll x, int y) {
if (p[y] >= x) return 0;
ll k = x <= Sqr ? ind1[x] : ind2[n / x];
int ans = (1ll * g2[k] - g1[k] + mod - (sp2[y] - sp1[y]) + mod) % mod;
for (int i = y + 1; i <= cnt && 1ll * p[i] * p[i] <= x; i++) {
ll pe = p[i];
for (int e = 1; pe <= x; e++, pe *= p[i]) {
ll xx = pe % mod;
ans = (ans + xx * (xx - 1) % mod * (S(x / pe, i) + (e != 1)) % mod) % mod;
}
}
return ans % mod;
}
int main() {
// freopen("in.in", "r", stdin);
ios::sync_with_stdio(false);
cin.tie(nullptr);
cin >> n;
Sqr = sqrt((long double)n);
pre(Sqr);
for (ll i = 1, j; i <= n; i = j + 1) {
j = n / (n / i);
w[++tot] = n / i;
int now = w[tot] % mod;
g1[tot] = 1ll * now * (now + 1) / 2 % mod - 1;
g2[tot] = 1ll * now * (now + 1) / 2 % mod * (2 * now + 1) % mod * inv3 % mod - 1;
if (n / i <= Sqr) ind1[n / i] = tot;
else ind2[n / (n / i)] = tot;
}
for (int i = 1; i <= cnt; i++) {
for (int j = 1; j <= tot && 1ll * p[i] * p[i] <= w[j]; j++) {
ll k = w[j] / p[i] <= Sqr ? ind1[w[j] / p[i]] : ind2[n / (w[j] / p[i])];
g1[j] = (g1[j] - 1ll * p[i] * (g1[k] - sp1[i - 1] + mod) % mod + mod) % mod;
g2[j] = (g2[j] - 1ll * p[i] * p[i] % mod * (g2[k] - sp2[i - 1] + mod) % mod + mod) % mod;
}
}
cout << (S(n, 0) + 1) % mod << "\n";
return 0;
}
时间复杂度
第一步的时间复杂度是 O ( n 0.75 ln n ) O(\frac{n^{0.75}}{\ln n}) O(lnnn0.75) 是没有争议的,有一些人底下写的是 log \log log 也是对的,因为只有常数倍的差别
对于第一重循环,因为素数的密度的是 1 ln n \frac{1}{\ln n} lnn1,所以可以提出一个 1 ln n \frac{1}{\ln n} lnn1
然后变为这个式子 O ( ∑ i = 1 n i + n i ) = O ( ∑ i = 1 n n i ) = O ( n ∑ i = 1 n 1 i ) = O ( n ∫ 1 n 1 x d x ) = O ( n n ∣ 1 n ) = O ( n 0.75 ) O(\sum_{i=1}^{\sqrt n}\sqrt i+\sqrt \frac{n}{i})=O(\sum_{i=1}^{\sqrt n}\sqrt \frac{n}{i})=O(\sqrt n\sum_{i=1}^{\sqrt n}\sqrt \frac{1}{i})=O(\sqrt n\int_{1}^{n} \frac{1}{\sqrt x}\, dx)=O(\sqrt n\sqrt n|_{1}^{\sqrt n})=O(n^{0.75}) O(∑i=1ni+in)=O(∑i=1nin)=O(n∑i=1ni1)=O(n∫1nx1dx)=O(nn∣1n)=O(n0.75)
第二步我在网上没有找到具体的证明,自己也感觉比较难证明,但是听说是一般情况是 O ( n 1 − δ ) O(n^{1-\delta}) O(n1−δ),但是如果 n ≤ 1 0 13 n\le 10^{13} n≤1013,就满足 O ( n 0.75 ln n ) O(\frac{n^{0.75}}{\ln n}) O(lnnn0.75),基本适合OI的范围。
学完 min_25 筛,就会发现前面那个求素数个数的问题实在是太简单了,连 S S S 函数都不用算,直接把 g g g 函数计算以下就行了,然后把多项式设置为 f ( x ) = 1 f(x)=1 f(x)=1 即可。