mydcwfy's Blog

Min_25 筛

Posted on
Symbols count in article: 5.9k Reading time ≈ 13 mins.

非常 NB 的一个筛法，虽然复杂度可能会趋近线性 $O(n ^ {1 - \epsilon})$，但是在 OI 范围内非常不错。

1. 算法介绍

要么是 $O(\dfrac{n ^ {\frac34}}{\log n})$，还有一种是 $O(n ^ {1 - \epsilon})$，反正可以在 $10^{10}$ 范围内跑 1s 左右。

2. 算法流程

感觉比杜教筛更好构造。

假设 $f(x)$ 的前缀和是我们要求的。

分为两步：先求出在当 $x \in P$（表示 $x$ 是质数）的 $f(x)$ 的前缀和，第二步通过质数的前缀和求出整个的前缀和。

1）求 $g(x) = f(x)[x \in P]$ 的前缀和

我们找到一个多项式 $g(x)$，使得 $x \in P$ 的时候 $f(x) = g(x)$。

定义 $\text{lpf}(x)$ 是 $x$ 的最小质因数，$p_x$ 表示第 $x$ 个质数。

先将多项式拆为多个 $x ^ k$，这样的话 $g(x)$ 就是一个完全积性函数了，即满足 $g(ab) = g(a) \times g(b)$。

我们定义 $h(n, k) = \sum_{i = 1}^n [i \in P \lor \text{lpf}(i) > p_k]g(i)$。

（显然状态数比 $n$ 还大）

我们考虑从 $h(…, k - 1)$ 转移到 $h(n, k)$。

考虑 $h(n, k - 1)$ 变为 $h(n, k)$，需要减去 $\text{lpf}(x) = p_k$ 的数。

如果我们将 $x$ 除以 $p_k$，又怎么样呢？

那么剩下的变为 $\text{lpf}(\dfrac x{p_k})\geq p_k$。至于为什么是 $\geq$，是因为可能有单个数可能有多个 $p_k$ 质因子。

这不和 $h(\dfrac n{p_k}, k - 1)$ 定义很像了吗！

对比一下：

其中为什么要加 $h(p_{k - 1}, p_{k - 1})$ 呢？因为 $h(\dfrac n{p_k}, k - 1)$ 包含了 $[1, \dfrac{n}{p_k}]$ 所有质数，而小于 $p_k$ 的质数我们前面的是没有计算到的，所以要加上 $[1, p_{k - 1}]$ 的所有质数的函数值，显然 $g(p_{k - 1}, p_{k - 1})$ 就可以满足要求。这样我们就可以得出：

又由于 $g$ 是完全积性函数（记得前面我们拆开多项式的目的吗？），我们可以得到这样的递推式：

按照这个递推，最后我们就可以得到所有质数位置的值的前缀和，即为 $h(n, k)$，其中 $p_k > \sqrt n$。

2）根据质数点的值求整个的前缀和

现在我们经过第一步的计算，已经知道了只有质数点的值 $g(x) = [x \in P]f(x)$ 的前缀，记为 $r(x)$。

我们定义 $S(n, k)$ 表示 $\displaystyle \sum_{i = 2}^n [\text{lpf}(i) > p_k]f(k)$。

类似刚刚推导 $h(n, k)$ 的方法，我们可以较为简单的求，注意是积性函数而不是完全积性函数。

我们考虑最小的质因数 $p_v > p_k$，然后枚举 $p_v$ 的指数。如是就可以从 $S(\dfrac n{p_v^e}, v)$ 转移，因为剩下的质因数都得大于 $p_v$。

注意我们的 $S(n, k)$ 没有定义 1，所以我们还要先加上 $> p_k$ 的质数的所有的和，这样我们就可以直接使用 $r(n) - h(p_k, k)$ 代替。

注意如果 $e \not= 1$，$p_v^e$ 本身就是合数，没有包含在 $r(n) - h(p_k, k)$ 中，所以我们要加上他们。

根据积性函数的定义，我们可以得到：

3. 代码实现

显然，如果一个数不是质数，那么 $\text{lpf}(n) \leq \sqrt n$，所以我们计算 $h(n, k)$ 的时候可以直接处理 $[1, \sqrt n]$ 的所有质数就可以了。

但是我们如果直接循环 $h(i, k)(i \in [1, n], p_k \leq \sqrt n)$，那么肯定暴力都跑不过。

可以发现，除了 $h(p_{k - 1}, k - 1)$，其余都是和 $n$ 有关。又有一个整除的小定理：

这说明，不管我们除什么，前面的 $h(n, k), h(\dfrac n{p_k}, k - 1)$ 都只会访问到 $\dfrac ni(i \in [1, n])$，而这个最多只有 $O(\sqrt n)$ 个取值。

这个实现可以直接整除分块，将可以访问到的值存到两个数组，可以表示为：

当然偷懒的话可以多一个 $\log$ 的 std::map，但是似乎也不难，直接写整除分块了（本身算法就有点卡，还是不要给自己添麻烦吧 qwq

1
2
3
4
5
6
7
8
9
10
11
using LL = long long;

inline int get_id(LL x){return x <= sq ? id1[x] : id2[n / x];}

for (LL l = 1, r, t; l <= n; l = r + 1)
{
	r = n / (n / l), t = n / r;
	a[++ tot] = t;//a 存的是这个位置存的是哪个数
	if (t <= sq) id1[t] = tot;
	else id2[n / t] = tot;
}

我们接着简单的分析复杂度：枚举质数 $O(\dfrac{\sqrt n}{\ln \sqrt n})$，枚举 $O(\sqrt n)$ 个数，前面时间复杂度为 $O(\dfrac n{\ln \sqrt n})$，显然跑不进 1s。

我们还有一个优化：在枚举 $h(n, k)$ 的时候，我们可以发现 $h(\dfrac n{p_k}, k - 1) - h(p_{k - 1}, k - 1)$ 有贡献时，肯定需要满足 $\dfrac{n}{p_k} \geq p_k$，否则 $h(\dfrac n{p_k}, k - 1)$ 就全是质数，并且只有 $[1, p_k)$ 中的质数，也就是 $h(p_{k - 1}, k - 1)$，那么后面贡献就是 0。所以我们只需枚举 $n$ 到 $p_k ^ 2$。

另外，$h(p_{k - 1}, k - 1)$ 最多也只有 $O(\dfrac{\sqrt n}{\ln \sqrt n})$ 个，我们直接预处理就可以了。这个我们记为 $sp_{k - 1}$。

（代码用的是 $g(n, k)$ 表示 $h(n, k)$，且只开了一维数组，因为访问的时候一定是标号小（即对应的 a[] 大）的从小的转移，可以用标号从小到大）

1
2
3
4
5
6
for (int i = 0; i < cnt; ++ i)
{
	LL le = (LL)prime[i] * prime[i];
	for (int j = 1; j <= tot && a[j] >= le; ++ j)
		g[j] -= g[get_id(a[j] / prime[i])] - sp[i - 1];
}

而对于后面是同样的，我们计算 $S(n, k)$ 时，如果 $p_v^2 > n$，只能枚举 $e = 1$，那么 $[e \not= 1] = 0$，$S(\dfrac{n}{p_v}, v)$ 由于 $\dfrac n{p_v} < p_v$ 所以就为 0，我们就无需枚举。

反正后面 $S(n, k)$ 是不需要记忆化的（跑得还挺快

最后注意 $S(n, 0) = \sum_{i = 2}^n f(i)$，返回这个，需要加 $f(1)$。

1
2
3
4
5
6
7
8
9
10
11
12
LL S(LL n, int i)
{
	if (prime[i] >= n) return 0;
	int now = get_id(n);
	LL res = (f[now] - sp[i] + Mod) % Mod;
	for (int nw = i + 1; nw < cnt && ((LL)prime[nw]) * prime[nw] <= n; ++ nw)
		for (LL now = prime[nw], k = 1; now <= n; k ++, now *= prime[nw]){
			LL t = now % Mod;
			res = (res + f(now) % Mod * ((S(n / now, nw) + (k != 1)) % Mod) % Mod) % Mod;
		}
	return res;
}

4. 例题

T1：区间素数个数

题目传送门 LOJ

观察到这个函数并不是积性函数，所以我们似乎没法使用杜教筛等算法求解。

观察标签可得，我们可以发现，$h(n, x)(p_x > \sqrt n)$ 就是 $\sum_{i = 1}^n[i \in P]g(i)$，我们直接让 $g(x) = 1$，即可求出 $[1, n]$ 的质数个数即可。

点击查看代码 
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
#include <iostream>
#include <cmath>

using LL = long long;
const int N = 1e6 + 10;
int id1[N], id2[N], tot, prime[N], cnt;
bool st[N];
LL g[N], n, sq, a[N], sp[N];

void sieve(int N)
{
	st[1] = 1;
	for (int i = 2; i < N; ++ i)
	{
		if (!st[i]) prime[cnt ++] = i;
		for (int j = 0; j < cnt && i * prime[j] < N; ++ j)
		{
			st[i * prime[j]] = true;
			if (i % prime[j] == 0) break;
		}
	}
}

inline int get_id(LL x){return x <= sq ? id1[x] : id2[n / x];}

//this problem has a times 0
void solveg()
{
	for (int i = 0; i < cnt; ++ i) sp[i] = i + 1;
	for (int i = 1; i <= tot; ++ i) g[i] = a[i] - 1;
	for (int i = 0; i < cnt; ++ i)
	{
		LL le = (LL)prime[i] * prime[i];
		for (int j = 1; j <= tot && a[j] >= le; ++ j)
			g[j] -= g[get_id(a[j] / prime[i])] - sp[i - 1];
	}
}

int main()
{
	std::cin >> n;
	sq = sqrt(n);
	sieve(sq + 1);
	for (LL l = 1, r; l <= n; l = r + 1)
	{
		r = n / (n / l);
		a[++ tot] = n / r;
		if (n / r <= sq) id1[n / r] = tot;
		else id2[n / (n / r)] = tot;
	}
	solveg();
	std::cout << g[get_id(n)] << std::endl;
	return 0;
}

T2：简单的函数

题目传送门 LOJ

首先有一个小学老师教的结论：偶质数只有 2。

那么，$f(i) = i\oplus 1 (i \in P)$ 就是：

我们就可以通过先筛 $g(i) = i,g(i) = 1$ 来得到 $f(i)(i \in P)$。

最后，我们计算 $\sum_{i = 1}^n f(i)$ 的时候，由于 $f(2) = 3$，所以我们在筛到 $2$ 的时候应该加上 2。

具体代码如下：

1
2
3
4
5
6
7
8
9
LL S(LL n, int x)//x 从 0  开始标号
{
	if (prime[x] >= n) return 0;
	LL res = (f[get_id(n)] - sp[x] + Mod + 2 * (x == -1 && n > 1)) % Mod;
	for (int nw = x + 1; (LL)prime[nw] * prime[nw] <= n && nw < cnt; ++ nw)
		for (LL k = 1, pk = prime[nw]; pk <= n; ++ k, pk *= prime[nw])
			res += (prime[nw] ^ k) * (S(n / pk, nw) + (k != 1)) % Mod;
	return res % Mod;
}

T3：DIVCNTK

题目传送门 SPOJ

要求 $\sum_{i = 1}^n d(i ^ k)\bmod 2^{64}$，其中 $d(x)$ 代表 $x$ 的约数个数。

我们可以简单的发现，$d(x ^ k)$ 是积性函数：$d(a^kb^k) = d(a^k)d(b^k) (a \perp b)$。

令 $f(x) = d(x ^ k)$，我们再尝试推一下 $f(p ^ c) = c(k + 1)$。发现这是一个简单的多项式，我们尝试拆开 $f(p^c) = c$，Min_25 筛即可。

代码不给了，留作练习题（逃