杜教篩
例題
我們先來看一道題目。
求 \( S(n) = \displaystyle\sum_{i=1}^n \mu(i) \) 的值。(\( 1\le n\le 10^{10} \))
由於 \( n \) 太大了,即使用線性篩也會 TLE。但是先回想一個引理:
\( \mu\ast 1 = \sum_{d\mid n}\mu(d)=[n=1] \)
對它取前綴和,會發現它很好算,也就是 \(\sum_{i=1}^n \mu\ast 1 = \sum_{i=1}^n [i=1] = 1\)。
因此若可以找到 \(\mu\ast 1\) 的前綴和,還有 \(\mu\) 的前綴和之間的關係,說不定就能夠加速計算 \(\mu\) 的前綴和。
而這就是大名鼎鼎的杜教篩。
杜教篩
杜教篩適用於以下情況:
給定函數 \( f \),要求 \( \displaystyle S(n)=\sum_{i=1}^n f(n) \)。
定理:
\( \displaystyle g(1)\times S(n)=\sum_{i=1}^n (g\ast f)(i) - \sum_{i=2}^n g(i)\times S\left(\left\lfloor\frac n i\right\rfloor\right) \)
證明如下:
\begin{align} \displaystyle \sum_{i=1}^n (g\ast f)(i) &= \sum_{i=1}^n \sum_{d\mid i} g(d)\times f\left(\frac{i}{d}\right) & \text{根據 Dirichlet 捲積的定義展開} \\ &= \sum_{d=1}^n \sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor} g(d)\times f(i) & \text{套路,將 sigma 的順序交換,並改變枚舉的變數} \\ &= \sum_{d=1}^n g(d) \sum_{i=1}^{\left\lfloor\frac{n}{d}\right\rfloor} f(i) & g(d) \text{ 與 } i \text{ 無關,將其移到前面} \\ &= \sum_{d=1}^n g(d)\times S(\left\lfloor\frac{n}{d}\right\rfloor) & \text{將 } \sum f \text{ 換成前綴和} \\ &= g(1)\times S(n) + \sum_{d=2}^n g(d)\times S(\left\lfloor\frac{n}{d}\right\rfloor) & \text{將 } d = 1 \text{ 的情況獨立出來} \end{align}
移項後,即可得到:
\begin{align} \displaystyle g(1)\times S(n)=\sum_{i=1}^n (g\ast f)(i) - \sum_{i=2}^n g(i)\times S\left(\left\lfloor\frac n i\right\rfloor\right) \end{align}
\( \blacksquare \)
至此可以發現,若我們有個快速計算 RHS 的方法,那只要再除掉 \(g(1)\),就能夠快速得到 \(S(n)\) 了。
當然的,\( g(1) \) 必須是非零,這點後面不再贅述。
寫成 pseudo code 可以得到:
def Find_S(n):
ans = sum of convolution(g, f)(1...n)
for i in range(2, n):
ans = ans - g(i) * Find_S(n / i)
return ans / g(1)
注意到,for
迴圈的部分,只要我們能夠快速求出 \( g \) 的前綴,那麼就可以數論分塊。
令 \( Sg \) 是 \( g \) 的前綴和。
def Find_S(n):
ans = sum of convolution(g, f)(1...n)
i = 1
while i <= n:
r = n / (n / i)
ans = ans - (Sg(r) - Sg(i - 1)) * Find_S(n / i)
i = r + 1
return ans / g(1)
因此,杜教篩的流程其實就是不斷透過上述定理還有數論分塊,來進行遞迴計算。
Find_S
是一個遞迴函式,接下來我們會討論在遞迴中,實際會走訪的狀態數是多少。
首先我們需要一個引理:
對任意正整數 \( n、a、b \),都有 \( \displaystyle\left\lfloor\frac{\left\lfloor\frac{n}{a}\right\rfloor}{b}\right\rfloor =\left\lfloor\frac{n}{ab}\right\rfloor \)。
證明:
令 \( \frac n a = x + y \),其中 \( x \) 是非負整數、\( y \) 是小數部分。
則 \[ \text{LHS} = \left\lfloor\frac x b\right\rfloor , \text{RHS} = \left\lfloor\frac{x + y}{b}\right\rfloor \]
假設 \( \text{LHS}\neq\text{RHS} \),那麼 \( \text{LHS} < \text{RHS} \)。
因為兩者皆是非負整數,所以在相異的假設下,至少相差一,由此可得:
\[ \left\lfloor\frac x b\right\rfloor \le \left\lfloor\frac{x + y}{b}\right\rfloor - 1 = \left\lfloor\frac{x + y - b}{b}\right\rfloor \Rightarrow x \le x+y-b \]
然而 \( y \) 是小數部分、\( b \) 是正整數,也就是 \( y < 1 \Rightarrow y-b < 0\),得到 \( x \) 比自己還要小,產生矛盾。
因此由反證法可知引理是正確的。
\( \blacksquare \)
接著我們需要定義一種集合:
\( \displaystyle R(n)=\left\{\left\lfloor\frac n i\right\rfloor :\ i=1,2,\cdots, n\right\} \)
翻譯成白話,\( R(n) \) 就是 \( n \) 除以所有比其小的正整數並取下高斯後,那些數字們所組成的集合。
而我們在數論分塊那章節時就知道,這個集合的大小是 \(\mathcal{O}(\sqrt{n})\)。
有了這個定義,我們就可以來看下面這個定理:
若 \( m\in R(n) \),則 \( R(m)\subseteq R(n) \)。
證明:
根據定義,\( m\in R(n) \) 代表存在某個正整數 \( i \) (\( 1\le i\le n \)) 使得 \( \left\lfloor\dfrac n i\right\rfloor = m \)。
因此
\[ R(m)=\left\{\left\lfloor\frac{\left\lfloor\frac{n}{i}\right\rfloor}{j}\right\rfloor :\ j=1, 2, \cdots, m \right\} \]
利用前述引理,可以簡化成:
\[ R(m)=\left\{\left\lfloor\frac{n}{ij}\right\rfloor :\ j=1, 2, \cdots, m \right\} \]
由此可見 \( R(m)\subseteq R(n) \)。
\( \blacksquare \)
有了這份定理,我們回顧 Find_S
。
Find_S
是遞迴函數,要計算 Find_S(n)
會需要遞迴到所有 \( i = 2\cdots n \) 的 Find_S(n / i)
。
注意到,也就是所有 \( i\in R(n)\backslash\{1\} \) 的 Find_S(n / i)
。
令 \( m = \frac n i \)。
計算 Find_S(m)
時也會需要用到所有 \( j\in R(m)\backslash\{1\} \) 的 Find_S(m / j)
。
根據定理可知,所需要的那些 \( j \),其實也會屬於 \( R(n)\backslash\{1\} \)。
也就是說,在計算 Find_S(n)
的遞迴中,所遇到的任何 Find_S(x)
,\( x \) 都會屬於 \( R(n) \)。
因為 \( |R(n)|=\mathcal{O}(\sqrt{n}) \),所以遞迴時所需的不同狀態數只有 \( \mathcal{O}(\sqrt{n}) \) 個,可以記憶化來計算。
def Find_S(n):
if Find_S(n) has been visited:
return ans[n]
ans[n] = sum of convolution(g, f)(1...n)
i = 1
while i <= n:
r = n / (n / i)
ans[n] = ans[n] - (Sg(r) - Sg(i - 1)) * Find_S(n / i)
i = r + 1
ans[n] = ans[n] / g(1)
return ans[n]
時間複雜度
上個章節討論了杜教篩的原理,並證明了會遞迴到的狀態只有 \( \mathcal{O}(\sqrt{n}) \),這裡我們要來討論時間複雜度。
在此我們有兩個假設:
- \( g\ast f \) 的前綴和很好計算
- \( g \) 的前綴和很好計算
方便起見,我們把它們當成可以 \( \mathcal{O}(1) \) 得到。
那麼問題就只剩下計算 \( \sum_{i=2}^n g(i)\times S\left(\left\lfloor\frac n i\right\rfloor\right) \) 的時間複雜度了。
我們需要所有 \( i\in R(n) \) 的 Find_S(n / i)
。
而計算 Find_S(x)
會需要用到所有 \( j\in R(x) \) 的 Find_S(x / j)
,也就是需要 \( \mathcal{O}(\sqrt{x}) \) 個狀態。
合併下來可知時間複雜度是:
\begin{align} \sum_{x\in R(n)} \sqrt{x} &= \sum_{x=1}^{\sqrt{n}} \left(\sqrt{x} + \sqrt{\frac{n}{x}}\right) & \text{按照 } \sqrt{n} \text{ 分成兩邊} \\ &=\int_1^\sqrt{n} \left(\sqrt{x} + \sqrt{\frac{n}{x}}\right) dx & \text{使用微積分來估算} \\ &=\left(\left.\frac{2}{3}\times x^{\frac{3}{2}} \right|_1^\sqrt{n}\right) + \left(\left.2\times\sqrt{n}\times x^{\frac{1}{2}}\right|_1^\sqrt{n}\right) & \text{透過微積分運算得到} \\ &=\mathcal{O}(n^{\frac{3}{4}}) & \end{align}
因此,純粹使用杜教篩的話,時間複雜度為 \( \mathcal{O}(n^{\frac{3}{4}}) \)。
更進一步,我們可以選定某個數字 \( k\ge \sqrt{n} \),將 \( S[1]\sim S[k] \) 使用線性篩來得到。
因此在 Find_S
中,我們只需要計算 \( R(n) \) 中大於 \( k \) 的元素即可。
因此時間複雜度會降為:
\begin{align} \sum_{x\in R(n) \\ x > k} \sqrt{x} &= \sum_{x=1}^{\frac{n}{k}} \sqrt{\frac{n}{x}} \\ &= \int_1^\frac{n}{k} \sqrt{\frac{n}{x}} dx \\ &= \left.2\times \sqrt{n} \times x^{\frac{1}{2}}\right|_1^\frac{n}{k} \\ &= \mathcal{O}\left(\frac{n}{\sqrt{k}}\right) \end{align}
結合線性篩,總時間複雜度是:\( \mathcal{O}\left(k + \frac{n}{\sqrt{k}}\right) \)。
最佳情況發生在 \( k = \frac{n}{\sqrt{k}} \),也就是 \( k = n^{\frac{2}{3}} \) 的時候。
因此杜教篩的時間複雜度,就進一步降低為 \( \mathcal{O}(n^\frac{2}{3}) \) 了。
例題
現在我們重新回顧最開始所講的例題。
求 \( S(n) = \displaystyle\sum_{i=1}^n \mu(i) \) 的值。(\(1\le n\le 10^{10}\))
由杜教篩的定理可知,我們需要找到一個函數 \( g \),使得:
- \( g\ast\mu \) 的前綴和很好計算
- \( g \) 的前綴和很好計算
而我們知道 \( (\mu\ast 1)(n) = [n=1]\),其前綴和很顯然,就是 \( \sum_{i=1}^n [i=1] = 1 \)。
\( 1(n) \) 的前綴和也是顯然,就是 \( \sum_{i=1}^n 1(i) = n \)。
因此可以將 \( g \) 函數選為 \( 1 \) 函數,因此可得:
\[ S(n) = 1 - \sum_{i=2}^n S\left(\left\lfloor\frac{n}{i}\right\rfloor\right) \]
利用這個遞迴式,寫成程式碼如下:
#include <bits/stdc++.h>
using namespace std;
using ll = long long;
ll n;
ll lim; // n^(2/3)
unordered_map<ll, ll> s; // prefix sum of mu
vector<ll> psm; // prefix sum of mu
void LinearSieve() {
vector<int> primes;
vector<int> pf(lim); // prime factor
for (int i = 2; i < lim; ++i) {
if (pf[i] == 0) {
primes.push_back(i);
pf[i] = i;
}
for (auto p : primes) {
if (i * p >= lim)
break;
pf[i * p] = p;
if (i % p == 0)
break;
}
}
vector<int> mu(lim);
psm.resize(lim);
mu[1] = psm[1] = 1;
for (int i = 2; i < lim; ++i) {
if (i % (pf[i] * pf[i]) != 0)
mu[i] = -mu[i / pf[i]];
psm[i] = psm[i - 1] + mu[i];
}
}
ll Find_S(ll n) {
if (n < lim)
return psm[n];
if (s.find(n) != s.end())
return s[n];
s[n] = 1;
for (ll i = 2, r = 0; i <= n; i = r + 1) {
r = n / (n / i);
s[n] -= (r - i + 1) * Find_S(n / i);
}
return s[n];
}
int main() {
cin >> n;
while ((__int128)lim * lim * lim <= (__int128)n * n)
++lim;
LinearSieve();
cout << Find_S(n) << "\n";
}
或者也可以使用迭代來實作:
#include <bits/stdc++.h>
using namespace std;
using ll = long long;
int lim;
unordered_map<ll, ll> s; // prefix sum of mu
vector<ll> psm; // prefix sum of mu
void LinearSieve() {
vector<int> primes;
vector<int> pf(lim); // prime factor
for (int i = 2; i < lim; ++i) {
if (pf[i] == 0) {
primes.push_back(i);
pf[i] = i;
}
for (auto p : primes) {
if (i * p >= lim)
break;
pf[i * p] = p;
if (i % p == 0)
break;
}
}
vector<int> mu(lim);
psm.resize(lim);
mu[1] = psm[1] = 1;
for (int i = 2; i < lim; ++i) {
if (i % (pf[i] * pf[i]) != 0)
mu[i] = -mu[i / pf[i]];
psm[i] = psm[i - 1] + mu[i];
}
}
int main() {
ll n;
cin >> n;
while ((__int128)lim * lim * lim <= (__int128)n * n)
++lim;
LinearSieve();
vector<ll> v;
for (ll i = 1, r = 0; i <= n; i = r + 1) {
r = n / (n / i);
v.push_back(i);
}
reverse(v.begin(), v.end());
for (auto i : v) {
if (n / i < lim)
continue;
ll m = n / i;
s[m] = 1;
for (ll j = 2, r = 0; j <= m; j = r + 1) {
r = m / (m / j);
if (m / j < lim)
s[m] -= (r - j + 1) * psm[m / j];
else
s[m] -= (r - j + 1) * s[m / j];
}
}
if (n == 1)
cout << psm[1] << "\n";
else
cout << s[n] << "\n";
}
接著來看第二道例題:
求 \( \displaystyle S(n) = \sum_{i=1}^n \varphi(i) \) 的值。(\(1\le n\le 10^{10}\))
同樣的,要使用杜教篩,就必須找到一個 \( g \) 使得 \( g\ast\varphi \) 和 \( g \) 本身的前綴都很容易計算。
回想數論函數章節所講的引理,選擇 \( 1 \) 當作 \( g \) 的話,可以得到:
- \( 1\ast\varphi = id \),\( \displaystyle\sum_{i=1}^n id(i)=\frac{n\times (n+1)}{2} \)。
- \( \displaystyle\sum_{i=1}^n 1(i) = n \)。
兩者的確都很容易得到。
代入杜教篩的定理中,可以得到:
\[ \displaystyle S(n) = \frac{n\times (n+1)}{2} - \sum_{i=2}^n S\left(\left\lfloor\frac{n}{i}\right\rfloor\right) \]
寫成程式碼如下:
#include <bits/stdc++.h>
using namespace std;
using ll = long long;
ll n;
int lim; // n^(2/3)
vector<__int128> psp; // prefix sum of phi
void LinearSieve() {
vector<int> primes;
vector<int> pf(lim); // prime factor
for (int i = 2; i < lim; ++i) {
if (pf[i] == 0) {
primes.push_back(i);
pf[i] = i;
}
for (auto p : primes) {
if (i * p >= lim)
break;
pf[i * p] = p;
if (i % p == 0)
break;
}
}
vector<int> phi(lim);
psp.resize(lim);
phi[1] = psp[1] = 1;
for (int i = 2; i < lim; ++i) {
int now = 1, res = i;
while (res % pf[i] == 0) {
res /= pf[i];
now *= pf[i];
}
phi[i] = phi[res] * now / pf[i] * (pf[i] - 1);
psp[i] = psp[i - 1] + phi[i];
}
}
unordered_map<ll, __int128> s;
__int128 Find_S(ll n) {
if (n < lim)
return psp[n];
if (s.find(n) != s.end())
return s[n];
s[n] = (__int128)n * (n + 1) / 2;
for (ll i = 2, r = 0; i <= n; i = r + 1) {
r = n / (n / i);
s[n] -= (r - i + 1) * Find_S(n / i);
}
return s[n];
}
int main() {
cin >> n;
while ((__int128)lim * lim * lim <= (__int128)n * n)
++lim;
LinearSieve();
__int128 x = Find_S(n);
string ans = "";
while (x)
ans.push_back(x % 10 + '0'), x /= 10;
reverse(ans.begin(), ans.end());
cout << ans << "\n";
}
或者寫成迭代的樣子:
#include <bits/stdc++.h>
using namespace std;
using ll = long long;
ll n;
int lim; // n^(2/3)
unordered_map<ll, __int128> s; // prefix sum of phi
vector<__int128> psp; // prefix sum of phi
void LinearSieve() {
vector<int> primes;
vector<int> pf(lim); // prime factor
for (int i = 2; i < lim; ++i) {
if (pf[i] == 0) {
primes.push_back(i);
pf[i] = i;
}
for (auto p : primes) {
if (i * p >= lim)
break;
pf[i * p] = p;
if (i % p == 0)
break;
}
}
vector<int> phi(lim);
psp.resize(lim);
phi[1] = psp[1] = 1;
for (int i = 2; i < lim; ++i) {
int now = 1, res = i;
while (res % pf[i] == 0) {
res /= pf[i];
now *= pf[i];
}
phi[i] = phi[res] * now / pf[i] * (pf[i] - 1);
psp[i] = psp[i - 1] + phi[i];
}
}
int main() {
ll n;
cin >> n;
while ((__int128)lim * lim * lim <= (__int128)n * n)
++lim;
LinearSieve();
vector<ll> v;
for (ll i = 1, r = 0; i <= n; i = r + 1) {
r = n / (n / i);
v.push_back(i);
}
reverse(v.begin(), v.end());
for (auto i : v) {
if (n / i < lim)
continue;
ll m = n / i;
s[m] = (__int128)m * (m + 1) / 2;
for (ll j = 2, r = 0; j <= m; j = r + 1) {
r = m / (m / j);
if (m / j < lim)
s[m] -= (r - j + 1) * psp[m / j];
else
s[m] -= (r - j + 1) * s[m / j];
}
}
__int128 x;
if (n == 1)
x = psp[1];
else
x = s[n];
string ans = "";
while (x)
ans.push_back(x % 10 + '0'), x /= 10;
reverse(ans.begin(), ans.end());
cout << ans << "\n";
}
另解,注意到:
\[ \sum_{i=1}^n \phi(i) = \sum_{i=1}^n\sum_{j=1}^i [\gcd(i, j)=1]= \frac{1}{2}\times(\sum_{i=1}^n\sum_{j=1}^n [\gcd(i, j)=1]+1) \]
第一個等式是根據 \( \varphi \) 的定義得到。
第二個等式是因為 \( (a, b) \)、\( (b, a) \) 這樣的 pair 會被重複算到兩次,加一是因為 \( i=1 \)、\( j=1 \) 的情況。
而對於右式:
\begin{align} \sum_{i=1}^n\sum_{j=1}^n [\gcd(i, j)=1] &= \sum_{i=1}^n\sum_{j=1}^n\sum_{d\mid \gcd(i, j)} \mu(d) & \text{套用在數論函數章節講過的引理} \\ &= \sum_{i=1}^n\sum_{j=1}^n\sum_{d\mid i \\ d\mid j} \mu(d) & \text{將 } \gcd \text{ 按照定義拆開} \\ &= \sum_{d=1}^n\sum_{i=1}^{\lfloor\frac n d\rfloor}\sum_{j=1}^{\lfloor\frac n d\rfloor} \mu(d) & \text{套路,交換 sigma 的順序} \\ &= \sum_{d=1}^n \mu(d) \sum_{i=1}^{\lfloor\frac n d\rfloor}\sum_{j=1}^{\lfloor\frac n d\rfloor} 1 & \mu(d) \text{ 和 } i \text{、} j \text{ 無關,移到前面}\\ &= \sum_{d=1}^n \mu(d) \left\lfloor\frac n d\right\rfloor^2 & \text{將後面兩個 sigma 整理好} \end{align}
這是典型數論分塊的形式,為此我們需要求 \( \mu \) 的前綴和。
回想數論分塊:
for (int i = 1, r = 0; i <= n; i = r + 1) {
r = n / (n / i);
ans += 1LL * (S[r] - S[i - 1]) * (n / i);
}
可以知道,要求的前綴和,會是每一塊的右界(\( i - 1 \) 是上一塊的右界),也就是 \( \sum_{i=1}^r \mu(i) \)。
注意到,\( r = \left\lfloor\frac{n}{\left\lfloor n / i \right\rfloor}\right\rfloor \),將 \( \lfloor n / i \rfloor \) 看成另一個整數 \( k \)。
這樣一來,\( r = \left\lfloor\frac n k \right\rfloor \),因此 \( r\in R(n) \)。
也就是說,要求的 \( \sum_{i=1}^r \mu(i) \),恰好和杜教篩中會求到的相同。
因此可以結合例題一還有數論分塊,來求出 \( \varphi \) 的前綴和。
Bottleneck 在於求 \( \mu \) 的前綴和,時間複雜度依然是 \( \mathcal{O}(n^\frac{2}{3}) \)。
洛谷 P213 即是上面兩道練習的模板題。
在此附上 AC code。
在這道題目中要注意,\( n + 1 \) 可能會導致 overflow。
附註:我只有測試過以下的 code,上面的 code 並沒有丟到 OJ 上測試過。
#include <bits/stdc++.h>
using namespace std;
using ll = long long;
int lim = (1 << 21);
int primeFactor[(1 << 21)];
vector<int> primes;
int mu[(1 << 21)];
ll phi[(1 << 21)];
unordered_map<int, int> mu_s;
unordered_map<int, ll> phi_s;
signed main() {
ios::sync_with_stdio(0);
cin.tie(0);
for (int i = 2; i < lim; ++i) {
if (!primeFactor[i]) {
primes.push_back(i);
primeFactor[i] = i;
}
for (auto p : primes) {
if (i * p >= lim)
break;
primeFactor[i * p] = p;
if (i % p == 0)
break;
}
}
mu[1] = 1;
phi[1] = 1;
for (int i = 2; i < lim; ++i) {
if (i % (1LL * primeFactor[i] * primeFactor[i]) == 0) {
mu[i] = 0;
phi[i] = phi[i / primeFactor[i]] * primeFactor[i];
} else {
mu[i] = -mu[i / primeFactor[i]];
phi[i] = phi[i / primeFactor[i]] * (primeFactor[i] - 1);
}
}
for (int i = 2; i < lim; ++i) {
mu[i] += mu[i - 1];
phi[i] += phi[i - 1];
}
int T;
cin >> T;
while (T--) {
int n;
cin >> n;
if (n < lim) {
cout << phi[n] << " " << mu[n] << "\n";
continue;
}
vector<int> R;
for (ll i = 1, r = 0; lim <= n / i; i = r + 1) {
r = n / (n / i);
R.push_back(i);
}
reverse(R.begin(), R.end());
for (auto i : R) {
ll m = n / i;
if (mu_s.find(m) != mu_s.end())
continue;
mu_s[m] = 1;
phi_s[m] = m * (m + 1) / 2LL;
for (ll j = 2, r = 0; j <= m; j = r + 1) {
r = m / (m / j);
if (m / j < lim) {
mu_s[m] -= (r - j + 1) * mu[m / j];
phi_s[m] -= (r - j + 1) * phi[m / j];
} else {
mu_s[m] -= (r - j + 1) * mu_s[m / j];
phi_s[m] -= (r - j + 1) * phi_s[m / j];
}
}
}
cout << phi_s[n] << " " << mu_s[n] << "\n";
}
}
練習題
求\( \displaystyle S(n)=\sum_{i=1}^n \sum_{j=1}^n ij\times \gcd(i,j) \) 的值。(\( n\le 10^{10} \))
(hint:先以數論函數章節所教的技巧進行化簡,變成前綴和的樣子後再做杜教篩。)