Luogu 4213 【模板】杜教篩(Sum)
當作杜教篩的筆記吧。
杜教篩
要求一個積性函式$f(i)$的字首和,現在這個東西並不是很好算,那麼我們考慮讓它捲上另外一個積性函式$g(i)$,使$(f * g)$的字首和變得方便計算,然後再反推出這個$f$函式的字首和。
$$\sum_{i = 1}^{n}(f * g)(i) = \sum_{i = 1}^{n}\sum_{d | i}g(d)f(\frac{i}{d}) = \sum_{d = 1}^{n}g(d)\sum_{i = 1}^{\left \lfloor \frac{n}{d} \right \rfloor}f(i) = \sum_{d = 1}^{n}g(d)S(\frac{n}{d}{})$$
把$g(1)S(n)$移過來
$$g(1)S(n) = \sum_{i = 1}^{n}(f * g)(i) - \sum_{i = 2}^{n}g(i)S(\left \lfloor \frac{n}{i} \right \rfloor)$$
這個式子就是杜教篩的精髓了。
我們可以先篩出$[1, \sqrt{n}]$區間內的該積性函式的字首和,然後再分塊遞迴求解$(\sqrt{n}, n]$中的該函式的字首和,可以做到$O(n^{\frac{2}{3}})$的優秀的複雜度(並不會這個複雜度的證明)。
應該用一個雜湊表存一下已經計算過的各個$S(n)$的值($unordered\_map$)。
接下來的問題就是考慮如何搭配出這個積性函式$g$了。
模板題
考慮如何計算$\mu$和$\varphi$。
我們知道$\mu * I = \epsilon$,那麼有
$$S(n) = \sum_{i = 1}^{n}\epsilon(i) - \sum_{i = 2}^{n}S(\left \lfloor \frac{n}{i} \right \rfloor)$$
滑稽吧,$\epsilon$的字首和還不是$1$。
我們又知道$\varphi * I = id$,那麼又有
$$S(n) = \sum_{i = 1}^{n}id(i) - \sum_{i = 2}^{n}S(\left \lfloor \frac{n}{i} \right \rfloor)$$
而$\sum_{i = 1}^{n}id(i) = \sum_{i = 1}^{n}i = \frac{i(i + 1)}{2}$。
解決了!
Code:
#include <cstdio> #include <cstring> #include <unordered_map> using namespace std; typedef long long ll; const int N = 5e6 + 5; const int Maxn = 5e6; int testCase, pCnt = 0, pri[N], mu[N], phi[N]; ll sumMu[N], sumPhi[N]; bool np[N]; unordered_map <int, ll> sMu, sPhi; template <typename T> inline void read(T &X) { X = 0; char ch = 0; T op = 1; for (; ch > '9'|| ch < '0'; ch = getchar()) if (ch == '-') op = -1; for (; ch >= '0' && ch <= '9'; ch = getchar()) X = (X << 3) + (X << 1) + ch - 48; X *= op; } inline void sieve() { mu[1] = 1, phi[1] = 1; for (int i = 2; i <= Maxn; i++) { if (!np[i]) pri[++pCnt] = i, phi[i] = i - 1, mu[i] = -1; for (int j = 1; j <= pCnt && i * pri[j] <= Maxn; j++) { np[i * pri[j]] = 1; if (i % pri[j] == 0) { phi[i * pri[j]] = phi[i] * pri[j]; mu[i * pri[j]] = 0; break; } phi[i * pri[j]] = phi[i] * phi[pri[j]]; mu[i * pri[j]] = -mu[i]; } } for (int i = 1; i <= Maxn; i++) { sumMu[i] = sumMu[i - 1] + mu[i]; sumPhi[i] = sumPhi[i - 1] + phi[i]; } } ll getPhi(int n) { if (n <= Maxn) return sumPhi[n]; if (sPhi[n]) return sPhi[n]; ll res = 1LL * n * (n + 1) / 2; for (int l = 2, r; l <= n; l = r + 1) { r = (n / (n / l)); res -= 1LL * (r - l + 1) * getPhi(n / l); } return sPhi[n] = res; } ll getMu(int n) { if (n <= Maxn) return sumMu[n]; if (sMu[n]) return sMu[n]; ll res = 1LL; for (int l = 2, r; l <= n; l = r + 1) { r = (n / (n / l)); res -= 1LL * (r - l + 1) * getMu(n / l); } return sMu[n] = res; } int main() { sieve(); read(testCase); for (int n; testCase--; ) { read(n); printf("%lld %lld\n", getPhi(n), getMu(n)); } return 0; }View Code
感覺時限特別急,能別開$long \ long$就別開。