杜教篩

例題

我們先來看一道題目。

求 $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：先以數論函數章節所教的技巧進行化簡，變成前綴和的樣子後再做杜教篩。）