Min_25 筛
简介
名字来源于其发明者的 \(ID\)
又名 \(\;Extended\;Eratosthenes\;Sieve\)（扩展欧拉筛）
设 \(f(n)\) 为一个积性函数，若可以 \(O(1)\) 求出 \(f(p^r)\)（\(p\) 取值为全体质数，下同；\(r\) 为正整数），且 \(f(p)\) 为一个关于 \(p\) 的多项式，则使用 \(Min\_25\) 筛可以以 \(O\left(\dfrac{n^{\tfrac34}}{\ln n}\right)\) 或 \(O(n^{1-\epsilon})\) 的时间复杂度、\(O(\sqrt n)\) 的空间复杂度求出 \(\sum_{i=1}^nf(i)\)
由于不需要构造，且时间复杂度（也考虑了常数）优于杜教筛，因此更加常用
符号
- \(a/b\)：\(\lfloor\frac ab\rfloor\)，即 \(a\) 除以 \(b\) 下取整
- \(isp(x)\)：\([(\exists d\ne 1,d\ne x),d\mid x]\)，即当 \(x\) 为质数时 \(isp(x)=1\)，否则 \(isp(x)=0\)
- \(p_x\)：第 \(x\) 个质数。特别地，\(p_0=1\)
- \(lpf(x)\)：\(x\) 的最小质因子（\(least\;prime\;factor\)）。特别地，\(lpf(1)=1\)
- \(F_p(n)\)：\(\sum_{i=2}^n[isp(i)]f(i)\)，即 \(2\sim n\) 中仅取质数的 \(x\) 的 \(f(x)\) 和
- \(F_k(n)\)（\(k\) 为数字，区别于 \(F_p(n)\)）：\(\sum_{i=2}^n[lpf(i)>p_k]f(i)\)，即 \(2\sim n\) 中取 最小质因子大于 \(p_k\) 的 \(x\) 的 \(f(x)\) 和
过程
根据定义有 \(F_0(n)=\sum_{i=2}^nf(i)\)。由于 \(f\) 为积性函数，因此 \(f(1)=1\)，故答案为 \(\sum_{i=1}^nf(i)=F_1(n)+f(1)=F_1(n)+1\)
求解 \(F_k(m)\)
考虑快速求解一个 \(F_k(m)\)
根据定义有
先拆出 \(i\) 为质数的项，这部分的总和为
剩下的 \(i\) 都是合数（根据定义），又可以分为两部分，只有一个质因数的（此时其次数一定大于 \(1\)）和有两个以上质因数的
对于前者，枚举这个质因子 \(p_i\) 和其次数 \(c\)，总和为
对于后者，类似地枚举最小质因子 \(p_i\) 和其次数 \(c\)。由于 \(f\) 是积性函数，其贡献可以拆成 \(f(p_i)^c\) 与剩余部分的乘积，总和为
（式中 \(p_i^2\le m\) 是因为合数的最小质因子的平方一定不大于其本身）
两者可以合并为
再加上前面的质数部分，可得
边界条件为 \(p_k\ge n\)，此时显然 \(F_k(n)=0\)
将 \(F_k(m)\) 转化为求解若干个 \(F_i(m/v)\;(i>k)\) 和两个 \(F_p(v)\)
由于求解原问题只要 \(F_0(n)\)，每次递归出的 \(m'\) 都是除以某数下取整，因此整个计算过程中用到的 \(F_k(m)\) 中，\(m\) 的取值只有 \(1,2,\cdots,\lfloor\sqrt n\rfloor,n/\lfloor\sqrt n\rfloor,n/(\lfloor\sqrt n\rfloor-1),\cdots,n/1\) 共 \(O(\sqrt n)\) 种
同时由其递归过程可以看出，实际进行计算的 \(F_k(m)\) 中，\(p_k\le m\) 且 \(p_k^2\le n\)，即 \(k\le \pi(\sqrt n)=O\left(\dfrac{\sqrt n}{\log n}\right)\)（更紧的上界为将 \(n\) 替换为调用它的 \(F_k(m)\)的 \(m\)）
因此 \(F_k(m)\) 递归求解时，用到的两个 \(F_p\) 中，前一个的参数取值有 \(1,2,\cdots,\lfloor\sqrt n\rfloor,n/\lfloor\sqrt n\rfloor,n/(\lfloor\sqrt n\rfloor-1),\cdots,n/1\) 共 \(O(\sqrt n)\) 种，后一个的参数不超过 \(\sqrt n\)
若已经处理出了所有需要的 \(F_p\) 且可以 \(O(1)\) 查询，则可以证明求解 \(F_k(n)\) 的时间复杂度为 \(O\left(\dfrac{n^{\tfrac34}}{\ln n}\right)\)（一说 \(O(n^{1-\epsilon})\)）
问题转化为在 \(O\left(\dfrac{n^{\tfrac34}}{\ln n}\right)\) 内求出所有 \(F_p(n/x)\)
求解 \(F_p(m)\)
考虑如何快速计算出一个 \(F_p(m)\)
由于保证了 \(f(p)\) 为关于 \(p\) 的多项式（设 \(f_p=\sum_{i=0}^wa_ip^i\)，其中 \(w\) 和 \(a_i\) 为常数），因此
因此问题转化为对于常数 \(v\)，求出 \(\large{\sum_{i\ge 1\land p_i\le m}p_i^v}\)，其中 \(m\in\{N/x\}\)
令 \(g(x)=x^v\)（显然为完全积性函数），\(G_k(m)=\sum_{i=2}^m[lpf(i)>p_k\lor isp(i)]g(i)\)，即 \(2\sim m\) 中的数，在第 \(k\) 轮（标准埃氏筛，即只在扫到质数时才筛一轮）埃氏筛结束后剩下数的 \(g\) 之和
显然第 \(pi(\lfloor\sqrt m\rfloor)\) 轮筛之后剩下的一定都是质数，不会再变化，因此 \(\large{F_p(m)=G_{pi(\lfloor\sqrt m\rfloor)}(m)}\)
问题转化为对于所有 \(m\in\{n/x\}\) 和所有 \(g(x)\)，求出 \(\large{G_{pi(\lfloor\sqrt m\rfloor)}(m)}\)
求解 \(G_k(m)\)
考虑计算 \(G_k(m)\)
当 \(k=0\) 时，\(G_k(m)\) 等于 \(2\sim m\) 中数的 \(g\) 之和，即 \(\sum_{i=2}^m g(i)=\sum_{i=2}^m i^v\)。一般 \(g\) 的次数在 \(0,1,2,3\) 之中，很少有更高次方的，因此几个常见求和公式就够用了
当 \(k>0\) 时，从 \(G_{k-1}(m)\) 转移过来，需要减去恰好在第 \(k\) 轮被删的 \(g\) 的总和
显然 \(p_k^2>m\) 时，第 \(k\) 轮不会删去任何数，即此时 \(G_k(m)=G_{k-1}(m)\)
否则在第 \(k\) 轮删去的是 \(lpf\) 等于 \(p_k\) 的合数
显然它们都含有质因子 \(p_k\)
由于 \(g\) 为完全积性函数，因此删去部分的 \(g\) 之和是 \(g(p_k)\) 乘以 满足 \(lpf(x)\ge p_k,x\le m/p_k\) 的 \(g(x)\) 之和
而 满足 \(lpf(x)\ge p_k,x\le m/p_k\) 的 \(g(x)\) 之和 等于 满足 \(lpf(x)\ge p_{k-1}\lor isp(x),x\le m/p_k\) 的 \(g(x)\) 之和 减去 满足 \(lpf(x)\ge p_{k-1}\lor isp(x),x\le p_{k-1}\) 的 \(g(x)\) 之和
前者等于 \(G_{k-1}(m/p_k)\)，后者等于 \(G_{k-1}(p_{k-1})\)
因此 \(G_k(m)=G_{k-1}(m)-g(p_k)(G_{k-1}(m/p_k)-G_{k-1}(p_{k-1}))\)
由于 \(G_k\) 的参数都是 \(N/x\)，一共 \(O(\sqrt n)\) 种，可以将同一个 \(g\) 的 \(\le \sqrt n\) 的和 \(>\sqrt n\) 的 \(G\) 通过映射存储到一起
实现时，可以对于每个需要的 \(g(x)=x^v\)，建立一个长为 \(2\sqrt n\) 的数组 \(g_v\)（将上文需要的 \(2\sqrt n\) 个位置映射到下标 \(1\sim2\sqrt n\)）和长为 \(\pi(\sqrt n)\) 的数组 \(gp_v\)
其中 \(\large{gp_v(x)=\sum_{i=1}^{p_x}[isp(i)]i^v}\)，且 \(x\le \pi(\sqrt n)\)
而 \(g_v(x)\) 初始存储上文中的 \(G_0(x)\)，通过 \(\pi(\sqrt x)\) 次递推到 \(G_{\pi(\sqrt x)}(x)\)
总计 \(\pi(\sqrt n)\) 次递推，每次递推只改变 \(p_k^2\le x\) 的 \(G_k(x)\) 位置，可以保证复杂度
总流程
- 求出 \(\sqrt N\) 以内的质数列表和 \(gp\) 的值
- 通过递推求出 \(G\)
- 根据 \(G\) 求出 \(F_p\) 在 \(\sqrt n\) 个位置处的值
- 根据 \(F_p(x)\) 的值，递归计算出 \(F_1(n)\)，得出答案
可以证明，时间复杂度为 \(O\left(\dfrac{n^{\frac34}}{\ln n}\right)\)
模板题：P5325 【模板】Min_25 筛 \(\quad\) 代码
\(Min\_25\) 筛模板：
#include <bits/stdc++.h>
#define M 1000000007
using namespace std;
namespace Min_25{
constexpr long long N_max = 1e10;
constexpr int sqN_mx = 1e5 + 10;
constexpr int _2sqN_mx = 2 * sqN_mx;
int add(int a, int b){a += b; return a >= M? a - M : a;}int sub(int a, int b){a -= b; return a < 0? a + M : a;}int mul(int a, int b){return 1ll * a * b % M;}
void sadd(int &x, int y){x += y;if (x >= M)x -= M;}void ssub(int &x, int y){x -= y;if (x < 0)x += M;}void dec(int &v){if (v == 0)v = M; --v;}
int fpw(int a, int b){int ret = 1;for (; b; b >>= 1, a = mul(a, a))if (b & 1)ret = mul(ret, a);return ret;}
long long n; int sqn;int prm[sqN_mx], L;
int /*gp0[_2sqN_mx],*/gp1[_2sqN_mx], gp2[_2sqN_mx]/*, gp3[_2sqN_mx], ...*/;//g on prime,^0,^1,^2,^3,...
int /*g0[_2sqN_mx], */g1[_2sqN_mx], g2[_2sqN_mx]/*, g3[_2sqN_mx], ...*/;//g,^0,^1,^2,^3,...
int id[_2sqN_mx];int &gid(long long v){return id[v <= sqn? v : n / v + sqn];}
/*--------------------------- need to change ---------------------------*/
/*func in this problem = p^k(p^k-1), so p^2-p^1 on primes*/
/**/int f(int, int, int pr){return mul(pr, sub(pr, 1));}//args are p,r,p^r, to calc f O(1), need to change on different function
/**/int Fp(long long v){int I = gid(v);return sub(g2[I], g1[I]);}//
/**/int Fpp(int k){return sub(gp2[k], gp1[k]);}//Fp(p_{k})
/*----------------------------------------------------------------------*/
void solprm(int n){
vector<bool> tgt(n + 1);
tgt[1] = 1;
for (int i = 1; i <= n; ++i){
if (!tgt[i]){
prm[++L] = i;
/*gp0[L] = add(gp0[L - 1], 1);*/
gp1[L] = add(gp1[L - 1], i);
gp2[L] = add(gp2[L - 1], mul(i, i));
/*gp3[L] = add(gp3[L - 1], mul(mul(i, i), i)); ...*/
}
for (int j = 1; j <= L && 1ll * prm[j] * i <= n; ++j){tgt[prm[j] * i] = 1;if (i % prm[j] == 0)break;}
}
}
int F(int k, long long n){
if (prm[k] >= n)return 0;int ret = sub(Fp(n), Fpp(k));
for (int i = k + 1; i <= L && 1ll * prm[i] * prm[i] <= n; ++i)
{long long Tp = prm[i];for (int c = 1; Tp <= n; ++c, Tp *= prm[i])ret = add(ret, mul(f(prm[i], c, Tp % M), F(i, n / Tp) + (c > 1)));}
return ret;
}
long long w[_2sqN_mx]; int tot;
int sol(long long N){
const int IV2 = fpw(2, M - 2);
const int IV6 = fpw(6, M - 2);
n = N;sqn = sqrtl(n);
solprm(sqn);
for (long long l = 1, r; l <= n; l = r + 1){
r = n / (n / l);w[++tot] = n / l;int v = w[tot] % M;gid(w[tot]) = tot;
/*g0[tot] = v;*/
g1[tot] = mul(mul(v, add(v, 1)), IV2);
g2[tot] = mul(mul(v, add(v, 1)), mul(add(mul(v, 2), 1), IV6));
/*g3[tot] = mul(g1[tot], g1[tot]); ...*/
/*dec(g0[tot]), */dec(g1[tot]), dec(g2[tot])/*, dec(g3[tot])*/;//remove 1
}//calc G(0,..)
for (int i = 1; i <= L; ++i){//G(i,..)
for (int k = 1; k <= tot && 1ll * prm[i] * prm[i] <= w[k]; ++k){
int id = gid(w[k] / prm[i]);
/*ssub(g0[k], mul(1, sub(g0[id], gp0[i - 1])));*/
ssub(g1[k], mul(prm[i], sub(g1[id], gp1[i - 1])));
ssub(g2[k], mul(mul(prm[i], prm[i]), sub(g2[id], gp2[i - 1])));
/*ssub(g3[k], mul(mul(mul(prm[i], prm[i]), prm[i]), sub(g3[id], gp3[i - 1]))); ...*/
}
}
return add(F(0, n), 1);//+f(1)
}
}
int main(){
ios::sync_with_stdio(false), cin.tie(nullptr), cout.tie(nullptr);
long long n;
cin >> n;
cout << Min_25::sol(n) << endl;
return 0;
}//refer https://www.luogu.com.cn/article/v5oxc6v1
/*
10 263
1000000000 710164413
*/
浙公网安备 33010602011771号