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\).
- 该函数为积性函数（个别非积性函数也可以）;
- 若 \(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)\)
该式可以递归求解(下文 P5325 代码 1)，也可以通过 \(S_{j}(i)\) 转移出 \(S_{j-1}(i)\) (P5325 代码 2)
式子的第一部分( \(\sum\limits_{e=1}^{p_j^e\le i}F(p_j^e)\) )可以快速地求和（暴力枚举 \(e\)）.
现在考虑第二部分：
令
由于对于 \(\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_{j,k}\)，就可以求和得到 \(T_j\)）
对于 \(T_k\) 考虑从 \(T_k(\lfloor \frac{n}{p_j}\rfloor)\) 转移. 具体而言：
当 \(p_{j+1}>\sqrt n\) 时，对 \(S_j\) 产生贡献的所有 \(F(x)\) 满足 \(x \in \mathbb{P}\).
因此，此时的
所以我们可以从满足 \(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)\).
\(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步完成:
- 筛出 \(\sqrt{n}\) 之内的素数并求它们的函数值的前缀和;
- 求出 \(S_{\text{pri}}\) 数组(第二部分);
- 求出 \(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);
}
