杜教篩

杜教篩被用於處理一類數論函數的前綴和問題。對於數論函數 \(f\)，杜教篩可以在低於線性時間的複雜度內計算 \(S(n)=\sum_{i=1}^{n}f(i)\)。

算法思想

我們想辦法構造一個 \(S(n)\) 關於 \(S\left(\left\lfloor\frac{n}{i}\right\rfloor\right)\) 的遞推式。

對於任意一個數論函數 \(g\)，必滿足：

\[ \begin{aligned} \sum_{i=1}^{n}(f * g)(i) & =\sum_{i=1}^{n}\sum_{d \mid i}g(d)f\left(\frac{i}{d}\right) \\ & =\sum_{i=1}^{n}g(i)S\left(\left\lfloor\frac{n}{i}\right\rfloor\right) \end{aligned} \]

其中 \(f*g\) 為數論函數 \(f\) 和 \(g\) 的狄利克雷卷積。

略證

\(g(d)f\left(\frac{i}{d}\right)\) 就是對所有 \(i\leq n\) 的做貢獻，因此變換枚舉順序，枚舉 \(d\),\(\frac{i}{d}\)（分別對應新的 \(i,j\)）

\[ \begin{aligned} \sum_{i=1}^n\sum_{d \mid i}g(d)f\left(\frac{i}{d}\right) & =\sum_{i=1}^n\sum_{j=1}^{\left\lfloor n/i \right\rfloor}g(i)f(j) \\ & =\sum_{i=1}^ng(i)\sum_{j=1}^{\left\lfloor n/i \right\rfloor}f(j) \\ & =\sum_{i=1}^ng(i)S\left(\left\lfloor\frac{n}{i}\right\rfloor\right) \end{aligned} \]

那麼可以得到遞推式：

\[ \begin{aligned} g(1)S(n) & = \sum_{i=1}^n g(i)S\left(\left\lfloor\frac{n}{i}\right\rfloor\right) - \sum_{i=2}^n g(i)S\left(\left\lfloor\frac{n}{i}\right\rfloor\right) \\ & = \sum_{i=1}^n (f * g)(i) - \sum_{i=2}^n g(i)S\left(\left\lfloor\frac{n}{i}\right\rfloor\right) \end{aligned} \]

假如我們可以構造恰當的數論函數 \(g\) 使得：

可以快速計算 \(\sum_{i=1}^n(f * g)(i)\)；
可以快速計算 \(g\) 的前綴和，以用數論分塊求解 \(\sum_{i=2}^ng(i)S\left(\left\lfloor\dfrac{n}{i}\right\rfloor\right)\)。

則我們可以在較短時間內求得 \(g(1)S(n)\)。

注意

無論數論函數 \(f\) 是否為積性函數，只要可以構造出恰當的數論函數 \(g\), 便都可以考慮用杜教篩求 \(f\) 的前綴和。

如考慮 \(f(n)=\mathrm{i}\varphi(n)\), 顯然 \(f\) 不是積性函數，但可取 \(g(n)=1\), 從而：

\[ \sum_{k=1}^n (f*g)(k)=\mathrm{i}\frac{n(n+1)}{2} \]

計算 \(\sum_{k\leq m} (f*g)(k)\) 和 \(\sum_{k \leq m} g(k)\) 的時間複雜度均為 \(O(1)\), 故可以考慮使用杜教篩。

時間複雜度

令 \(R(n)=\left\{\left\lfloor \dfrac{n}{k} \right\rfloor: k=2,3,\dots,n\right\}\), 我們有如下引理：

引理

對任意的 \(m\in R(n)\)，我們有 \(R(m)\subseteq R(n)\).

證明

令 \(m=\left\lfloor\dfrac{n}{x}\right\rfloor\), 任取 \(\left\lfloor\dfrac{m}{y}\right\rfloor \in R(m)\), 由整除分塊/數論分塊的引理 1 可知：

\[ \left\lfloor\dfrac{m}{y}\right\rfloor=\left\lfloor\dfrac{n}{xy}\right\rfloor\in R(n). \]

設計算 \(\sum_{i=1}^n(f * g)(i)\) 和 \(\sum_{i=1}^n g(i)\) 的時間複雜度均為 \(O(1)\). 由引理可知：使用記憶化之後，每個 \(S(k)~(k\in R(n))\) 均只會計算一次。

由整除分塊/數論分塊的引理 2 可知 \(|R(n)|=2\sqrt{n}+\Theta(1)\). 設計算 \(S(n)\) 的時間複雜度為 \(T(n)\), 則：

\[ \begin{aligned} T(n) & = \sum_{k\in R(n)} T(k)\\ & = \Theta(\sqrt n)+\sum_{k=1}^{\lfloor\sqrt n\rfloor} O(\sqrt k)+\sum_{k=2}^{\lfloor\sqrt n\rfloor} O\left(\sqrt{\dfrac{n}{k}}\right)\\ & = O\left(\int_{0}^{\sqrt n} \left(\sqrt{x} + \sqrt{\dfrac{n}{x}}\right) \mathrm{d}x\right)\\ & = O\left(n^{3/4}\right). \end{aligned} \]

若我們可以預處理出一部分 \(S(k)\), 其中 \(k=1,2,\dots,m\)，\(m\geq \lfloor\sqrt n\rfloor\). 設預處理的時間複雜度為 \(T_0(m)\), 則此時的 \(T(n)\) 為：

\[ \begin{aligned} T(n) & = T_0(m)+\sum_{k\in R(n);k>m} T(k)\\ & = T_0(m)+\sum_{k=1}^{\lfloor n/m \rfloor} O\left(\sqrt{\dfrac{n}{k}}\right)\\ & = O\left(T_0(m)+\int_{0}^{n/m} \sqrt{\dfrac{n}{x}} \mathrm{d}x\right)\\ & = O\left(T_0(m)+\dfrac{n}{\sqrt m}\right). \end{aligned} \]

若 \(T_0(m)=O(m)\)（如線性篩），由均值不等式可知：當 \(m=\Theta\left(n^{2/3}\right)\) 時，\(T(n)\) 取得最小值 \(O\left(n^{2/3}\right)\).

偽證一例

設計算 \(S(n)\) 的複雜度為 \(T(n)\), 則有：

\[ T(n)=\Theta\left(\sqrt{n}\right)+O\left(\sum_{i=2}^{\lfloor\sqrt{n}\rfloor} T\left(\left\lfloor\frac{n}{i}\right\rfloor\right)\right) \]

\[ \begin{aligned} T\left(\left\lfloor\frac{n}{i}\right\rfloor\right) & = \Theta\left(\sqrt{\frac{n}{i}}\right)+O\left(\sum_{j=2}^{\lfloor\sqrt{n/i}\rfloor} T\left(\left\lfloor\frac{n}{ij}\right\rfloor\right)\right) \\ & = O\left(\sqrt{\frac{n}{i}}\right) \end{aligned} \]

Note

\(O\left(\sum_{j=2}^{\lfloor\sqrt{n/i}\rfloor} T\left(\left\lfloor\dfrac{n}{ij}\right\rfloor\right)\right)\) 視作高階無窮小，從而可以捨去。

故：

\[ \begin{aligned} T(n) & = \Theta\left(\sqrt{n}\right)+O\left(\sum_{i=2}^{\lfloor\sqrt{n}\rfloor} \sqrt{\frac{n}{i}}\right) \\ & = O\left(\sum_{i=1}^{\lfloor\sqrt{n}\rfloor} \sqrt{\frac{n}{i}}\right) \\ & = O\left(\int_{0}^{\sqrt{n}}\sqrt{\frac{n}{x}}\mathrm{d}x\right) \\ & = O\left(n^{3/4}\right) \end{aligned} \]

Bug

問題在於「視作高階無窮小，從而可以捨去」這一處。我們將 \(T\left(\left\lfloor\dfrac{n}{i}\right\rfloor\right)\) 代入 \(T(n)\) 的式子裏，有：

\[ \begin{aligned} T(n) & = \Theta\left(\sqrt{n}\right)+O\left(\sum_{i=2}^{\lfloor\sqrt{n}\rfloor} \sqrt{\frac{n}{i}}\right)+O\left(\sum_{i=2}^{\lfloor\sqrt{n}\rfloor}\sum_{j=2}^{\lfloor\sqrt{n/i}\rfloor} T\left(\left\lfloor\frac{n}{ij}\right\rfloor\right)\right)\\ & = O\left(\sqrt{n}+\int_{0}^{\sqrt{n}}\sqrt{\frac{n}{x}}\mathrm{d}x\right)+O\left(\sum_{i=2}^{\lfloor\sqrt{n}\rfloor}\sum_{j=2}^{\lfloor\sqrt{n/i}\rfloor} T\left(\left\lfloor\frac{n}{ij}\right\rfloor\right)\right)\\ & = O\left(n^{3/4}\right)+O\left(\sum_{i=2}^{\lfloor\sqrt{n}\rfloor}\sum_{j=2}^{\lfloor\sqrt{n/i}\rfloor} T\left(\left\lfloor\frac{n}{ij}\right\rfloor\right)\right)\\ \end{aligned} \]

我們考慮 \(\displaystyle\sum_{i=2}^{\lfloor\sqrt{n}\rfloor}\sum_{j=2}^{\lfloor\sqrt{n/i}\rfloor} T\left(\left\lfloor\frac{n}{ij}\right\rfloor\right)\) 這部分，不難發現：

\[ \begin{aligned} \sum_{i=2}^{\lfloor\sqrt{n}\rfloor}\sum_{j=2}^{\lfloor\sqrt{n/i}\rfloor} T\left(\left\lfloor\frac{n}{ij}\right\rfloor\right) & = \Omega\left(\sum_{i=2}^{\lfloor\sqrt{n}\rfloor} T\left(\left\lfloor\frac{n}{i}\cdot\left\lfloor\sqrt\frac{n}{i}\right\rfloor^{-1}\right\rfloor\right)\right) \\ & = \Omega\left(\sum_{i=2}^{\lfloor\sqrt{n}\rfloor} T\left(\left\lfloor\sqrt\frac{n}{i}\right\rfloor\right)\right) \end{aligned} \]

由於沒有引入記憶化，因此上式中的 \(T\left(\left\lfloor\sqrt{\dfrac{n}{i}}\right\rfloor\right)\) 仍然是 \(\Omega\left(\left(\dfrac{n}{i}\right)^{1/4}\right)\) 的，進而所謂的「高階無窮小」部分是不可以捨去的。

實際上杜教篩的亞線性時間複雜度是由記憶化保證的。只有使用了記憶化之後才能保證不會出現那個多重求和的項。

例題

問題一

P4213【模板】杜教篩（Sum）

求 \(S_1(n)= \sum_{i=1}^{n} \mu(i)\) 和 \(S_2(n)= \sum_{i=1}^{n} \varphi(i)\) 的值，\(1\leq n<2^{31}\).

莫比烏斯函數前綴和歐拉函數前綴和

我們知道：

\[ \epsilon = [n=1] = \mu * 1 = \sum_{d \mid n} \mu(d) \]

\[ \begin{aligned} S_1(n) & =\sum_{i=1}^n \epsilon (i)-\sum_{i=2}^n S_1 \left(\left\lfloor \frac n i \right\rfloor\right) \\ & = 1-\sum_{i=2}^n S_1\left(\left\lfloor \frac n i \right\rfloor\right) \end{aligned} \]

時間複雜度的推導見時間複雜度一節。

對於較大的值，需要用 map/unordered_map 存下其對應的值，方便以後使用時直接使用之前計算的結果。

當然也可以用杜教篩求出 \(\varphi (x)\) 的前綴和，但是更好的方法是應用莫比烏斯反演。

莫比烏斯反演杜教篩

\[ \begin{aligned} \sum_{i=1}^n \sum_{j=1}^n [\gcd(i,j)=1] & =\sum_{i=1}^n \sum_{j=1}^n \sum_{d \mid i,d \mid j} \mu(d) \\ & =\sum_{d=1}^n \mu(d) {\left\lfloor \frac n d \right\rfloor}^2 \end{aligned} \]

由於題目所求的是 \(\sum_{i=1}^n \sum_{j=1}^i [\gcd(i,j)=1]\), 所以我們排除掉 \(i=1,j=1\) 的情況，並將結果除以 \(2\) 即可。

觀察到，只需求出莫比烏斯函數的前綴和，就可以快速計算出歐拉函數的前綴和了。時間複雜度 \(O\left(n^{\frac 2 3}\right)\).

求 \(S(n)=\sum_{i=1}^n\varphi(i)\).

同樣的，\(\varphi * 1=\operatorname{id}\), 從而：

\[ \begin{aligned} S(n) & =\sum_{i=1}^n i - \sum_{i=2}^n S\left(\left\lfloor\frac{n}{i}\right\rfloor\right) \\ & =\frac{1}{2}n(n+1) - \sum_{i=2}^n S\left(\left\lfloor\frac{n}{i}\right\rfloor\right) \end{aligned} \]

代碼實現

#include <algorithm>
#include <cstdio>
#include <cstring>
#include <map>
using namespace std;
const int maxn = 2000010;
long long T, n, pri[maxn], cur, mu[maxn], sum_mu[maxn];
bool vis[maxn];
map<long long, long long> mp_mu;

long long S_mu(long long x) {  // 求mu的前缀和
  if (x < maxn) return sum_mu[x];
  if (mp_mu[x]) return mp_mu[x];  // 如果map中已有该大小的mu值，则可直接返回
  long long ret = (long long)1;
  for (long long i = 2, j; i <= x; i = j + 1) {
    j = x / (x / i);
    ret -= S_mu(x / i) * (j - i + 1);
  }
  return mp_mu[x] = ret;  // 路径压缩，方便下次计算
}

long long S_phi(long long x) {  // 求phi的前缀和
  long long ret = (long long)0;
  long long j;
  for (long long i = 1; i <= x; i = j + 1) {
    j = x / (x / i);
    ret += (S_mu(j) - S_mu(i - 1)) * (x / i) * (x / i);
  }
  return (ret - 1) / 2 + 1;
}

int main() {
  scanf("%lld", &T);
  mu[1] = 1;
  for (int i = 2; i < maxn; i++) {  // 线性筛预处理mu数组
    if (!vis[i]) {
      pri[++cur] = i;
      mu[i] = -1;
    }
    for (int j = 1; j <= cur && i * pri[j] < maxn; j++) {
      vis[i * pri[j]] = true;
      if (i % pri[j])
        mu[i * pri[j]] = -mu[i];
      else {
        mu[i * pri[j]] = 0;
        break;
      }
    }
  }
  for (int i = 1; i < maxn; i++)
    sum_mu[i] = sum_mu[i - 1] + mu[i];  // 求mu数组前缀和
  while (T--) {
    scanf("%lld", &n);
    printf("%lld %lld\n", S_phi(n), S_mu(n));
  }
  return 0;
}

問題二

「LuoguP3768」簡單的數學題

大意：求

\[ \sum_{i=1}^n\sum_{j=1}^ni\cdot j\cdot\gcd(i,j)\pmod p \]

其中 \(n\leq 10^{10},5\times 10^8\leq p\leq 1.1\times 10^9\),\(p\) 是質數。

利用 \(\varphi * 1=\operatorname{id}\) 做莫比烏斯反演化為：

\[ \sum_{d=1}^nF^2\left(\left\lfloor\frac{n}{d}\right\rfloor\right)\cdot d^2\varphi(d) \]

其中 \(F(n)=\dfrac{1}{2}n(n+1)\)

對 \(\sum_{d=1}^nF\left(\left\lfloor\dfrac{n}{d}\right\rfloor\right)^2\) 做數論分塊，\(d^2\varphi(d)\) 的前綴和用杜教篩處理：

\[ f(n)=n^2\varphi(n)=(\operatorname{id}^2\varphi)(n) \]

\[ S(n)=\sum_{i=1}^nf(i)=\sum_{i=1}^n(\operatorname{id}^2\varphi)(i) \]

需要構造積性函數 \(g\)，使得 \(f\times g\) 和 \(g\) 能快速求和。

單純的 \(\varphi\) 的前綴和可以用 \(\varphi * 1\) 的杜教篩處理，但是這裏的 \(f\) 多了一個 \(\operatorname{id}^2\)，那麼我們就卷一個 \(\operatorname{id}^2\) 上去，讓它變成常數：

\[ S(n)=\sum_{i=1}^n\left(\left(\operatorname{id}^2\varphi\right) * \operatorname{id}^2\right)(i)-\sum_{i=2}^n\operatorname{id}^2(i)S\left(\left\lfloor\frac{n}{i}\right\rfloor\right) \]

化一下卷積：

\[ \begin{aligned} ((\operatorname{id}^2\varphi)* \operatorname{id}^2)(i) & =\sum_{d \mid i}\left(\operatorname{id}^2\varphi\right)(d)\operatorname{id}^2\left(\frac{i}{d}\right) \\ & =\sum_{d \mid i}d^2\varphi(d)\left(\frac{i}{d}\right)^2 \\ & =\sum_{d \mid i}i^2\varphi(d)=i^2\sum_{d \mid i}\varphi(d) \\ & =i^2(\varphi*1)(i)=i^3 \end{aligned} \]

再化一下 \(S(n)\):

\[ \begin{aligned} S(n) & =\sum_{i=1}^n\left((\operatorname{id}^2\varphi)* \operatorname{id}^2\right)(i)-\sum_{i=2}^n\operatorname{id}^2(i)S\left(\left\lfloor\frac{n}{i}\right\rfloor\right) \\ & =\sum_{i=1}^ni^3-\sum_{i=2}^ni^2S\left(\left\lfloor\frac{n}{i}\right\rfloor\right) \\ & =\left(\frac{1}{2}n(n+1)\right)^2-\sum_{i=2}^ni^2S\left(\left\lfloor\frac{n}{i}\right\rfloor\right) \\ \end{aligned} \]

分塊求解即可。

代碼實現

#include <cmath>
#include <cstdio>
#include <map>
using namespace std;
const int N = 5e6, NP = 5e6, SZ = N;
long long n, P, inv2, inv6, s[N];
int phi[N], p[NP], cnt, pn;
bool bp[N];
map<long long, long long> s_map;

long long ksm(long long a, long long m) {  // 求逆元用
  long long res = 1;
  while (m) {
    if (m & 1) res = res * a % P;
    a = a * a % P, m >>= 1;
  }
  return res;
}

void prime_work(int k) {  // 线性筛phi，s
  bp[0] = bp[1] = 1, phi[1] = 1;
  for (int i = 2; i <= k; i++) {
    if (!bp[i]) p[++cnt] = i, phi[i] = i - 1;
    for (int j = 1; j <= cnt && i * p[j] <= k; j++) {
      bp[i * p[j]] = 1;
      if (i % p[j] == 0) {
        phi[i * p[j]] = phi[i] * p[j];
        break;
      } else
        phi[i * p[j]] = phi[i] * phi[p[j]];
    }
  }
  for (int i = 1; i <= k; i++)
    s[i] = (1ll * i * i % P * phi[i] % P + s[i - 1]) % P;
}

long long s3(long long k) {  // 立方和
  return k %= P, (k * (k + 1) / 2) % P * ((k * (k + 1) / 2) % P) % P;
}

long long s2(long long k) {  // 平方和
  return k %= P, k * (k + 1) % P * (k * 2 + 1) % P * inv6 % P;
}

long long calc(long long k) {  // 计算S(k)
  if (k <= pn) return s[k];
  if (s_map[k]) return s_map[k];  // 对于超过pn的用map离散存储
  long long res = s3(k), pre = 1, cur;
  for (long long i = 2, j; i <= k; i = j + 1)
    j = k / (k / i), cur = s2(j),
    res = (res - calc(k / i) * (cur - pre) % P) % P, pre = cur;
  return s_map[k] = (res + P) % P;
}

long long solve() {
  long long res = 0, pre = 0, cur;
  for (long long i = 1, j; i <= n; i = j + 1) {
    j = n / (n / i);
    cur = calc(j);
    res = (res + (s3(n / i) * (cur - pre)) % P) % P;
    pre = cur;
  }
  return (res + P) % P;
}

int main() {
  scanf("%lld%lld", &P, &n);
  inv2 = ksm(2, P - 2), inv6 = ksm(6, P - 2);
  pn = (long long)pow(n, 0.666667);  // n^(2/3)
  prime_work(pn);
  printf("%lld", solve());
  return 0;
}  // 不要为了省什么内存把数组开小,会卡80

參考資料

任之洲，2016，《積性函數求和的幾種方法》，2016 年信息學奧林匹克中國國家隊候選隊員論文
杜教篩的時空複雜度分析 - riteme.site

本页面最近更新：，更新历史
发现错误？想一起完善？在 GitHub 上编辑此页！
本页面贡献者：hsfzLZH1, sshwy, StudyingFather, Marcythm
本页面的全部内容在 CC BY-SA 4.0 和 SATA 协议之条款下提供，附加条款亦可能应用