Min_25筛

Min_25 筛

tags: 数学 数论

---

简介

Min_25 筛法用于求解一些积性函数的前缀和.
具体来说，对于积性函数的前缀和函数 \(S\)，Min_25 筛可以求出 \(S(n)\)，并同时求出 \(S(\lfloor\frac{n}{2}\rfloor),S(\lfloor\frac{n}{3}\rfloor),\cdots,S(1)\).

适用条件

设要求的函数为 \(F\).

1. 该函数为积性函数（个别非积性函数也可以）;
2. 若 \(p\) 为质数，\(F(p^e)\) 可以表示为一个可以快速求得的多项式.

记号

\(\mathbb{P}\) 表示质数集
\(\text{lpf}(n)\) 表示 \(n\) 的最小质因数
\(p_i\) 表示从小到大第 \(i\) 个质数(定义\(p_0 = 1\))

\(S_j(n) = \sum\limits_{i = 1,\text{lpf}(i)>p_j}^{n}F(i)\)

\(S_{\text{pri}}(n)=\sum\limits_{i=2,i\in \mathbb{P}}^{n}F(i)\)

求解方法

第一部分

我们发现答案 \(\sum\limits_{i=1}^n F(i)=S_0(n)+F(1)\)

\[\begin{align*} S_j(i) &= \sum\limits_{x = 1,\text{lpf}(x)>p_j}^{i}F(x)\\ &= \sum_{x = j+1}^{p_x \le i}\sum_{e=1}^{p_x^e \le i}F(p_x^e)(1+\sum_{y = 2,\text{lpf}(y)>p_x}^{\lfloor \frac{i}{p_i} \rfloor}F(y)) & \text{提出最小质因子} \\ &= \sum_{x = j+1}^{p_x \le i}\sum_{e=1}^{p_x^e \le i}F(p_x^e)([e>1]+\sum_{y = 2,\text{lpf}(y)>p_x}^{\lfloor \frac{i}{p_x} \rfloor}F(y))+\sum_{i=p_j+1,i\in \mathbb{P}}^{i}F(y) & \text{提出质数} \\ &= \sum_{x = j+1}^{p_x \le \sqrt{i}}\sum_{e=1}^{p_x^e \le i}F(p_x^e)([e>1]+\sum_{y = 2,\text{lpf}(y)>p_x}^{\lfloor \frac{i}{p_x} \rfloor}F(y))+\sum_{i=p_j+1,i\in \mathbb{P}}^{i}F(y) & p_x\text{上界变为}\sqrt i \\ &= \sum_{x = j+1}^{p_x \le \sqrt{i}}\sum_{e=1}^{p_x^e \le i}F(p_x^e)([e>1]+\sum_{y = 2,\text{lpf}(y)>p_x}^{\lfloor \frac{i}{p_x} \rfloor}F(y))+(S_{\mathtt{pri}}(i)-S_{\mathtt{pri}}(p_j)) & \text{代入}S_\text{pri} \\ &= \sum_{x = j+1}^{p_x \le \sqrt{i}}\sum_{e=1}^{p_x^e \le i}F(p_x^e)([e>1]+S_x(\lfloor \frac{i}{p_x^e}\rfloor))+(S_{\mathtt{pri}}(i)-S_{\mathtt{pri}}(p_j)) & \text{代入}S_x \end{align*} \]

该式可以递归求解(下文 P5325 代码 1)，也可以通过 \(S_{j}(i)\) 转移出 \(S_{j-1}(i)\) (P5325 代码 2)

\[\begin{align*} S_{j-1}(i)-S_{j}(i) &=\sum_{x = j}^{p_x \le \sqrt{i}}\sum_{e=1}^{p_x^e \le i}F(p_x^e)([e>1]+S_x(\lfloor \frac{i}{p_x^e}\rfloor))+(S_{\mathtt{pri}}(i)-S_{\mathtt{pri}}(p_{j-1}))\\&-\sum_{x = j+1}^{p_x \le \sqrt{i}}\sum_{e=1}^{p_x^e \le i}F(p_x^e)([e>1]+S_x(\lfloor \frac{i}{p_x^e}\rfloor))+(S_{\mathtt{pri}}(i)-S_{\mathtt{pri}}(p_j))\\ &=\sum_{e=1}^{p_j^e\le i}F(p_j^e)([e>1]+S_j(\lfloor \frac{i}{p_j^e}\rfloor))+(S_{\text{pri}}(p_j)-S_{\text{pri}}(p_{j-1})) & \text{第一层}\sum\text{只剩下}p_j \\ &=\sum_{e=1}^{p_j^e\le i}F(p_j^e)([e>1]+S_j(\lfloor \frac{i}{p_j^e}\rfloor))+F(p_j) & \text{消去}S_\text{pri}\\ &=\sum_{e=2}^{p_j^e\le i}F(p_j^e)+\sum_{e=1}^{p_j^e\le i}F(p_j^e)S_j(\lfloor \frac{i}{p_j^e}\rfloor)+F(p_j) & \text{拆开括号} \\ &=\sum_{e=1}^{p_j^e\le i}F(p_j^e)+\sum_{e=1}^{p_j^e\le i}F(p_j^e)S_j(\lfloor \frac{i}{p_j^e}\rfloor) & \text{合并} \end{align*} \]

式子的第一部分( \(\sum\limits_{e=1}^{p_j^e\le i}F(p_j^e)\) )可以快速地求和（暴力枚举 \(e\)）.

现在考虑第二部分：
令

\[T(n)=\sum_{e=1}^{p_j^e\le n}F(p_j^e)S_j(\lfloor \frac{n}{p_j^e}\rfloor) \]

由于对于 \(\forall x=p^e,p\in\mathbb{P}\)，有 \(F(x)\) 是多项式（不妨设 \(F(x)=a_0+a_1x^1+a_2x^2+\cdots+a_mx^m\)），那么我们把 \(T_j(i)\) 中的 \(F(p_j^e)\) 拆开：

设

\[T_{k}(i)=\sum_{e=1}^{p_j^e\le i}(p_j^e)^kS_j(\lfloor \frac{i}{p_j^e}\rfloor) \]

有

\[T(n)=\sum_{k=0}^{m}a_kT_{k}(n) \]

（我们只要求出每个 \(T_{j,k}\)，就可以求和得到 \(T_j\)）

对于 \(T_k\) 考虑从 \(T_k(\lfloor \frac{n}{p_j}\rfloor)\) 转移. 具体而言：

\[\begin{align*} T_k(\lfloor \frac{n}{p_j} \rfloor) &=\sum_{e=1}^{p_j^e\le\lfloor \frac{n}{p_j}\rfloor} (p_j^e)^kS_j(\lfloor \frac{\lfloor \frac{n}{p_j} \rfloor}{p_j^e}\rfloor)\\ &=\sum_{e=1}^{p_j^{e+1}\le n} (p_j^e)^kS_j(\lfloor \frac{n}{p_j^{e+1}}\rfloor)\\ &=\sum_{e=2}^{p_j^{e}\le n} (p_j^{e-1})^kS_j(\lfloor \frac{n}{p_j^{e}}\rfloor) & e\text{变为}e+1\\ \end{align*} \]

\[\begin{align*} T_k(n) &=\sum_{e=1}^{p_j^e\le n}(p_j^e)^kS_j(\lfloor \frac{n}{p_j^e}\rfloor)\\ &=p_j^kS_j(\lfloor\frac{n}{p_j}\rfloor)+\sum_{e=2}^{p_j^e\le n}(p_j^e)^kS_j(\lfloor \frac{n}{p_j^e}\rfloor) & \text{移出}e=1 \\ &=p_j^kS_j(\lfloor\frac{n}{p_j}\rfloor)+\sum_{e=2}^{p_j^e\le n}p_j^k(p_j^{e-1})^kS_j(\lfloor \frac{n}{p_j^e}\rfloor) & \text{移出} p_j^k\\ &=p_j^kS_j(\lfloor\frac{n}{p_j}\rfloor)+p_j^k\sum_{e=2}^{p_j^e\le n}(p_j^{e-1})^kS_j(\lfloor \frac{n}{p_j^e}\rfloor) & \text{移出} p_j^k\\ &=p_j^kS_j(\lfloor\frac{n}{p_j}\rfloor)+p_j^kT_k(\lfloor \frac{n}{p_j^e}\rfloor) & \text{代入}T_k\\ \end{align*} \]

当 \(p_{j+1}>\sqrt n\) 时，对 \(S_j\) 产生贡献的所有 \(F(x)\) 满足 \(x \in \mathbb{P}\).
因此，此时的

\[S_j(n)=S_\text{pri}(n)-S_\text{pri}(p_j) \]

所以我们可以从满足 \(p_j\le \sqrt{n}\) 的最大的 \(j\) 算起，通过转移式一直转移到 \(j=0\)，从而求出 \(S_0(n)\).

时间复杂度：

质数密度估计为 \(O(\frac{n}{\ln{n}})\)
由于 \(\lfloor\frac{n}{i}\rfloor\) 只有 \(O(\sqrt n)\) 个值，时间复杂度为 \(O(\sum\limits_{i=1}^{\sqrt{n}}\frac{\sqrt{i}}{\ln\sqrt{i}}+\sum\limits_{i=1}^{\sqrt{n}}\frac{\sqrt{\lfloor\frac{n}{i}\rfloor}}{\ln\sqrt{\lfloor\frac{n}{i}\rfloor}})=O(\frac{n^{\frac{3}{4}}}{\log n})\)

第二部分

第一部分中，求 \(S_j(i)\) 时用到了 \(S_{\text{pri}}(i)\) 和 \(S_{\text{pri}}(p_j)\)，其中 \(i\) 的取值为 \(n, \lfloor\frac{n}{2}\rfloor,\lfloor\frac{n}{3}\rfloor,\cdots,1\)，\(j\) 满足 \(p_j\le\sqrt{n}\).

后者可以通过线性筛筛出 \(\sqrt{n}\) 内的质数算出，而前者的计算方法如下.

首先把 \(F(\mathbb{P})\) 拆成形如 \(\mathbb{P}^k\) 的若干项，对每一项分别求解.

\(S_{\text{pri}}(i)=\sum\limits_{x=1,x\in\mathbb{P}}^{i}F(i)=\sum\limits_{x=1,x\in\mathbb{P}}^{i}x^k\)

设 \(S'_j(i)=\sum\limits_{x=1,x\in \mathbb{P}\operatorname{or}\text{lpf}(i)>p_j}^{i}x^k\)，即埃氏筛法中筛完第 \(j\) 轮后 \(1...n\) 中剩下的数的函数值之和.

同样的，我们尝试从 \(S'_{j-1}(i)\) 转移到 \(S'_j(i)\).

\[\begin{align*} S'_{j-1}(i)-S'_j(i)&=\sum_{x=1,\text{lpf}(x)=p_j\operatorname{and} x\not=p_j}^i x^k\\ &=p_j^k\sum_{x=1,\text{lpf}(x)\ge p_j}^{\lfloor\frac{i}{p_j}\rfloor} x^k\\ &=p_j^k(S'_{j-1}(\lfloor\frac{i}{p_j}\rfloor)-S_{\text{pri}}(p_{j-1})) \end{align*} \]

\(S'_0(i)=\sum\limits_{x=2}^{i}x^k\)，可以快速求得.

根据埃氏筛法可知，若 \(j\) 取满足 \(p_j\le\sqrt{i}\) 时的最大值，\(S_{\text{pri}}(i)=S'_{j}(i)\).

对于 \(i\)，发现有用的取值也是 \(n, \lfloor\frac{n}{2}\rfloor,\lfloor\frac{n}{3}\rfloor,\cdots,1\).

我们把 \(j\) 从 \(0\) 一直转移到满足 \(p_j^2\le n\) 时的最大值，即可求出 \(i\) 取 \(n, \lfloor\frac{n}{2}\rfloor,\lfloor\frac{n}{3}\rfloor,\cdots,1\) 时的所有 \(S_{\text{pri}}(i)\).

时间复杂度与优化后的第一部分是一样的.

总结

该算法分3步完成:

1. 筛出 \(\sqrt{n}\) 之内的素数并求它们的函数值的前缀和;
2. 求出 \(S_{\text{pri}}\) 数组(第二部分);
3. 求出 \(S\) 数组(第一部分).

其中第一步时间复杂度 \(O(\sqrt{n})\)，第二部和第三部均为 \(O(\frac{n^{\frac{3}{4}}}{\log n})\)，所以总时间复杂度为 \(O(\frac{n^{\frac{3}{4}}}{\log n})\).

例题

以洛谷模板题 P5325 为例.

由题知，\(F(p^e) = (p^e)^2-p^e\)，即可拆为两项处理.

代码如下:

/*
本代码用DP转移的方式计算S，常数较大，但能在求出S(n)的同时求出S(n/2),S(n/3),...,S(1)
*/
#include <bits/stdc++.h>
using namespace std;
const long long mod = 1000000007;
const int maxn = 200000;
long long n;
int rt;
int pri[maxn + 5];
bool vis[maxn + 5];
int pri_cnt;
long long prisum1[maxn + 5], prisum2[maxn + 5];
long long g1[maxn + 5], g2[maxn + 5];
long long t1[maxn + 5], t2[maxn + 5];
long long l[maxn + 5];
long long s[maxn + 5];
int id1[maxn + 5], id2[maxn + 5];
long long num[maxn + 5];
void sieve(int l) {
	for (int i = 2; i <= l; i++) {
		if (!vis[i])
			pri[++pri_cnt] = i;
		for (int j = 1; i * pri[j] <= l; j++) {
			vis[i * pri[j]] = true;
			if (i % pri[j] == 0)
				break;
		}
	}
}
long long find(long long t) {
	return t = (t <= rt) ? id1[t] : id2[n / t];
}
int cnt = 1;
void calcg() { // 计算Spri
	for (int i = 1; i <= pri_cnt; i++) {
		prisum1[i] = (prisum1[i - 1] + pri[i]) % mod;
		prisum2[i] = (prisum2[i - 1] + ((long long)pri[i] * pri[i]) % mod) % mod;
	}
	for (long long l = 1, r; l <= n; l = r + 1, cnt++) {
		r = n / (n / l);
		num[cnt] = n / l;
		long long tmp = num[cnt] % mod;
		g1[cnt] = (((tmp + 1) * tmp / 2ll) % mod + mod - 1) % mod;
		g2[cnt] = (((((tmp + 1) * tmp % mod) * (2ll * tmp + 1)) % mod) * 166666668 + mod - 1) % mod;
		if (num[cnt] <= rt)
			id1[num[cnt]] = cnt;
		else
			id2[l] = cnt;
	}
	for (int j = 1; j <= pri_cnt; j++) {
		for (int i = 1; i <= cnt && (long long)pri[j] * pri[j] <= num[i]; i++) {
			long long t = find(num[i] / pri[j]);
			g1[i] = (g1[i] - pri[j] * (g1[t] - prisum1[j - 1] + mod) % mod + mod) % mod;
			g2[i] = (g2[i] - ((long long)pri[j] * pri[j] % mod) * (g2[t] - prisum2[j - 1] + mod) % mod + mod) % mod;
		}
	}
}
void calcs() { // 计算S
	memset(s, -1, sizeof(s));
	for (int j = pri_cnt; j > 0; j--) {
		int top = 1;
		while (top <= cnt && (long long)pri[j] * pri[j] <= num[top])
			top++;
		top--;
		long long p = pri[j];
		l[top + 1] = 0;
		for (int i = top; i > 0; i--) {
			if (s[i] == -1)
				s[i] = (g2[i] - prisum2[j] - g1[i] + prisum1[j] + 2ll * mod) % mod;
			l[i] = l[i + 1];
			while (num[i] >= p) {
				l[i] = (l[i] + (long long)(p % mod) * (p % mod) % mod - p + mod) % mod;
				p *= pri[j];
			}
			long long t = find(num[i] / pri[j]);
			if (num[i] / pri[j] < (long long)pri[j] * pri[j]) {
				t1[i] = (long long)pri[j] * (g2[t] - prisum2[j] - g1[t] + prisum1[j] + 2ll * mod) % mod;
				t2[i] = ((long long)pri[j] * pri[j] % mod) * (g2[t] - prisum2[j] - g1[t] + prisum1[j] + 2ll * mod) % mod;
			} else {
				t1[i] = (long long)pri[j] * (t1[t] + s[t]) % mod;
				t2[i] = ((long long)pri[j] * pri[j] % mod) * (t2[t] + s[t]) % mod;
			}
		}
		for (int i = 1; i <= top; i++)
			s[i] = (s[i] + l[i] + t2[i] - t1[i] + mod) % mod;
	}
}
int main() {
	scanf("%lld", &n);
	rt = sqrt(n);
	sieve(rt);
	calcg();
	calcs();
	printf("%lld", max(0ll, s[find(n)]) + 1);
}

/*
本代码用递归的方式计算S，常数较小
*/
#include <bits/stdc++.h>
using namespace std;
const long long mod = 1000000007;
const int maxn = 200000;
long long n;
int rt;
int pri[maxn + 5];
bool vis[maxn + 5];
int pri_cnt;
long long prisum1[maxn + 5], prisum2[maxn + 5];
long long g1[maxn + 5], g2[maxn + 5];
int id1[maxn + 5], id2[maxn + 5];
long long num[maxn + 5];
void sieve(int l) {
	for (int i = 2; i <= l; i++) {
		if (!vis[i])
			pri[++pri_cnt] = i;
		for (int j = 1; i * pri[j] <= l; j++) {
			vis[i * pri[j]] = true;
			if (i % pri[j] == 0)
				break;
		}
	}
}
void calcg() {
	for (int i = 1; i <= pri_cnt; i++) {
		prisum1[i] = (prisum1[i - 1] + pri[i]) % mod;
		prisum2[i] = (prisum2[i - 1] + ((long long)pri[i] * pri[i]) % mod) % mod;
	}
	int cnt = 1;
	for (long long l = 1, r; l <= n; l = r + 1, cnt++) {
		r = n / (n / l);
		num[cnt] = n / l;
		long long tmp = num[cnt] % mod;
		g1[cnt] = (((tmp + 1) * tmp / 2ll) % mod + mod - 1) % mod;
		g2[cnt] = (((((tmp + 1) * tmp % mod) * (2ll * tmp + 1)) % mod) * 166666668 + mod - 1) % mod;
		if (num[cnt] <= rt)
			id1[num[cnt]] = cnt;
		else
			id2[l] = cnt;
	}
	for (int j = 1; j <= pri_cnt; j++) {
		for (int i = 1; i <= cnt && (long long)pri[j] * pri[j] <= num[i]; i++) {
			long long t = num[i] / pri[j];
			t = (t <= rt) ? id1[t] : id2[n / t];
			g1[i] = (g1[i] - pri[j] * (g1[t] - prisum1[j - 1] + mod) % mod + mod) % mod;
			g2[i] = (g2[i] - ((long long)pri[j] * pri[j] % mod) * (g2[t] - prisum2[j - 1] + mod) % mod + mod) % mod;
		}
	}
}
long long f(long long x, int y) {
	if (pri[y] >= x)
		return 0;
	int t = (x <= rt) ? id1[x] : id2[n / x];
	long long res = (g2[t] - prisum2[y] - g1[t] + prisum1[y] + 2ll * mod) % mod;
	for (int i = y + 1; i <= pri_cnt && (long long)pri[i] * pri[i] <= x; i++) {
		long long num = pri[i];
		for (int j = 1; num <= x; j++, num *= pri[i]) {
			long long tmp = num % mod;
			res = (res + (tmp * (tmp - 1 + mod) % mod) * (f(x / num, i) + (int)(j > 1)) % mod) % mod;
		}
	}
	return res;
}
int main() {
	scanf("%lld", &n);
	rt = sqrt(n);
	sieve(rt);
	calcg();
	printf("%lld", f(n, 0) + 1);
}