Du's Sieve

The Du Sieve is used to handle a class of prefix sum problems for arithmetic functions. For an arithmetic function \(f\), the Du Sieve can compute \(S(n)=\sum_{i=1}^{n}f(i)\) in sublinear time complexity.

Algorithm Idea

We aim to construct a recurrence relation for \(S(n)\) in terms of \(S\left(\left\lfloor\frac{n}{i}\right\rfloor\right)\).

For any arithmetic function \(g\), the following must hold:

\[ \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} \]

where \(f*g\) is the Dirichlet convolution of the arithmetic functions \(f\) and \(g\).

Brief Proof

\(g(d)f\left(\frac{i}{d}\right)\) contributes to all \(i\leq n\), so we change the order of summation, enumerating \(d\) and \(\frac{i}{d}\) (corresponding to new \(i,j\) respectively):

\[ \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} \]

Then we can obtain the recurrence:

\[ \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} \]

If we can construct a suitable arithmetic function \(g\) such that:

\(\sum_{i=1}^n(f * g)(i)\) can be computed quickly;
The prefix sum of \(g\) can be computed quickly, to solve \(\sum_{i=2}^ng(i)S\left(\left\lfloor\dfrac{n}{i}\right\rfloor\right)\) using number-theoretic blocking.

Then we can compute \(g(1)S(n)\) in relatively short time.

Note

Regardless of whether the arithmetic function \(f\) is multiplicative, as long as a suitable arithmetic function \(g\) can be constructed, the Du Sieve can be considered for computing the prefix sum of \(f\).

For example, consider \(f(n)=\mathrm{i}\varphi(n)\). Clearly \(f\) is not multiplicative, but we can take \(g(n)=1\), so that:

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

Computing \(\sum_{k\leq m} (f*g)(k)\) and \(\sum_{k \leq m} g(k)\) both have time complexity \(O(1)\), so the Du Sieve can be considered.

Time Complexity

Let \(R(n)=\left\{\left\lfloor \dfrac{n}{k} \right\rfloor: k=2,3,\dots,n\right\}\). Using the properties of number-theoretic blocking, we know that for any \(m\in R(n)\), \(R(m)\subseteq R(n)\). That is, after using memoization, we only need to compute \(S(k)\) once for all \(k\in R(n)\) to obtain the values for \(R(n)\). The number of such points is \(|R(n)|=O(\sqrt{n})\).

Assume computing \(\sum_{i=1}^n(f * g)(i)\) and \(\sum_{i=1}^n g(i)\) both have time complexity \(O(1)\). Let the time complexity of computing \(S(n)\) be \(T(n)\). Then:

\[ \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} \]

If we can precompute a portion of \(S(k)\), where \(k=1,2,\dots,m\), \(m\geq \lfloor\sqrt n\rfloor\). Let the precomputation time complexity be \(T_0(m)\). Then \(T(n)\) becomes:

\[ \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} \]

If \(T_0(m)=O(m)\) (e.g., linear sieve), by the AM-GM inequality, when \(m=\Theta\left(n^{2/3}\right)\), \(T(n)\) achieves its minimum \(O\left(n^{2/3}\right)\).

An Example of Flawed Proof

Let the complexity of computing \(S(n)\) be \(T(n)\). Then:

\[ 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} \]

Here, \(O\left(\sum_{j=2}^{\lfloor\sqrt{n/i}\rfloor} T\left(\left\lfloor\dfrac{n}{ij}\right\rfloor\right)\right)\) is treated as a higher-order infinitesimal and thus can be discarded. Hence:

\[ \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

The issue lies in the step "treated as a higher-order infinitesimal and thus can be discarded". Substituting \(T\left(\left\lfloor\dfrac{n}{i}\right\rfloor\right)\) into the expression for \(T(n)\) gives:

\[ \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} \]

Consider the part \(\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)\). It is not hard to see:

\[ \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} \]

Since memoization is not introduced, \(T\left(\left\lfloor\sqrt{\dfrac{n}{i}}\right\rfloor\right)\) in the above expression is still \(\Omega\left(\left(\dfrac{n}{i}\right)^{1/4}\right)\). Therefore, the so-called "higher-order infinitesimal" part cannot be discarded.

In fact, the sublinear time complexity of the Du Sieve is guaranteed by memoization. Only after using memoization can we ensure that the multiple summation term does not appear.

Examples

Problem 1

P4213【Template】Du Sieve (Sum)

Compute \(S_1(n)= \sum_{i=1}^{n} \mu(i)\) and \(S_2(n)= \sum_{i=1}^{n} \varphi(i)\), where \(1\leq n<2^{31}\).

Prefix Sum of Möbius FunctionPrefix Sum of Euler's Totient Function

We know:

\[ \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} \]

The derivation of time complexity is in the Time Complexity section.

For larger values, a map/unordered_map should be used to store the computed results for future reuse.

The Du Sieve can also be used to find the prefix sum of \(\varphi (x)\), but a better method is to apply Möbius inversion.

Möbius InversionDu Sieve

\[ \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} \]

Since the problem asks for \(\sum_{i=1}^n \sum_{j=1}^i [\gcd(i,j)=1]\), we exclude the case \(i=1,j=1\) and divide the result by \(2\).

Observe that we only need the prefix sum of the Möbius function to quickly compute the prefix sum of Euler's totient function. Time complexity is \(O\left(n^{\frac 2 3}\right)\).

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

Similarly, \(\varphi * 1=\operatorname{id}\), so:

\[ \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} \]

Code Implementation

#include <cstring>
#include <iostream>
#include <map>
using namespace std;
constexpr 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) {  // Calculate prefix sum of mu
  if (x < MAXN) return sum_mu[x];
  if (mp_mu[x]) return mp_mu[x];  // If the mu value for this size already exists in the map, return directly
  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;  // Path compression for future calculations
}

long long S_phi(long long x) {  // Calculate prefix sum of 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() {
  cin.tie(nullptr)->sync_with_stdio(false);
  cin >> T;
  mu[1] = 1;
  for (int i = 2; i < MAXN; i++) {  // Linear sieve preprocessing for mu array
    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];  // Calculate prefix sum of mu array
  while (T--) {
    cin >> n;
    cout << S_phi(n) << ' ' << S_mu(n) << '\n';
  }
  return 0;
}

Problem 2

「LuoguP3768」Simple Math Problem

Summary: Compute

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

where \(n\leq 10^{10},5\times 10^8\leq p\leq 1.1\times 10^9\), and \(p\) is a prime.

Using \(\varphi * 1=\operatorname{id}\) for Möbius inversion transforms it into:

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

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

Perform number-theoretic blocking on \(\sum_{d=1}^nF\left(\left\lfloor\dfrac{n}{d}\right\rfloor\right)^2\), and handle the prefix sum of \(d^2\varphi(d)\) using the Du Sieve:

\[ 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) \]

We need to construct a multiplicative function \(g\) such that \(f\times g\) and \(g\) can be summed quickly.

The prefix sum of \(\varphi\) alone can be handled by the Du Sieve with \(\varphi * 1\), but here \(f\) has an extra \(\operatorname{id}^2\), so we convolve with \(\operatorname{id}^2\) to turn it into a constant:

\[ 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) \]

Simplify the convolution:

\[ \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} \]

Further simplify \(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} \]

Solve using blocking.

Code Implementation

// Do not skimp on array size to save memory, it will cause TLE at 80
#include <cmath>
#include <iostream>
#include <map>
using namespace std;
constexpr 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) {  // For calculating modular inverse
  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) {  // Linear sieve for phi, s
  bp[0] = bp[1] = true, 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]] = true;
      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) {  // Sum of cubes
  return k %= P, (k * (k + 1) / 2) % P * ((k * (k + 1) / 2) % P) % P;
}

long long s2(long long k) {  // Sum of squares
  return k %= P, k * (k + 1) % P * (k * 2 + 1) % P * inv6 % P;
}

long long calc(long long k) {  // Calculate S(k)
  if (k <= pn) return s[k];
  if (s_map[k]) return s_map[k];  // For values exceeding pn, use map for discrete storage
  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() {
  cin.tie(nullptr)->sync_with_stdio(false);
  cin >> 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);
  cout << solve();
  return 0;
}

References

Ren Zhizhou, 2016, "Several Methods for Summing Multiplicative Functions", 2016 Candidate Papers for the Chinese National Team of the Informatics Olympiad.
Time and Space Complexity Analysis of the Du Sieve - riteme.site

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