min_25 筛介绍
我们考虑这样一个问题。
其中 (1 le n le 10^{10})
其中 (f(i)) 是一个奇怪的函数、并不像 (μ(i),φ(i),iφ(i))那样具有那么好的性质。但是满足以下条件：
- 若 (p)为质数，则 (f(p))是一个关于 (p)的多项式，比如 (μ(p)=−1,φ(p)=p−1).
- 若 (p)为质数，(e)为正整数，则 (f(pe))可以很快求出。（通常是 (O(1))）
- (f(n))为积性函数。
什么是积性函数：对于所有互质的 (a) 和 (b) ，总有 (f(ab)=f(a)f(b)) ，则称 (f(x)) 为积性函数。
然后就可以使用 min_25 筛了。（顾名思义是 min_25 发明的）
首先，我们需要知道 min_25 筛是埃氏筛的一个拓展，它的思想很大一部分借助于埃氏筛。
回想一下埃氏筛，我们是每次将最小质因子为 (P_i) 的合数筛去，剩下的就是质数。
我们知道这些最小质因子至多为 (sqrt{n})，所以合数可以通过枚举最小质因子来计算，质数我们则使用另外的方法。
首先我们看质数的怎么做。
根据条件 (1)，我们知道 (f(i))是一个多项式，这样的话我们可以按照次数将 (f(i))拆成 (i^k) 之和，因为 (i^k) 是一个完全积性函数（很快就有用的）。
为了计算这个，我们需要引入一个辅助数组 (g(n,j))。（鬼知道是怎么想到的）
其中 (minp(i))表示 (i)的最小质因子，所以
既然我们要使用质数，所以我们可以先用欧拉筛把所有 (le sqrt{n}) 的质数筛出来，同时还要预处理 (∑^j_{i=1}[iin Prime]i^k).
我们考虑 dp 计算。既然是埃氏筛，我们就要在 (g(n,j−1))中最小质因子为 (P_j)的合数筛去。
我们假设 (P_j^2 le n)，否则肯定不行。
首先，由于它是完全积性函数，所以 (P_j) 可以直接提出来，剩下的减去 (i≤⌊frac{n}{P_j}⌋)中的数就可以了。
这些数中要求质因子 (≥P_j)，所以是 (g(⌊frac{n}{P_j}⌋,j−1))，但是这里面质数被重复计算了，所以要减去里面的质数。
但是这样的话是 (O(n∗|P|))的，时间和空间都承受不了。但是我们发现我们可以使用一个优化。
我们发现 (g(n,j))中 (n)只有 (O(sqrt{n}))种取值，因为每次递归的时候是 (n)变为 (⌊frac{n}{P_j}⌋)，而我们发现
所以 (n)只会变为 (⌊nx⌋)，于是我们就直接 “手动” 离散化，这个可以看代码。
然后 (g(n,j))的第二维也可以滚动数组滚掉。所以时间 (O(sqrt{n}∗|P|))，空间 (sqrt{n}) .
预处理部分终于结束了，接下来我们考虑计算答案，首先我们还是需要一个辅助数组。
像上面说的一样，分质数和合数两类计算。
前面两项指的是质数的部分，后面的和式是枚举合数的最小质因子 (P_k)和它的次数 (e)，
这个跟 (g)不同，(S)是要按照第二维倒着计算的，但是我们也可以使用递归的方法来计算。
(S(n,0))就是最终答案。
要注意的是，(1)既不是质数也不是合数，所以最后要加上。
至于上面说的那个手动离散化，我们要开两个数组 (id1)和 (d2)，分别记录 (≤ sqrt{n}) 和 (>sqrt{n})的部分的数值的编号。这样就不用 map 了，可以省掉一个 log。
至于时间复杂度？我也不知道，总之跟杜教筛差不多，甚至有时候比杜教筛还要快。
有人说是什么 (Oleft(frac{n^{frac{3}{4}}}{log{n}} ight)) ，或者什么 (Oleft(n^{1 - epsilon} ight)) 。总之差不多就可以了。
min_25 筛代码实现：
#include<bits/stdc++.h>
#define Rint register LL
using namespace std;
typedef long long LL;
const int N = 1000003, mod = 1e9 + 7, inv2 = 500000004, inv3 = 333333336;
LL n, Sqr, pri[N], tot, pre1[N], pre2[N], ind1[N], ind2[N], g1[N], g2[N], w[N], cnt;
bool notp[N];
inline void init(LL m){
notp[0] = notp[1] = true;
for(Rint i = 2;i <= m;i ++){
if(!notp[i]){
pri[++ tot] = i;
pre1[tot] = (pre1[tot - 1] + i) % mod;
pre2[tot] = (pre2[tot - 1] + i * i) % mod;
}
for(Rint j = 1;j <= tot && i * pri[j] <= m;j ++){
notp[i * pri[j]] = true;
if(!(i % pri[j])) break;
}
}
}
inline LL S(LL x, int y){
if(pri[y] >= x) return 0;
LL k = (x <= Sqr) ? ind1[x] : ind2[n / x];
LL ans = (g2[k] - g1[k] + pre1[y] - pre2[y] + mod + mod) % mod;
for(Rint i = y + 1;i <= tot && pri[i] * pri[i] <= x;i ++){
LL pe = pri[i];
for(Rint e = 1;pe <= x;e ++, pe *= pri[i]){
LL xx = pe % mod;
ans = (ans + xx * (xx - 1) % mod * (S(x / pe, i) + (e > 1))) % mod;
}
}
return ans % mod;
}
int main(){
scanf("%lld", &n);
Sqr = sqrt(n);
init(Sqr);
for(Rint i = 1, last;i <= n;i = last + 1){
last = n / (n / i);
w[++ cnt] = n / i;
LL xx = w[cnt] % mod;
g1[cnt] = (xx * (xx + 1) / 2 + mod - 1) % mod;
g2[cnt] = (xx * (xx + 1) / 2 % mod * (2 * xx + 1) % mod * inv3 % mod + mod - 1) % mod;
if(n / i <= Sqr) ind1[w[cnt]] = cnt;
else ind2[last] = cnt;
}
for(Rint i = 1;i <= tot;i ++)
for(Rint j = 1;j <= cnt && pri[i] * pri[i] <= w[j];j ++){
LL k = (w[j] / pri[i] <= Sqr) ? ind1[w[j] / pri[i]] : ind2[n / (w[j] / pri[i])];
g1[j] -= pri[i] * (g1[k] - pre1[i - 1] + mod) % mod;
g2[j] -= pri[i] * pri[i] % mod * (g2[k] - pre2[i - 1] + mod) % mod;
if(g1[j] < 0) g1[j] += mod;
if(g2[j] < 0) g2[j] += mod;
}
printf("%lld", (S(n, 0) + 1) % mod);
}
讲那么多，写道题练手吧
UOJ188#【UR #13】Sanrd
题目链接：UOJ
这道题，也算是 min_25 的一个基础应用吧。。。
我们要求
，其中 (f(i))表示 (i)的次大质因子。
按照套路，我们设
所以
素数个数用 min_25 筛可以很快求出来。
#include<bits/stdc++.h>
#define Rint register LL
using namespace std;
typedef long long LL;
const int N = 1000003;
LL Sqr, pri[N], tot, g[N], w[N], id1[N], id2[N], cnt;
bool notp[N];
inline void init(int m){
notp[0] = notp[1] = true;
for(Rint i = 2;i <= m;i ++){
if(!notp[i]) pri[++ tot] = i;
for(Rint j = 1;j <= tot && i * pri[j] <= m;j ++){
notp[i * pri[j]] = true;
if(!(i % pri[j])) break;
}
}
}
inline LL solve(LL n, LL x, int y){
if(x <= 1) return 0;
LL ans = 0;
for(Rint k = y + 1;k <= tot && pri[k] * pri[k] <= x;k ++){
for(Rint pe = pri[k];pe * pri[k] <= x;pe *= pri[k]){
LL kk = (x / pe <= Sqr) ? id1[x / pe] : id2[n / (x / pe)];
ans += solve(n, x / pe, k) + pri[k] * (g[kk] - k + 1);
}
}
return ans;
}
inline LL solve(LL n){
Sqr = sqrt(n);
tot = cnt = 0;
init(Sqr);
for(Rint i = 1, last;i <= n;i = last + 1){
w[++ cnt] = n / i;
last = n / w[cnt];
g[cnt] = w[cnt] - 1;
if(w[cnt] <= Sqr) id1[w[cnt]] = cnt;
else id2[last] = cnt;
}
for(Rint i = 1;i <= tot;i ++)
for(Rint j = 1;j <= cnt && pri[i] * pri[i] <= w[j];j ++){
LL k = (w[j] / pri[i] <= Sqr) ? id1[w[j] / pri[i]] : id2[n / (w[j] / pri[i])];
g[j] -= g[k] - i + 1;
}
return solve(n, n, 0);
}
int main(){
LL l, r;
scanf("%lld%lld", &l, &r);
printf("%lld
", solve(r) - solve(l - 1));
}
参考
OI wiki：https://oi-wiki.org/math/min-25/
新版min25筛（O(n^(2/3))）详解：https://zhuanlan.zhihu.com/p/60378354
LaTeX数学公式大全：https://www.luogu.com.cn/blog/IowaBattleship/latex-gong-shi-tai-quan