【学习笔记】Min_25筛

源自一个日本老哥神仙级别的构造,原链接在这里:$min\_25’s~blog$

然后 $min25$ 筛大概是对积性函数求和的一种方法。他求和的函数有以下特征:

1、$\boldsymbol{F}$ 是整值积性函数

2、广义上讲,需要 ${i\in \mathbb{P}},\boldsymbol{F}(i),\boldsymbol{F}(i^k)$ 这俩东西可以很简单地计算

3、然而一般情况下,都是 $f$ 这东西可以化作一个低次多项式。

嗯,第三条是大部分情况下都会把原函数拆分为 $\sum _{i\in\mathbb{P}} i^k$ 这种形式,因为这种形式的和比较好求。然而其实这么说来十分的naive,毕竟 $\varphi,\mu$ 也很好求。

大概思想就是神仙构造,比杜教筛神仙许多的构造。当然也有人管这个叫“扩展埃拉托色尼筛”…不过也确实很形象。

啊,学了好久啊…总是学不会呢,再聪明一点吧。

$1$ 构造一个 $g(n,j)$

本节先假设 $f$ 为 完全积性函数

嗯,这东西是这么一种形式 :

其中 $[]$ 还是艾佛森括号,$prime_j$ 表示线性筛出来的质数集合中第 $j$ 大的元素,$\mathrm{minf}(i)$ 表示 $i$ 的最小素因子。

直观地来说,$g(n,j)$ 表示的就是 $n$ 以内的在埃筛算法进行第 $i$ 轮后尚未被筛去的数的 $f$ 的和。

然后思考怎么转移,假设枚举到了 $j$:

1、如果当前 $prime_j^2>n$ ,那么对于当前的 $prime$ 不可能存在更多的 $\mathrm{minf}(i)> prime_j$ 。所以此时 $g(n,j)=g(n,j-1)$
2、如果当前 $prime_j^2\leq n$,即 $\lfloor \frac{n}{prime_j}\rfloor \geq prime_j$ 那么考虑应该减去那些最小素因子为 $prime_j$ 的数的 $f$ 值。 然而这个值发现本质上可以通过 $g$ 数组本身来刻画。即:

(1)首先,根据上文假设,$f$ 应该是一个完全积性函数,所以可以直接提出来。
(2)第二项的前一半是所有剔除掉 $prime_j$ 这一项之后,最小质因子 $>prime_{j-1}$ 即 $\geq prime_j$ 的数字的和,因为这些数字没有比 $prime_j$ 更小的质因子,所以要减去;后面是要减去前一项里面顺带计算了的 $<prime_j$ 的质数,因为显然 $1\sim prime_{j-1}$ 不可能存在以比 $prime_{j-1}$ 大的质数为最小质因子的数,所以累加的只有质数;而据定义这一部分需要倍计算进去,但是他们 $ \times ~prime_j$ 之后的数显然最小质因子不是 $prime_j$ 。

所以最终的递推式,为了简洁用 $p_j$ 替代 $prime_j$:

…鬼知道我为了理解这部分打了多久的摆花了多少时间和精力233

不过观察整s个 $g$ 的递推,可以发现只跟 $f$ 在质数处的取值有关。所以如果对于一个积性函数 $\boldsymbol{F}$, 我们可以将其拆分成几个完全积性函数 $\boldsymbol{F’,F’’….}$ 在 $x\in \mathbb{P}$ 处的乘积或者是加和,就可以实现快速计算 $g$。

这也正是上文中对于 $\boldsymbol{F}$ 的限制。如果 $\boldsymbol{F}$ 是个简单的低次多项式,那么就可以按照幂次拆分成几项分别求和再相加;或者更广义一点, $\boldsymbol{F}$ 可以拆分成许多在质数处以求得的函数,那么同样适用于 $min25$ 筛

$2$ 求和

终于要求和啦!

首先令 $\mathrm{S}(n,j)=\sum_{i}^{n} f(i) [\mathrm{minf}(i)>prime_j]$ ,那么最终结果就是 $\mathrm{S}(n,0)+f(1)$ 。最后的 $1$ 是因为 $1$ 没有质因子。

(ps : 然而似乎 $f(1)\equiv 1$?因为毕竟 $f(x)$ 是积性函数,那么就一定有 $f(x)\cdot f(1)=f(x)$……除了 $f(x)\equiv0$ 这种情况)

那么就可以知道 $\rm S$ 可以这么推:

前半部分是所有质数的贡献,同样由于不存在比最大的质数 $p_{\max}$ 更大的质因子,所以第一项就是统计了所有质数;后一项则是根据定义减掉了比 $j$ 小的质数的贡献。

后半部分计算合数。这个地方考虑枚举合数的最小质因子(此处显然 $>prime_j$,即枚举 $k>j$),然后利用积性函数的性质提取。值得注意的是最后有一个艾佛森括号,$[e>1]$,其意义是前文提到过的“ $1$ 没有质因子”,而当 $e=1$ 时,$prime_k^1\in \{prime\}$ , 是质数, 在前半部分计算过了。

嗯,于是这东西就可以暴力 $dfs$ 了。稳得很。

$3$ 复杂度

大家都说是 $\frac{n^{\frac{3}{4}}}{\ln n}$,然而神仙说其实应该是 $O(n^{1-\epsilon})$,是一种亚线性筛法,但感觉用在 $n\leq 1e10$ 还是可以 $1.5s$ 以内出结果的。

然后有个小细节,就是发现算 $g$ 这个东西的时候,空间开不下。然而其实发现

这东西说明了,$n$ 除以某个数得到的值其实是有限的。于是就可以直接当 $k>\sqrt n$ 时,用第一张表存 $\lfloor \frac{n}{k}\rfloor$ ;当 $k\leq \sqrt n$ 时,用第二张表存 $k$。同时 $g$ 的第二维也可以直接滚掉。发现这样空间复杂度就做到了 $O(\sqrt n)$ 。

$4$ 例题

LG5325 【模板】Min_25筛

给定积性函数 $f$,满足 $p\in \mathbb{P},f(p^k)=p^k(p^k-1)$ 。求和

$n\leq 10^{10}$。

发现这就是比较套路的多项式拆成单项式,所以就直接拆成平方项和一次项做就好了。

然后递推 $g$ 之前要先把 $g(n,0)$ 算出来,这东西完全可以 $\sqrt n$ 根据定义直接算。

。。。还有,注意取模

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
LL g1[N], g2[N] ;
LL n, s1[N], s2[N] ;
LL dex[2][N] ; int tot ;
const int P = 1000000007 ;
const LL Inv6 = 166666668ll ;
const LL Inv2 = 500000004ll ;
int s, pr[N], cnt, chk[N] ; LL Id[N] ;

void pre_sieve(int w){
chk[0] = chk[1] = 1 ;
for (int i = 2 ; i <= w ; ++ i){
if (!chk[i])
pr[++ cnt] = i,
s1[cnt] = (s1[cnt - 1] + i) % P,
s2[cnt] = (s2[cnt - 1] + 1ll * i * i % P) % P ;
for (int j = 1 ; j <= cnt ; ++ j){
if (i * pr[j] > w) break ;
chk[i * pr[j]] = 1 ;
if (i % pr[j] == 0) break ;
}
}
}
LL calc1(LL x){
x %= P ;
return x * (x + 1) % P * Inv2 % P ;
}
LL calc2(LL x){
x %= P ;
return x * (x + 1) % P * (2 * x + 1) % P * Inv6 % P ;
}
LL S(LL p, int q){
if (pr[q] >= p) return 0 ;
LL id = (p <= s ? dex[0][p] : dex[1][n / p]) ;
LL ans = ((g2[id] - g1[id] - s2[q] + s1[q]) % P + P) % P ;
//上一句本质上就是
for (int k = q + 1 ; k <= cnt ; ++ k){
if (1ll * pr[k] * pr[k] > p) break ; LL o, t, pq ;
for (o = 1, pq = pr[k] ; pq <= p ; ++ o, pq = 1ll * pr[k] * pq)
t = pq % P,
ans = (ans + t * (t - 1) % P * (S(p / pq, k) % P + (o > 1))) % P ;
}
return ans ;
}
int main(){
cin >> n,
s = sqrt(n) + 1,
pre_sieve(s) ; LL l, r, w ;
for (l = 1 ; l <= n ; l = r + 1){
r = n / (n / l),
Id[++ tot] = n / r ;
g1[tot] = (calc1(n / l) - 1 + P) % P ;
g2[tot] = (calc2(n / l) - 1 + P) % P ;
if (n / l > s)
dex[1][l] = tot ;
else dex[0][n / l] = tot ;
}
for (int i = 1 ; i <= cnt ; ++ i){
w = 1ll * pr[i] * pr[i] ; LL now ;
for (int j = 1 ; j <= tot && w <= Id[j] ; ++ j){
now = Id[j] / pr[i] ;
now = (now <= s) ? dex[0][now] : dex[1][n / now] ;
(g1[j] -= 1ll * pr[i] * (g1[now] - s1[i - 1] + P) % P) %= P ;
(g2[j] -= 1ll * pr[i] * pr[i] % P * (g2[now] - s2[i - 1] + P) % P) %= P ;
g1[j] += (g1[j] < 0) ? P : 0 ;
g2[j] += (g2[j] < 0) ? P : 0 ;
}
}
/*for (int i = 1 ; i <= tot ; ++ i)
printf("%lld %lld\n", g1[i], g2[i]) ; */
// cout << 233 << endl ;
cout << (S(n, 0) + 1) % P << endl ; return 0 ;
}

LG4213 【模板】杜教筛

求 $\sum \varphi(i)$ 和 $\sum \mu(i)$ .

$n< 2^{31}$

发现 $p\in \mathbb{P}$ ,$\varphi(p)=p-1,\varphi(p^k)=p^{k-1}(p-1)$,可以快速求,于是可以用一个 $g$ 来推 $f(x)=x$ 的前缀质数和,一个 $g$ 来推 $1$ ;同样对于 $p\in \mathbb{P},\mu(p)=-1,\mu(p^{k})=0 $ 所以这东西也可以通过递推 $1$ 来实现。

代码方面,虽然比自己实现的杜教筛要快,但是比其他人的杜教筛要慢不少。。。awsl

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
#include <cmath>
#include <cstdio>
#include <iostream>

using namespace std ;

#define N 200010
#define rg register
#define LL long long

using namespace std ;

int n ; LL sp[N] ;
LL gmu[N], gphi[N] ;
LL dex[2][N] ; int T, tot ;
int s, pr[N], cnt, chk[N] ; LL Id[N] ;

void pre_sieve(int w){
chk[0] = chk[1] = 1 ;
for (rg int i = 2 ; i <= w ; ++ i){
if (!chk[i])
pr[++ cnt] = i,
sp[cnt] = sp[cnt - 1] + i ;
for (rg int j = 1 ; j <= cnt ; ++ j){
if (i * pr[j] > w) break ;
chk[i * pr[j]] = 1 ;
if (i % pr[j] == 0) break ;
}
}
}
inline LL calc1(LL x){
return x * (x + 1) / 2 ;
}
LL Sphi(LL p, int q){
if (pr[q] >= p) return 0 ;
int id = (p <= s ? dex[0][p] : dex[1][n / p]) ;
rg LL ans = (gphi[id] - gmu[id]) - (sp[q] - q) ;
//上一句本质上就是g(n,|\{prime \}|) - g(prime_{j-1},j-1)
for (int k = q + 1 ; k <= cnt ; ++ k){//之后就是枚举更大的素因子(枚举合数)
if (1ll * pr[k] * pr[k] > p) break ; LL o, t, pq ;
for (o = 1, pq = pr[k] ; pq <= p ; ++ o, pq = 1ll * pr[k] * pq)
ans = (ans + (pq / pr[k]) * (pr[k] - 1) * (Sphi(p / pq, k) + (o > 1))) ;
}
return ans ;
}
LL Smu(LL p, int q){
if (pr[q] >= p) return 0 ;
int id = (p <= s ? dex[0][p] : dex[1][n / p]) ;
LL ans = -gmu[id] + q ;
for (int k = q + 1 ; k <= cnt ; ++ k){
if (1ll * pr[k] * pr[k] > p) break ;
ans += -(Smu(p / pr[k], k)) ;
}
return ans ;
}
int main(){
cin >> T ;
pre_sieve(N - 1) ;
while (T --){
scanf("%d", &n) ; tot = 0 ;
rg LL l, r, w, t ; s = sqrt(n) ;
for (l = 1 ; l <= n ; l = r + 1){
//先用整除分块预处理g(n,0)
t = n / l ;
r = n / t ;
Id[++ tot] = t ;//(保存每一项,最多只有2*sqrt(n)个)
gmu[tot] = t - 1 ;//预处理1,即在j=0时,就是该函数在所有的素数的取值-1
gphi[tot] = calc1(t) - 1 ;//预处理n
if (n / l > s)
dex[1][l] = tot ;
else dex[0][n / l] = tot ;//空间节省
}
for (int i = 1 ; i <= cnt ; ++ i){
w = 1ll * pr[i] * pr[i] ; LL now ;
for (int j = 1 ; j <= tot && Id[j] >= w ; ++ j){
now = Id[j] / pr[i] ;
now = (now <= s) ? dex[0][now] : dex[1][n / now] ;
gmu[j] -= gmu[now] - (i - 1), //按照递推式求g1和g2
gphi[j] -= 1ll * pr[i] * (gphi[now] - sp[i - 1]) ;
}
}
printf("%lld %lld\n", Sphi(n, 0) + 1, Smu(n, 0) + 1) ;
}
}

$5$ 一点后话

在写题的过程中发现了另一种写法。大概就是原本 $\rm S$ 的定义是这样的:

然后他们推导的是这样的:

这样最后的答案就是 $\mathrm{S}(n,1)+1$ ,在写的时候也需要后移一位。