AIOJ

The Du sieve is used to compute prefix sums of a class of number-theoretic functions. For a number-theoretic function $f$ , it can compute $S(n)=\sum_{i=1}^{n}f(i)$ in sublinear time.

Algorithmic Idea

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

For any number-theoretic function $g$ , we have

\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 number-theoretic functions $f$ and $g$ .

Sketch of Proof

The term $g(d)f\left(\frac{i}{d}\right)$ contributes to every $i\leq n$ . Therefore, changing the enumeration order and enumerating $d$ and $\frac{i}{d}$ (corresponding to the new $i,j$ ) gives

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

Thus we 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}

Suppose that we can construct a suitable number-theoretic function $g$ such that:

$\sum_{i=1}^n(f * g)(i)$ can be computed quickly;
the prefix sums of $g$ can be computed quickly, so that $\sum_{i=2}^ng(i)S\left(\left\lfloor\dfrac{n}{i}\right\rfloor\right)$ can be solved using quotient grouping.

Then we can obtain $g(1)S(n)$ in a short time.

Note

Regardless of whether the number-theoretic function $f$ is multiplicative, as long as a suitable number-theoretic 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$ , giving

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

The time complexities of computing $\sum_{k\leq m} (f*g)(k)$ and $\sum_{k \leq m} g(k)$ are both $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\}$ . By the property of quotient grouping, for any $m\in R(n)$ , we have $R(m)\subseteq R(n)$ . In other words, with memoization, it is enough to compute $k\in R(n)$ once for every $S(k)$ to obtain the values of $R(n)$ . The number of these points is $|R(n)|=O(\sqrt{n})$ .

Suppose the time complexities of computing $\sum_{i=1}^n(f * g)(i)$ and $\sum_{i=1}^n g(i)$ are both $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 preprocess some values $S(k)$ for $k=1,2,\dots,m$ , where $m\geq \lfloor\sqrt n\rfloor$ , and the preprocessing time complexity is $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)$ , as with a linear sieve, the arithmetic mean-geometric mean inequality shows that $m=\Theta\left(n^{2/3}\right)$ is minimized at $T(n)$ , giving $O\left(n^{2/3}\right)$ .

An Incorrect Proof

Suppose the complexity of computing $S(n)$ is $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 discarded. Therefore,

\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 problem lies in treating the term as a higher-order infinitesimal and discarding it. 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 term $\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 difficult to see that

\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 has not been introduced, $T\left(\left\lfloor\sqrt{\dfrac{n}{i}}\right\rfloor\right)$ in the expression above is still $\Omega\left(\left(\dfrac{n}{i}\right)^{1/4}\right)$ , so the supposedly higher-order infinitesimal term cannot be discarded.

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

Examples

Problem 1

P4213 [Template] Du Sieve (Sum)

Find the values of $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 the Möbius Function

We know that

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

See the section on time complexity for the derivation.

For large values, the corresponding values should be stored in a map/unordered_map, so that previously computed results can be used directly later.

Prefix Sum of the Euler Totient Function

The Du sieve can of course also be used to find the prefix sum of $\varphi (x)$ , but Möbius inversion gives a better method.

Möbius Inversion

\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$ .

It can be seen that once the prefix sum of the Möbius function is known, the prefix sum of the Euler totient function can be computed quickly. The time complexity is $O\left(n^{\frac 2 3}\right)$ .

Du Sieve

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}

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) {  // Compute the prefix sum of mu.
  if (x < MAXN) return sum_mu[x];
  if (mp_mu[x]) return mp_mu[x];  // Return the existing value if mu for this size is already in the map.
  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 makes later computation faster.
}

long long S_phi(long long x) {  // Compute the 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++) {  // Precompute the mu array with a linear sieve.
    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];  // Compute the prefix sum of the mu array.
  while (T--) {
    cin >> n;
    cout << S_phi(n) << ' ' << S_mu(n) << '\n';
  }
  return 0;
}

Problem 2

Luogu P3768, “A Simple Math Problem”

In essence, 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 prime.

Using $\varphi * 1=\operatorname{id}$ and Möbius inversion, this becomes

\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)$ .

Apply quotient grouping to $\sum_{d=1}^nF\left(\left\lfloor\dfrac{n}{d}\right\rfloor\right)^2$ , and use the Du sieve to handle the prefix sum of $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)

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 applying the Du sieve to $\varphi * 1$ , but here $f$ has an additional $\operatorname{id}^2$ . We therefore convolve by an $\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)

Simplifying the convolution gives

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

Simplifying $S(n)$ once more:

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

The problem can then be solved by quotient grouping.

Implementation

// Do not make the arrays smaller just to save memory; it will fail on test 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) {  // Used to compute inverses.
  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 and 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) {  // Compute S(k).
  if (k <= pn) return s[k];
  if (s_map[k]) return s_map[k];  // Store values above pn discretely in the 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() {
  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

Zhizhou Ren, 2016, Several Methods for Summing Multiplicative Functions, a paper by a candidate for the Chinese National Olympiad in Informatics team.
Analysis of the time and space complexity of the Du sieve - riteme.site