杜教筛是一种在非线性时间内求积性函数前缀和的算法。
前置知识:
积性函数:
数论函数:定义域为正整数的函数。
积性函数:对于任意互质的整数$a,b$,有$f(ab) = f(a)f(b)$的数论函数
完全积性函数:对于任意整数$a,b$,有$f(ab) = f(a)f(b)$的数论函数
常见的积性函数:$\phi, \mu, \sigma, d$(分别表示:欧拉函数,莫比乌斯函数,因子和函数,因子个数函数)
常见的完全积性函数:$\epsilon, I, id$($\epsilon(n) = [n = 1], I(n) = 1, id(n) = n$)
狄利克雷卷积:
设$f,g$是两个数论函数,它们的狄利克雷卷积为:$(f*g)(n) = \sum_{d|n}f(d)g(\frac{n}{d})$
运算律:满足交换律,结合律
单位元:$\epsilon$($f * \epsilon = f$)
几个重要的性质:
- $\mu * I = \epsilon$
- $\phi * I = id$
- $\mu * id = \phi$
利用狄利克雷卷积,我们能很容易的证明莫比乌斯反演:
若:
$$g(n) = \sum_{d|n}f(d) => g = f * I$$
则:
$$f(n) = \sum_{d|n}\mu(d)g(\frac{n}{d}) => \mu * g = \mu * I * f = \epsilon * f = f$$
杜教筛:
下面就是正题了,我们现在要求积性函数$f(i)$的前缀和,设$\sum_{i=1}^{n} f(i) = S(n)$
再找一个积性函数$g$,考虑$g * f$的前缀和:
$$\sum_{i=1}^{n} (g * f)(i)$$
$$= \sum_{i=1}^{n} \sum_{d|i} g(d) * f(\frac{i}{d})$$
$$= \sum_{d=1}^{n} g(d) \sum_{i=1}^{\frac{n}{d}} f(\frac{di}{d})$$
$$= \sum_{d=1}^{n} g(d) S(\frac{n}{d})$$
上面这个式子启发我们去考虑$g(1)S(n)$,容易得到:
$$g(1)S(n) = \sum_{i=1}^{n} g(i) S(\frac{n}{i}) - \sum_{j=2}^{n} g(j) S(\frac{n}{j})$$
即:
$$g(1)S(n) = \sum_{i=1}^{n} (g * f)(i) - \sum_{j=2}^{n} g(j) S(\frac{n}{j})$$
至此,我们已经得到了杜教筛的核心式子。
假如我们选取一个合适的积性函数$g$,使得我们能够快速求解$g$的前缀和以及$g * f$的前缀和,那么我们就可以利用数论分块+递归求解上面的式子。
1 ll GetSum(int n) { // 算 f 前缀和的函数 2 ll ans = f_g_sum(n); // 算 f * g 的前缀和 3 // 以下这个 for 循环是数论分块 4 for(ll l = 2, r; l <= n; l = r + 1) { // 注意从 2 开始 5 r = (n / (n / l)); 6 ans -= (g_sum(r) - g_sum(l - 1)) * GetSum(n / l); 7 // g_sum 是 g 的前缀和 8 // 递归 GetSum 求解 9 } return ans; 10 }
可以证明,这样做的复杂度是$O(n^{\frac{3}{4}})$的。
我们还可以进一步优化,即先线性筛出前$m$个答案,然后再利用杜教筛求解。
当$m = n^{\frac{2}{3}}$时,复杂度是$O(n^{\frac{2}{3}})$
模板题:杜教筛
我们先考虑$\mu$的前缀和,注意到$\mu * I = \epsilon$,于是我们取$g = I, f = \mu$。
那$\epsilon, I$的前缀和都是可以$O(1)$ 求的,直接套用杜教筛的式子即可。
然后是求$\phi$的前缀和,因为有$\phi * I = id$,所以取$g = I, f = \phi$。
这里要注意一下,当$n = 2^{31} - 1$的时候可能会在计算过程中爆$ll$,转换成$ull$运算即可。
1 #include<bits/stdc++.h> 2 #include <tr1/unordered_map> 3 using namespace std; 4 using namespace std::tr1; 5 6 typedef long long ll; 7 typedef unsigned long long ull; 8 9 const int inf = 0x7FFFFFFF; 10 const int MAXN = 7000000; 11 12 int cnt, prim[MAXN + 10]; 13 bool v[MAXN + 10]; 14 ll mu[MAXN], phi[MAXN]; 15 unordered_map<int, ll> ansmu, ansphi; 16 17 void init() 18 { 19 v[1] = mu[1] = phi[1] = 1; 20 for(int i = 2; i <= MAXN; ++i) 21 { 22 if(!v[i]) prim[++cnt] = i, mu[i] = -1, phi[i] = i - 1; 23 for(int j = 1; j <= cnt && i * prim[j] <= MAXN; ++j) 24 { 25 v[i * prim[j]] = 1; 26 if(i % prim[j]) mu[i * prim[j]] = -mu[i], phi[i * prim[j]] = phi[i] * phi[prim[j]]; 27 else {mu[i * prim[j]] = 0, phi[i * prim[j]] = phi[i] * prim[j]; break;} 28 } 29 } 30 for(int i = 2; i <= MAXN; ++i) mu[i] += mu[i - 1], phi[i] += phi[i - 1]; 31 } 32 33 ll sum_mu(int n) 34 { 35 if(n <= MAXN) return mu[n]; 36 if(ansmu[n]) return ansmu[n]; 37 ll ret = 0; 38 for(int l = 2, r = 0; r < inf && l <= n; l = r + 1) 39 { 40 r = n / (n / l); 41 ret += (r - l + 1) * sum_mu(n / l); 42 } 43 return ansmu[n] = 1ll - ret; 44 } 45 46 ll sum_phi(int n) 47 { 48 if(n <= MAXN) return phi[n]; 49 if(ansphi[n]) return ansphi[n]; 50 ll ret = 0; 51 for(int l = 2, r = 0; r < inf && l <= n; l = r + 1) 52 { 53 r = n / (n / l); 54 ret += (r - l + 1) * sum_phi(n / l); 55 } 56 return ansphi[n] = (ull)n * (n + 1ull) / 2ull - ret; 57 } 58 59 int main() 60 { 61 init(); int T, n; 62 scanf("%d", &T); 63 while(T--) 64 { 65 scanf("%d", &n); 66 printf("%lld %lld\n", sum_phi(n), sum_mu(n)); 67 } 68 return 0; 69 }