Min_25筛

Min_25 筛

tags: 数学 数论

---

简介

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

适用条件

设要求的函数为 (F).

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

记号

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

(S_j(n) = sumlimits_{i = 1, ext{lpf}(i)>p_j}^{n}F(i))

(S_{ ext{pri}}(n)=sumlimits_{i=2,iin mathbb{P}}^{n}F(i))

求解方法

第一部分

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

[egin{align*} S_j(i) &= sumlimits_{x = 1, ext{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, ext{lpf}(y)>p_x}^{lfloor frac{i}{p_i} floor}F(y)) & ext{提出最小质因子} \ &= 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, ext{lpf}(y)>p_x}^{lfloor frac{i}{p_x} floor}F(y))+sum_{i=p_j+1,iin mathbb{P}}^{i}F(y) & ext{提出质数} \ &= 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, ext{lpf}(y)>p_x}^{lfloor frac{i}{p_x} floor}F(y))+sum_{i=p_j+1,iin mathbb{P}}^{i}F(y) & p_x ext{上界变为}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, ext{lpf}(y)>p_x}^{lfloor frac{i}{p_x} floor}F(y))+(S_{mathtt{pri}}(i)-S_{mathtt{pri}}(p_j)) & ext{代入}S_ ext{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} floor))+(S_{mathtt{pri}}(i)-S_{mathtt{pri}}(p_j)) & ext{代入}S_x end{align*} ]

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

[egin{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} floor))+(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} floor))+(S_{mathtt{pri}}(i)-S_{mathtt{pri}}(p_j))\ &=sum_{e=1}^{p_j^ele i}F(p_j^e)([e>1]+S_j(lfloor frac{i}{p_j^e} floor))+(S_{ ext{pri}}(p_j)-S_{ ext{pri}}(p_{j-1})) & ext{第一层}sum ext{只剩下}p_j \ &=sum_{e=1}^{p_j^ele i}F(p_j^e)([e>1]+S_j(lfloor frac{i}{p_j^e} floor))+F(p_j) & ext{消去}S_ ext{pri}\ &=sum_{e=2}^{p_j^ele i}F(p_j^e)+sum_{e=1}^{p_j^ele i}F(p_j^e)S_j(lfloor frac{i}{p_j^e} floor)+F(p_j) & ext{拆开括号} \ &=sum_{e=1}^{p_j^ele i}F(p_j^e)+sum_{e=1}^{p_j^ele i}F(p_j^e)S_j(lfloor frac{i}{p_j^e} floor) & ext{合并} end{align*} ]

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

现在考虑第二部分：
令

[T(n)=sum_{e=1}^{p_j^ele n}F(p_j^e)S_j(lfloor frac{n}{p_j^e} floor) ]

由于对于 (forall x=p^e,pinmathbb{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^ele i}(p_j^e)^kS_j(lfloor frac{i}{p_j^e} floor) ]

有

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

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

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

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

[egin{align*} T_k(n) &=sum_{e=1}^{p_j^ele n}(p_j^e)^kS_j(lfloor frac{n}{p_j^e} floor)\ &=p_j^kS_j(lfloorfrac{n}{p_j} floor)+sum_{e=2}^{p_j^ele n}(p_j^e)^kS_j(lfloor frac{n}{p_j^e} floor) & ext{移出}e=1 \ &=p_j^kS_j(lfloorfrac{n}{p_j} floor)+sum_{e=2}^{p_j^ele n}p_j^k(p_j^{e-1})^kS_j(lfloor frac{n}{p_j^e} floor) & ext{移出} p_j^k\ &=p_j^kS_j(lfloorfrac{n}{p_j} floor)+p_j^ksum_{e=2}^{p_j^ele n}(p_j^{e-1})^kS_j(lfloor frac{n}{p_j^e} floor) & ext{移出} p_j^k\ &=p_j^kS_j(lfloorfrac{n}{p_j} floor)+p_j^kT_k(lfloor frac{n}{p_j^e} floor) & ext{代入}T_k\ end{align*} ]

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

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

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

时间复杂度：

质数密度估计为 (O(frac{n}{ln{n}}))
由于 (lfloorfrac{n}{i} floor) 只有 (O(sqrt n)) 个值，时间复杂度为 (O(sumlimits_{i=1}^{sqrt{n}}frac{sqrt{i}}{lnsqrt{i}}+sumlimits_{i=1}^{sqrt{n}}frac{sqrt{lfloorfrac{n}{i} floor}}{lnsqrt{lfloorfrac{n}{i} floor}})=O(frac{n^{frac{3}{4}}}{log n}))

第二部分

第一部分中，求 (S_j(i)) 时用到了 (S_{ ext{pri}}(i)) 和 (S_{ ext{pri}}(p_j))，其中 (i) 的取值为 (n, lfloorfrac{n}{2} floor,lfloorfrac{n}{3} floor,cdots,1)，(j) 满足 (p_jlesqrt{n}).

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

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

(S_{ ext{pri}}(i)=sumlimits_{x=1,xinmathbb{P}}^{i}F(i)=sumlimits_{x=1,xinmathbb{P}}^{i}x^k)

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

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

[egin{align*} S'_{j-1}(i)-S'_j(i)&=sum_{x=1, ext{lpf}(x)=p_joperatorname{and} x ot=p_j}^i x^k\ &=p_j^ksum_{x=1, ext{lpf}(x)ge p_j}^{lfloorfrac{i}{p_j} floor} x^k\ &=p_j^k(S'_{j-1}(lfloorfrac{i}{p_j} floor)-S_{ ext{pri}}(p_{j-1})) end{align*} ]

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

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

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

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

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

总结

该算法分3步完成:

1. 筛出 (sqrt{n}) 之内的素数并求它们的函数值的前缀和;
2. 求出 (S_{ ext{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);
}