【学习笔记】杜教筛

大概是一种比较快的筛法?

用于积性函数求和。比如给定了一个函数 $\boldsymbol{\rm f} $ ,考虑构造出这样的一个 $\boldsymbol{\rm g}$ ,使得在如此转化之后:

得出来的:

如果其中 $\bf f*g$ 比较容易计算,就可以达到快速计算 $\mathbf{S}$ 的目的了。

似乎很简单?因为这种题的难点在于构造,构造完了就可以分块做了

1
2
3
4
5
6
ll solve(ll n){
ll s = sum(1, n, f * g) ;
for (int l = 1, r ; l <= n ; l = r + 1)
r = n / (n / l), s -= solve(n / l) * sum(1, l, g) ;
return (s /= g(1)) ;
}

然而这么做的复杂度不对。考虑这样一个引理:

证明大概是把 $z$ 设成带余除法的标准形式,然后搞一下就可以了。此处略去。

于是发现,其实每次要计算的 $n’$ ,都是某个 $\lfloor \frac{n}{x}\rfloor$ ,因为假设这次是 $\lfloor \frac{n}{x}\rfloor$ ,下次是 $\lfloor \frac{\lfloor \frac{n}{x}\rfloor}{y}\rfloor=\lfloor \frac{n}{x\cdot y}\rfloor$ 。所以 $\sqrt n$ 为上界,可以直接通过 $\lfloor \frac{n}{x}\rfloor $ 记忆化所有的 $x$ 。

同时,发现对于 $x$ 的调用,都只会是

这些数,且每次计算 $\mathbf{S}$ 的复杂度都是 $\sqrt n$ 。所以最终复杂度是:

其实就是两步近似。

然而可以做到更优,即发现对于一些很小的 $n’$ ,完全可以直接预处理出来。所以假设对于前 $p(0<p<1)$ 的 $\mathbf S$ 直接预处理,那么复杂度会变成:

发现此时取 $p=\frac{2}{3}$ 为最优,所以复杂度 $O(n^{\frac{2}{3}})$。

草,据说这东西是 xudyh 在考场上构造出来的,这也太神了吧…

$1$ 例题

LG4213 【模板】杜教筛(Sum)

$1\leq n\leq 2\times 10^9$

考虑杜教筛。

考虑如何构造 $\varphi$ 的 $g$ 。发现有 $\sum_{d|n}\varphi(d)=n$,即 $\varphi \mathbf{1}=\bf Id$,那么 $\bf fg$ 就是 $\frac{n\cdot (n+1)}{2}$。

考虑如何构造 $\mu$ 的 $g$ 。发现有 $\sum_{d|n}\mu(d)=[n=1]$,即 $\mu \bf 1 = ϵ$,那么 $\bf fg$ 就是 $[n=1]$。

然后就没有然后了。去年 $6$ 月我写的是不带记忆化的,慢的很。今天重写加了记忆化,复杂度看起来还可以?但是不知道为啥被 min_25 吊着锤。(其实也就快一点

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
#define rr register  
#define MAXS 100010
#define MAXN 5200010
#define LL long long

using namespace std ;
int T, N, S ;
bool chk[MAXN] ;
LL phi[MAXN], Mu[MAXN] ;
LL Smu[MAXS], Sphi[MAXS] ;
int pr[MAXN], t, i, j, x ;

void init(){
chk[0] = chk[1] = 1 ;
Mu[1] = 1, phi[1] = 1 ;
for (i = 2 ; i < MAXN ; ++ i){
if (!chk[i])
pr[++ t] = i, Mu[i] = -1, phi[i] = i - 1 ;
for (j = 1 ; j <= t && pr[j] * i < MAXN ; ++ j){
chk[i * pr[j]] = 1 ;
if (i % pr[j] == 0) {
phi[i * pr[j]] = phi[i] * pr[j] ;
break ;
}
else
Mu[i * pr[j]] = - Mu[i],
phi[i * pr[j]] = phi[i] * phi[pr[j]] ;
}
}
for (i = 1 ; i < MAXN ; ++ i) Mu[i] = Mu[i - 1] + Mu[i], phi[i] = phi[i - 1] + phi[i] ;
}
inline LL get_Mu(int p){
if (p < MAXN) return Mu[p] ;
if (p <= S && Smu[p]) return Smu[p] ;
if (p > S && Smu[N / p]) return Smu[N / p] ;

LL ret = 1LL ; rr int l, r ;
for (l = 2 ; l <= p ; l = r + 1)
r = (p / (p / l)), ret -= (r - l + 1) * get_Mu(p / l) ;
return (p > S ? (Smu[N / p] = ret) : (Smu[p] = ret));
}
inline LL get_phi(int p){
if (p < MAXN) return phi[p] ;
if (p <= S && Sphi[p]) return Sphi[p] ;
if (p > S && Sphi[N / p]) return Sphi[N / p] ;

LL ret = 1ll * (1 + p) * p / 2LL ; rr int l, r ;
for (l = 2 ; l <= p ; l = r + 1)
r = (p / (p / l)), ret -= (r - l + 1) * get_phi(p / l) ;
return (p > S ? (Sphi[N / p] = ret) : (Sphi[p] = ret));
}
int main(){
cin >> T, init() ;
while (T --){
cin >> N, S = sqrt(N) ;
memset(Smu, 0, sizeof(Smu)) ;
memset(Sphi, 0, sizeof(Sphi)) ;
cout << get_phi(N) << " " << get_Mu(N) << endl ;
}
}

$3$ 高级应用

高级应用:我不是很会的应用(

$p\in \mathbb{P},1\leq n\leq 10^{10}$ .

首先是朴素反演:

然后发现前面都是完全积性,后面那个平方里面还带着求和不是很好搞…

这个地方有个很牛的引理:

这东西大概是能归纳出来的。似乎还没有什么比较有趣的证明约等于只能选择背过

然后就是

发现都比较容易算,并且都积性,于是就可以直接杜教筛了。

上句话在扯 $p$。后面那一坨根本没法筛。于是考虑整除分块套一个杜教筛。这样复杂度还是 $n^{\frac{2}{3}}$ !!这点很重要!!因为本质上我的整除分块还是在枚举某些 $\lfloor\frac{\lfloor\frac{n}{l}\rfloor}{x}\rfloor$。于是该记下来的不多不少,还是记下来了,复杂度不变。

代码实现有点细节,还需要多学习多积累。

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
#include <map>
#include <cstdio>
#include <iostream>

#define MAXN 100100
#define MAXNN 10000100
#define LL long long

using namespace std ;

map<LL, LL> pq ;
bool vis[MAXNN], chk[MAXN] ;
LL Inv6, Inv4, cnt, p[MAXNN] ;
LL P ; LL N, Sp[MAXN], phi[MAXNN] ;

LL expow(LL a, int b){
LL res = 1 ;
while (b){
if (b & 1)
(res *= a) %= P ;
(a *= a) %= P, b >>= 1 ;
}
return res ;
}
LL calc(LL l, LL r){
l -- ; l %= P, r %= P ;
LL r1 = r * (r + 1) % P * (2 * r + 1) % P * Inv6 % P ;
LL r2 = l * (l + 1) % P * (2 * l + 1) % P * Inv6 % P ;
return ((r1 - r2) % P + P) % P ;
}
LL calcc(LL x){
LL ans ; x %= P ;
ans = x * x % P ;
(ans *= Inv4) %= P ;
(ans *= (x + 1) % P * (x + 1) % P) %= P ;
return ans ;
}
LL solve(LL x){
if (x < MAXNN)
return phi[x] ;
LL y = N / x ;
if (chk[y]) return Sp[y] ;
chk[y] = 1 ; LL &ans = Sp[y] ;
ans = calcc(x) ;
for (LL l = 2, r ; l <= x ; l = r + 1){
r = x / (x / l),
(ans -= calc(l, r) * solve(x / l) % P) %= P ;
if (ans < 0) ans += P ;
}
// cout << ans << endl ;
return Sp[y] ;
}
LL work(LL x){
LL res = 0 ;
for (LL l = 1, r ; l <= x ; l = r + 1){
r = x / (x / l) ;
(res += (solve(r) - solve(l - 1)) * calcc(x / l) % P) %= P ;
}
return (res % P + P) % P ;
}
int main(){
cin >> P >> N ; phi[1] = vis[0] = vis[1] = 1 ;
Inv4 = expow(4, P - 2), Inv6 = expow(6, P - 2) ;
for (int i = 2 ; i < MAXNN ; ++ i){
if (!vis[i])
p[++ cnt] = i, phi[i] = i - 1 ;
for (int j = 1 ; j <= cnt ; ++ j){
if (p[j] * i >= MAXNN) break ;
vis[p[j] * i] = 1 ;
if (i % p[j] == 0) {
phi[i * p[j]] = phi[i] * p[j] % P ;
break ;
}
else phi[i * p[j]] = phi[i] * phi[p[j]] % P ;
}
}
for (int i = 1 ; i < MAXNN ; ++ i)
(phi[i] = 1ll * i * i % P * phi[i] + phi[i - 1]) %= P ;
//cout << phi[8000000] << endl ;
printf("%lld\n", work(N)) ; return 0 ;
}