标签：25 frac 看不懂 Min 积性 sum 合数 minp 质数

前言

由于其由 Min_25 发明并最早开始使用，故称「Min_25 筛」。

• 从此种筛法的思想方法来说，其又被称为「Extended Eratosthenes Sieve」。

其可以在 \(O(\frac{n^{\frac{3}{4}}}{logn})\) 的时间复杂度下解决一类 积性函数 的前缀和问题。
要求： \(f(p)\) 是一个关于 \(p\) 的项数较少的多项式或可以快速求值； \(f(p^{c})\) 可以快速求值。
摘自 OI wiki

解析

例题

积性函数 \(f(p^k) = p^k(p^k - 1)\)
求：\(\sum_{i = 1} ^ n f(i)\)

准备 \(1\)

规定 \(p\) 代表质数, \(minp(i)\) 代表 \(i\) 的最小质因子（其中 \(minp(p) = p\)），则：
\(\sum_{i = 1} ^ n f(i) = \sum_{p \leq n} f(p) + \sum_{i \leq n \ and \ i \notin p} f(i)\)
这一步表示将答案拆成质数的答案和合数的答案。
那么对于合数部分：
根据积性函数的性质，则有：
\(f(i) = \prod_{k = 1} ^ m f(p_k ^ {c_k})\)。
提尽 \(i\) 的最小质因子 \(p_e\)，将 \(i\) 分成两部分：
\(p_e ^ {c_e} \times \frac{i}{p_e ^ {c_e}}\)
则上式可变为：
\(f(i) = f(p_e ^ {c_e}) \times f(\frac{i}{p_e ^ {c_e})})\)
那么合数部分的答案可转化为：
\(\sum_{p ^ e \leq n} f(p ^ e) \sum_{1 \leq i \leq n \ and \ p ^ e \mid i \ and \ minp(\frac{i}{p ^ e}) > p} f(\frac{i}{p ^ e})\)
代表的意思是：
枚举质数 \(p\) 及它的指数 \(e\)，
在可以提的合数 \(i\) 中提出 \(p ^ e\)，将合数 \(i\) 转化为 \(p ^ e \times \frac{i}{p ^ e}\)，即上式乘法分配律展开的结果。
然后再判断如下条件：

• \(p\) 是否是 \(i\) 的最小质因子
• \(p\) 是否提尽

这两个条件均可通过 \(minp(\frac{i}{p ^ e})\) 是否大于 \(p\) 来判断：
对于 \(minp(\frac{i}{p ^ e}) < p\) 的 \(i\)，那么不满足条件 \(1\)。
对于 \(minp(\frac{i}{p ^ e}) = p\) 的 \(i\)，那么不满足条件 \(2\)。
这些均可通过对 \(i\) 分解质因数来理解。
那么对于 \(f(p ^ e)\) 由题意可以直接求，对于 \(f(\frac{i}{p ^ e})\) 则是同样的求法，直到 \(\frac{i}{p ^ e}\) 是 \(p_k ^ {c_k}\) 可以直接求。
而上述的约束则保证了不会多求漏求。

准备 \(2\)

所以我们现在需要求出所有 \(f(p ^ k)\) 的值。
因为 \(n\) 的范围很大，所以我们没法线性筛求出 \(f\) 的前缀和。
考虑 DP。。。
设 \(h_k(i) = i ^ k\)，显然 \(h_k\) 是一个完全积性函数。
设 \(g_k(n,j) = \sum_{i = 1} ^ n [i \in Prime \ or \ minp(i) > p_j] h_k(i)\)，其中 \(p_j\) 表示第 \(j\) 个质数。
换句话说 \(g_k(n,j)\) 也就是 \(1\) ~ \(n\) 中质数或者最小质因数大于 \(p_j\) 的 \(i\) 的 \(h_k(i)\) 的前缀和。
\(g_k\) 的作用在下文，但也可以先抛开问题本身，看懂 \(g_k\) 的求法。
考虑转移。。。
对于 \(g_k(n, j - 1)\) 根据定义，其中包含两部分：

• \(g_k(n, j)\)
• 最小质因数等于 \(p_j\) 的 \(h_k(i)\) 的和。

所以可由 \(g_k(n, j - 1)\) 减去多余的部分得到 \(g_k(n, j)\)。
所以有状态转移方程如下：
\(g_k(n,j) = g_k(n, j - 1) - p_{j} ^ k (g(\frac{n}{p_{j} ^ k}, j - 1) - g(p_j, j - 1))\)
从状态转移方程可以看出：
最小质因数等于 \(p_j\) 的 \(h_k(i)\) 的和即为 \(g_k(n, j - 1) - p_{j} ^ k (g(\frac{n}{p_j}, j - 1) - g(p_{j - 1}, j - 1))\)。
原因如下：
根据定义不难得到 \(g_k(p_j, j - 1)\) 即为前 \(j - 1\) 个质数的 \(k\) 次方和。
而对于 \(g_k(\frac{n}{p_j}, j - 1)\)，包含了两部分：

• \(1\) ~ \(p_{j - 1}\) 的质数和
• 最小质因数大于等于 \(p_j\) 的合数和，及 \(p_{j}\) ~ \(\frac{n}{p_j}\) 的质数和

注意第二点，如果我们将这些最小质因数大于等于 \(p_j\) 的合数或者后一部分的质数同时乘上 \(p_j\) 那么我们得到的是不是都是最小质因数等于 \(p_j\) 的合数，这点应该很好理解。
那也就是说得到了我们想要筛去的合数。
因为我们只知道上述两部分的和 \(g_k(\frac{n}{p_j}, j - 1)\)，所以在乘的过程中第一部分也会被乘上一个 \(p_j\)，那么我们应该把这一部分质数减去。
也就是 \(g_k(p_j, j - 1)\)。
又因为 \(h_k\) 是完全积性函数，而我们谈论的和是 \(h_k\) 的累加，根据乘法分配律，我们可以对 \(g_k\) 运用完全积性函数的性质（\(f(ab) = f(a)f(b)\)）。
那么要得到最小质因数等于 \(p_j\) 的 \(h_k(i)\) 的和，我们可以直接给我们的第二部分乘上 \(h_k(p_j)\) 的值，也就是 \(p_j ^ k\)，得到我们要求的函数值。
所以关于状态转移方程的解释就结束了。
关于初始化：
显然对于 \(g_k(n, 0) = \sum_{i = 1} i ^ k\)

求解

而我们得到的 \(g_k\) 有什么用呢？
设 \(g(n, j) = \sum_{i = 1} ^ n [i \in Prime \ or \ minp(i) > p_j] f(i)\)。
因为 \(f(p) = p(p - 1) = p ^ 2 - p\)。
那么将满足条件的 \(f(i)\) 累加，就得到两部分：

• \(i ^ 2\) 的和
• \(i\) 的和的相反数

如果令上文的 \(k\) 分别等于 \(1, 2\)。
那么就有： \(g(n, j) = g_2(n, j) - g_1(n, j)\)。
如果设 \(p_x\) 为小于等于 \(\sqrt{n}\) 的最大的质数，那么 \(1\) ~ \(n\) 的质数的 \(f\) 的和 \(g(n, x)\) 记为 \(g(n)\)。
设 \(S(n, j) = \sum_{minp(i) > p_j} f(i)\)。
显然我们也需要 DP 。。。
同样的，考虑 \(S(n, j)\) 的组成：

• 质数部分：大于 \(p_j\) 小于等于 \(n\) 的质数
• 合数部分

设 \(sum_k(x)\) 表示小于等于 \(p_x\) 的质数的 \(k\) 次方的和。
那么质数部分就等于：
\(g(n) - （sum_2(j) - sum_1(j)）\)
对于合数部分：
我们模仿 \(g_k\) 的求法
因为 \(f\) 不是完全积性函数，而是积性函数，我们不能像求 \(g_k\) 那样只提一个 \(p_j\)。
所以我们要把 \(p_j\) 提尽。
得到：
\(\sum_{k > j \ and \ p_k ^ e \leq n} f(p_k ^ e)(S(\frac{n}{p_k ^ e} + [e \neq 1]))\)
利用 \(S(\frac{n}{p_k ^ e}\) 的思想和求 \(g_k\) 的时候完全一样，这里就不在赘述了，实在不懂的同学可以手推一下。
后面加上一个 \([e \neq 1]\) 是因为如果 \(e > 1\)，那么\(p ^ e\) 在我们提 \(p_k ^ e\) 时没有算进去，而 \(e = 1\) 时，在质数部分就算了，就不用了算进去。

下标的离散化

因为 \(n = 1e10\)，直接开肯定不行。
由性质 \(\lfloor \frac{\lfloor \frac{x}{a} \rfloor}{b} \rfloor = \lfloor \frac{x}{ab} \rfloor\)，在递归的时候，无论我们把 \(n\) 除以几，得到的都是 \(n\) 除以某个数，所以我们可以直接只保存 \(n / x\) 的值。
来自 wucstdio 大佬

代码

#include<cstdio>
#include<cmath>
#include<cstring>

using namespace std;

typedef long long LL;

const LL mod = 1e9 + 7, inv2 = 500000004, inv3 = 333333336;
const int N = 1e5 + 5;

LL n, Prime[N + 5], num, sp1[N + 5], sp2[N + 5], tot =  0, w[3 * N], g1[3 * N], g2[3 * N], sqr;
int id1[N + 5], id2[N + 5];
bool notPrime[N + 5];

void init() {
	memset(notPrime, false, sizeof(notPrime));
	num = 0;
	for(int i = 2; i <= N; i++) {
		if(!notPrime[i]) {
			Prime[++num] = i;
			sp1[num] = (sp1[num - 1] + i) % mod;
			sp2[num] = (sp2[num - 1] + 1ll * i * i % mod) % mod;
		}
		for(int j = 1; j <= num && Prime[j] * i <= N; j++) {
			notPrime[i * Prime[j]] = true;
			if(i % Prime[j] == 0) break;
		}
	}
}

LL s(LL x, int y) {
	if(Prime[y] >= x) return 0;
	LL k = x <= sqr ? id1[x] : id2[n / x];
	LL ans = (g2[k] - g1[k] + mod - (sp2[y] - sp1[y]) + mod) % mod;
	for(int i = y + 1; i <= num && Prime[i] * Prime[i] <= x; i++) {
		LL P = Prime[i];
		for(int e = 1; P <= x; e++, P = P * Prime[i]) {
			LL xx = P % mod;
			ans = (ans + xx * (xx - 1) % mod * (s(x / P, i) + (e != 1))) % mod;
		}
	}
	return ans % mod;
}

int main() {
	init();
	scanf("%lld", &n);
	sqr = sqrt(n);
	for(LL l = 1, r; l <= n; l = r + 1) {
		r = (n / (n / l));
		w[++tot] = n / l;
		g1[tot] = w[tot] % mod;
		g2[tot] = g1[tot] * (g1[tot] + 1) / 2 % mod * (2 * g1[tot] + 1) % mod * inv3 % mod - 1;
		if(g2[tot] < 0) g2[tot] += mod;
		g1[tot] = g1[tot] * (g1[tot] + 1) / 2 % mod - 1;
		if(g1[tot] < 0) g1[tot] += mod;
		if(w[tot] <= sqr) id1[w[tot]] = tot;
			else id2[n / w[tot]] = tot;
	}
	for(int i = 1; i <= num; i++)
		for(int j = 1; j <= tot && Prime[i] * Prime[i] <= w[j]; j++) {
			LL k = w[j] / Prime[i];
			k = k <= sqr ? id1[k] : id2[n / k];
			g1[j] -= Prime[i] * (g1[k] - sp1[i - 1] + mod) % mod;
			g2[j] -= Prime[i] * Prime[i] % mod * (g2[k] - sp2[i - 1] + mod) % mod;
			g1[j] %= mod, g2[j] %= mod;
			if(g1[j] < 0) g1[j] += mod;
			if(g2[j] < 0) g2[j] += mod;
		}
	printf("%lld", (s(n, 0) + 1) % mod); 
	return 0;
}

标签：25,frac,看不懂,Min,积性,sum,合数,minp,质数
来源： https://www.cnblogs.com/sjzyh/p/15000939.html

本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享；
2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关；
3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关；
4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除；
5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。