Skip to content

Dirichlet Convolution

This article introduces Dirichlet convolution and Dirichlet generating functions.

Dirichlet Convolution

The Dirichlet convolution of arithmetic functions \(f(n)\) and \(g(n)\), denoted as \(f \ast g\), is defined as the arithmetic function

\[ (f \ast g)(n) = \sum_{k\mid n}f(k)g\left(\dfrac{n}{k}\right) = \sum_{k\ell=n}f(k)g(\ell). \]

Dirichlet convolution is an important operation on arithmetic functions. Many properties of arithmetic functions are revealed through this operation.

Examples
  1. The unit function \(\varepsilon\) is the Dirichlet convolution of the Möbius function \(\mu\) and the constant function \(1\):
\[ \varepsilon=\mu \ast 1 \iff\varepsilon(n)=\sum_{d\mid n}\mu(d). \]
  1. The divisor count function \(\tau\) is the Dirichlet convolution of the constant function \(1\) with itself:
\[ \tau=1 \ast 1 \iff \tau(n)=\sum_{d\mid n}1. \]
  1. The divisor sum function \(\sigma\) is the Dirichlet convolution of the identity function \(\mathrm{id}\) and the constant function \(1\):
\[ \sigma=\mathrm{id} \ast 1 \iff\sigma(n)=\sum_{d\mid n}d. \]
  1. Euler's totient function \(\varphi\) is the Dirichlet convolution of the identity function \(\mathrm{id}\) and the Möbius function \(\mu\):
\[ \varphi=\mathrm{id}\ast \mu \iff\varphi(n)=\sum_{d\mid n}d\cdot\mu\left(\frac{n}{d}\right). \]

Möbius inversion transforms identities of arithmetic functions by utilizing \(\varepsilon=\mu \ast 1\).

Properties

Dirichlet convolution possesses a series of algebraic properties.

Theorem

Let \(f,g,h\) be arithmetic functions. Then:

  1. Commutativity: \(f\ast g=g\ast f\).
  2. Associativity: \((f\ast g)\ast h=f\ast(g\ast h)\).
  3. Distributivity: \((f+g)\ast h = f\ast h + g\ast h\).
  4. Identity element: \(f\ast\varepsilon = \varepsilon \ast f = f\), where \(\varepsilon(n) = [n=1]\) is the convolution identity, and \([\cdot]\) is the Iverson bracket.
  5. Inverse: If and only if \(f(1)\neq 0\), there exists \(g\) such that \(f\ast g=g\ast f=\varepsilon\), and \(g\) is called the Dirichlet inverse of \(f\), denoted as \(f^{-1}\). Moreover, the inverse \(g\) satisfies the recurrence formula
\[ g(n) = \dfrac{\varepsilon(n) - \sum_{k\ell = n,~k\neq 1}f(k)g(\ell)}{f(1)}. \]
Proof

To verify commutativity, direct calculation shows

\[ (f\ast g)(n) = \sum_{k\ell=n}f(k)g(\ell) = (g\ast f)(n). \]

To verify associativity, direct calculation shows

\[ ((f\ast g)\ast h)(n) = \sum_{k\ell m = n}f(k)g(\ell)h(m) = (f\ast (g\ast h))(n). \]

To verify distributivity, direct calculation shows

\[ \begin{aligned} ((f+g)\ast h)(n) &= \sum_{k\ell = n}(f(k) + g(k))h(\ell) \\ &= \sum_{k\ell=n}f(k)h(\ell) + \sum_{k\ell=n}g(k)h(\ell) = (f\ast h+g\ast h)(n). \end{aligned} \]

To verify that \(\varepsilon(n)\) is the identity element, direct calculation shows

\[ (f\ast\varepsilon)(n) = \sum_{k\ell = n}f(k)\varepsilon(\ell) = f(n). \]

The second equality holds because \(\varepsilon(\ell)\) is nonzero only when \(\ell=1\), i.e., \(k=n\).

Finally, we need to prove that \(f^{-1}\) exists if and only if \(f(1)\neq 0\). For any \(f\), suppose there exists \(g\) such that \(f\ast g=\varepsilon\). This implies

\[ (f\ast g)(n) = \sum_{k\ell = n}f(k)g(\ell) = \varepsilon(n). \]

This actually gives a system of equations for the values of \(g(n)\), from which \(g(n)\) can be directly solved. In particular, when \(n=1\), the equation becomes \(f(1)g(1)=1\), so the existence of \(g\) requires at least \(f(1)\neq 0\). As long as \(f(1)\neq 0\), we can directly solve

\[ g(n) = \dfrac{\varepsilon(n) - \sum_{k\ell = n,~k\neq 1}f(k)g(\ell)}{f(1)}. \]

This can be used to recursively compute the values of \(g(n)\). Therefore, the inverse \(g\) exists if and only if \(f(1)\neq 0\).

In the language of abstract algebra, these algebraic properties indicate that the set of all arithmetic functions forms a commutative ring under pointwise addition and Dirichlet convolution, and its set of invertible elements consists precisely of those functions that take nonzero values at \(n=1\). This ring is called the Dirichlet ring.

Multiplicative functions are a special class of arithmetic functions. They are closed under Dirichlet convolution and Dirichlet inversion.

Theorem

Let \(f,g\) be multiplicative functions. Then, \(f\ast g\) is also multiplicative. Moreover, the inverse \(f^{-1}\) must exist and is also multiplicative.

Proof

For the first point, let \(h=f\ast g\). Direct verification shows that for \(n_1\perp n_2\), we have

\[ \begin{aligned} h(n_1)h(n_2) &= \left(\sum_{k_1\ell_1=n_1}f(k_1)g(\ell_1)\right)\left(\sum_{k_2\ell_2=n_2}f(k_2)g(\ell_2)\right)\\ &= \sum_{k_1\ell_1=n_1,~k_2\ell_2=n_2}f(k_1)f(k_2)g(\ell_1)g(\ell_2)\\ &= \sum_{k\ell = n_1n_2}f(k)g(\ell) \\ &= h(n_1n_2). \end{aligned} \]

The logic behind the third equality's change of summation order is: when \(k\) runs over the divisors of \(n_1n_2\), the prime factors of \(k\) can be divided into two categories based on whether they belong to \(n_1\) or \(n_2\). Multiplying the prime factors (counting multiplicity) from each category gives \(k_1\) and \(k_2\), which will run over the divisors of \(n_1\) and \(n_2\), respectively. Conversely, from the divisors \(k_1\) and \(k_2\) of \(n_1\) and \(n_2\), we can always obtain the divisor \(k=k_1k_2\) of \(n_1n_2\).

For the second point, let \(g=f^{-1}\) and consider mathematical induction. First, \(g(1)=1/f(1)=1\). At this point, the recurrence formula for the inverse can be written as

\[ g(n) = \varepsilon(n) - \sum_{k\ell = n,~k\neq 1} f(k)g(\ell). \]

Therefore, for \(n_1\perp n_2\) and \(n_1n_2 > 1\), we have

\[ \begin{aligned} g(n_1n_2) &= -\sum_{k\ell=n_1n_2,~k\neq 1}f(k)g(\ell) \\ &= -\sum_{k_1\ell_1=n_1,~k_2\ell_2=n_2,~k_1k_2\neq 1}f(k_1)f(k_2)g(\ell_1)g(\ell_2)\\ &= f(1)f(1)g(n_1)g(n_2) - \sum_{k_1\ell_1=n_1,~k_2\ell_2=n_2}f(k_1)f(k_2)g(\ell_1)g(\ell_2)\\ &= g(n_1)g(n_2) - \left(\sum_{k_1\ell_1=n_1}f(k_1)g(\ell_1)\right)\left(\sum_{k_2\ell_2=n_2}f(k_2)g(\ell_2)\right) \\ &= g(n_1)g(n_2) - \varepsilon(n_1)\varepsilon(n_2)\\ &= g(n_1)g(n_2). \end{aligned} \]

The second equality uses the induction hypothesis: for \(\ell_1\ell_2 < n_1n_2\) and \(\ell_1\perp\ell_2\), the condition \(g(\ell_1\ell_2)=g(\ell_1)g(\ell_2)\) holds.

In the language of abstract algebra, the set of all multiplicative functions forms a subgroup of the multiplicative group of the Dirichlet ring.

More special are completely multiplicative functions.

Theorem

Let \(\alpha\) be a completely multiplicative function, and \(f,g\) be arithmetic functions. Then:

  1. Distributivity: \((\alpha f)\ast(\alpha g) = \alpha\cdot(f\ast g)\).
  2. Inverse: \((\alpha f)^{-1}=\alpha f^{-1}\), provided \(f^{-1}\) exists.
  3. A multiplicative function \(f\) is completely multiplicative if and only if \(f^{-1}=\mu f\), where \(\mu\) is the Möbius function.
Proof

For the first point, direct verification shows

\[ \begin{aligned} ((\alpha f)\ast(\alpha g))(n) &= \sum_{k\ell = n}(\alpha f)(k)(\alpha g)(\ell) \\ &= \sum_{k\ell = n}\alpha(k)f(k)\alpha(\ell)g(\ell) \\ &= \sum_{k\ell = n}\alpha(n)f(k)g(\ell) \\ &= \alpha(n)(f\ast g)(n). \end{aligned} \]

The third equality uses the property of completely multiplicative functions: \(\alpha(n)=\alpha(k)\alpha(\ell)\) holds for all \(n=k\ell\).

For the second point, using the first point we have

\[ (\alpha f)\ast(\alpha f^{-1}) = \alpha(f\ast f^{-1}) = \alpha\varepsilon = \varepsilon. \]

The last equality only uses \(\alpha(1)=1\). By the definition of inverse, \((\alpha f)^{-1}=\alpha f^{-1}\).

For the third point, using the second point and \(1^{-1}=\mu\), if \(f\) is completely multiplicative, then

\[ f^{-1} = (1f)^{-1} = 1^{-1}\cdot f = \mu f. \]

Here, \(1\) is the constant function. Conversely, if \(f\) is multiplicative and \(f^{-1}=\mu f\), then to prove \(f\) is completely multiplicative, it suffices to prove that for all primes \(p\) and \(e\in\mathbf N_+\), \(f(p^e)=f(p)^e\) holds. We prove this by induction on \(e\in\mathbf N_+\). The base case \(e=1\) is trivial. For any \(e > 1\), applying the inverse recurrence formula, we have

\[ \begin{aligned} f^{-1}(p^e) &= -\sum_{i=1}^{e}f(p^i)f^{-1}(p^{e-i}) \\ &= -\sum_{i=1}^{e}f(p^i)\mu(p^{e-i})f(p^{e-i})\\ &= -f(p^e)f(1)\mu(1) - f(p^{e-1})\mu(p)f(p)\\ &= -f(p^e) + f(p)^e. \end{aligned} \]

The last equality uses the induction hypothesis \(f(p^{e-1})=f(p)^{e-1}\). Applying \(f^{-1}=\mu f\), we get

\[ f^{-1}(p^e) = \mu(p^e)f(p^e) = 0. \]

Substituting into the previous equation yields

\[ f(p^e) = f(p)^e. \]

Thus, the induction step holds. The proposition is proved.

In the language of abstract algebra, if \(\alpha\) is completely multiplicative, the mapping \(f\mapsto \alpha f\) is an endomorphism of the Dirichlet ring.

Dirichlet Generating Functions

Closely related to Dirichlet convolution are Dirichlet generating functions.

The Dirichlet generating function (DGF) of an arithmetic function \(f(n)\)—that is, the sequence \(\{f(n)\}\)—is defined as the formal Dirichlet series:

\[ F(s) = \sum_{n=1}^{\infty}\dfrac{f(n)}{n^s}. \]

Here \(s\) is a formal variable. In common Dirichlet generating functions, \(s\) can often be regarded as a complex variable, allowing discussion of the analytic properties of Dirichlet series, but this is beyond the scope of competitive programming.

The product of Dirichlet generating functions corresponds to the Dirichlet convolution of the respective arithmetic functions:

Theorem

For arithmetic functions \(f,g\) and their Dirichlet generating functions \(F,G\), the generating function of their Dirichlet convolution \(f\ast g\) equals \(F\cdot G\).

Proof

Direct verification:

\[ \begin{aligned} F(s)G(s) &= \left(\sum_{k=1}^\infty\dfrac{f(k)}{k^s}\right)\left(\sum_{\ell=1}^\infty\dfrac{g(\ell)}{\ell^s}\right)= \sum_{k=1}^\infty\sum_{\ell=1}^\infty\dfrac{f(k)g(\ell)}{(k\ell)^s}\\ &= \sum_{n=1}^{\infty}\dfrac{\sum_{k\ell = n}f(k)g(\ell)}{n^s} = \sum_{n=1}^\infty\dfrac{(f\ast g)(n)}{n^s}. \end{aligned} \]

Utilizing the correspondence between Dirichlet convolution and the product of Dirichlet generating functions, we can understand the properties of Dirichlet convolution from the perspective of Dirichlet generating functions. Since the multiplication of formal Dirichlet series satisfies commutativity, associativity, and distributivity over addition, the Dirichlet convolution operation on arithmetic functions satisfies the same algebraic properties.

Euler Product

The special nature of multiplicative functions is also reflected in their Dirichlet generating functions. Due to the unique factorization theorem, the generating function \(F(s)\) of a multiplicative function \(f(n)\) can be written as:

\[ \begin{aligned} F(s) &= \sum_{n=1}^{\infty}\dfrac{f(n)}{n^s} = \sum_{n=1}^{\infty}\prod_{p\in\mathbf P}\dfrac{f(p^{e})}{p^{es}} = \prod_{p\in\mathbf P}\sum_{e=0}^{\infty}\dfrac{f(p^e)}{p^{es}}\\ &= \prod_{p\in\mathbf P}\left(1 + \dfrac{f(p)}{p^s} + \dfrac{f(p^2)}{p^{2s}} + \dfrac{f(p^3)}{p^{3s}} + \cdots\right). \end{aligned} \]

This means \(F(s)\) can be decomposed into a product of factors \(F_p(s)\), where each \(F_p(s)\) corresponds to an arithmetic function that may only be nonzero at powers of \(p\). This infinite product is also called the Euler product. If \(F(s)\) and \(G(s)\) can both be decomposed similarly, then their product can as well; translating this observation to arithmetic functions yields that the Dirichlet convolution of multiplicative functions remains multiplicative.

Furthermore, if \(f(n)\) is also completely multiplicative, then \(f(p^e)=f(p)^e\), and the above expression can be further simplified:

\[ F(s) = \prod_{p\in\mathbf P}\sum_{e=0}^{\infty}\dfrac{f(p)^e}{p^{es}} = \prod_{p\in\mathbf P}\left(1-\dfrac{f(p)}{p^s}\right)^{-1}. \]

Unlike multiplicative functions, the form of the Dirichlet generating function of a completely multiplicative function is not closed under multiplication. Therefore, the Dirichlet convolution and Dirichlet inverse of completely multiplicative functions are not necessarily completely multiplicative, but they are always multiplicative.

Examples
  1. The unit function \(\varepsilon(n)\) is completely multiplicative. Its Dirichlet generating function is the constant function with respect to the indeterminate \(s\):

    \[ E(s) = \sum_{n=1}^{\infty}\dfrac{\varepsilon(n)}{n^s} = 1. \]
  2. The constant function \(1(n)\) is completely multiplicative. Its Dirichlet generating function is the Riemann zeta function:

    \[ I(s) = \sum_{n=1}^{\infty}\dfrac{1}{n^s} = \prod_{p\in\mathbf P}\dfrac{1}{1-p^{-s}} = \zeta(s). \]
  3. The Möbius function \(\mu(n)\) is the Dirichlet inverse of the constant function. Its Dirichlet generating function is the reciprocal of \(\zeta(s)\):

    \[ M(s) = \sum_{n=1}^{\infty}\dfrac{\mu(n)}{n^s} = \prod_{p\in\mathbf P}(1-p^{-s}) = \dfrac{1}{\zeta(s)}. \]
  4. The power function \(\operatorname{id}_k(n)=n^k\) is completely multiplicative. In particular, when \(k=0\), it is the constant function; when \(k=1\), it is the identity function. Its Dirichlet generating function is:

    \[ I_k(s) = \sum_{n=1}^{\infty}\dfrac{n^k}{n^s} = \prod_{p\in\mathbf P}\dfrac{1}{1-p^{k-s}} = \zeta(s-k). \]
  5. Euler's totient function \(\varphi(n)\) is multiplicative. Its Dirichlet generating function is:

    \[ \begin{aligned} \Phi(s) &= \prod_{p\in\mathbf P}\left(1+\dfrac{p-1}{p^s} + \dfrac{p(p-1)}{p^{2s}} + \dfrac{p^2(p-1)}{p^{3s}}+\cdots\right)\\ &= \prod_{p\in\mathbf P}\left(\dfrac{1}{1-p^{1-s}}-\dfrac{1}{p^s}\dfrac{1}{1-p^{1-s}}\right) = \prod_{p\in\mathbf P}\dfrac{1-p^{-s}}{1-p^{1-s}} = \dfrac{\zeta(s-1)}{\zeta(s)}. \end{aligned} \]

    Combining with the Dirichlet function expression for the power function yields \(\mathrm{id} = \varphi\ast 1\).

  6. The divisor function \(\sigma_k(n)=\sum_{d\mid n}d^k\) is multiplicative. Its Dirichlet generating function is:

    \[ \begin{aligned} \Sigma_k(s) &= \prod_{p\in\mathbf P}\left(1+\dfrac{1+p^k}{p^s}+\dfrac{1+p^k+p^{2k}}{p^{2s}}+\dfrac{1+p^k+p^{2k}+p^{3k}}{p^{3s}}+\cdots\right) \\ &= \prod_{p\in\mathbf P}\dfrac{1}{1-p^k}\left((1-p^k)+\dfrac{1-p^{2k}}{p^s}+\dfrac{1-p^{3k}}{p^{2s}}+\dfrac{1-p^{4k}}{p^{3k}}+\cdots\right)\\ &= \prod_{p\in\mathbf P}\dfrac{1}{1-p^k}\left(\dfrac{1}{1-p^{-s}} - \dfrac{p^k}{1-p^{k-s}}\right)\\ &= \prod_{p\in\mathbf P}\dfrac{1}{(1-p^{-s})(1-p^{k-s})} = \zeta(s-k)\zeta(s). \end{aligned} \]

    Combining with the Dirichlet expression for the power function yields \(\sigma_k = \mathrm{id}_k\ast 1\). This is exactly the definition of \(\sigma_k\).

  7. The indicator function for square-free numbers \(u(n)=|\mu(n)|\) is multiplicative. Its Dirichlet generating function is:

    \[ U(s) = \prod_{p\in\mathbf P}(1+p^{-s}) = \prod_{p\in\mathbf P}\dfrac{1-p^{-2s}}{1-p^{-s}} = \dfrac{\zeta(s)}{\zeta(2s)}. \]

Applications

Dirichlet generating functions can be used to express multiplicative functions as Dirichlet convolutions.

For example, in the process of the Du's sieve, to compute the prefix sum of a multiplicative function \(f\), we need to find another multiplicative function \(g\) such that both \(f\ast g\) and \(g\) can have their prefix sums computed quickly. This process can be derived using Dirichlet generating functions.

Taking the example from the Du's sieve section Luogu P3768 Simple Math Problem, we need to construct an arithmetic function \(g(n)\) for \(f(n)=n^2\varphi(n)\) satisfying the above conditions. Since \(f\) is multiplicative, its Dirichlet generating function is

\[ F(s) = \prod_{p\in\mathbf P}\left(1 + \sum_{k=1}^{\infty}\dfrac{p^{3k-1}(p-1)}{p^{ks}}\right) = \prod_{p\in\mathbf P}\dfrac{1-p^{2-s}}{1-p^{3-s}} = \dfrac{\zeta(s-3)}{\zeta(s-2)}. \]

Comparing with the Dirichlet generating function of the power function, we see that taking \(g = \mathrm{id}_2\) yields \(f \ast g = \mathrm{id}_3\). Both can have their prefix sums computed quickly.

Computation of Dirichlet Convolution

This section discusses the computation of Dirichlet convolution, i.e., given sequences \(\{f(k)\}_{k=1}^n\) and \(\{g(k)\}_{k=1}^n\), solving for the first several terms \(\{h(k)\}_{k=1}^n\) of the Dirichlet convolution \(h=f\ast g\). Depending on the properties of the functions involved, the algorithm's complexity varies.

General Case

If \(f,g,h\) have no special properties, then Dirichlet convolution can only be computed using its definition:

\[ h(n) = \sum_{k\ell = n}f(k)g(\ell). \]

Enumerate \(k\) and \(\ell\), and accumulate the contribution \(f(k)g(\ell)\) to \(h(k\ell)\). The enumeration complexity is

\[ O\left(\sum_{k=1}^{n}\dfrac{n}{k}\right) = O(n\log n). \]

A reference implementation is as follows:

Reference Implementation
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
// Compute the Dirichlet convolution h = f * g.
auto dirichlet_convolute(const std::vector<int>& f, const std::vector<int>& g) {
  int n = f.size() - 1;
  std::vector<int> h(n + 1);
  for (int k = 1; k <= n; ++k) {
    for (int d = 1; k * d <= n; ++d) {
      h[k * d] += f[k] * g[d];
    }
  }
  return h;
}

Convolution with a Multiplicative Function

If \(g\) is multiplicative, then the Euler product can be used to accelerate the computation of Dirichlet convolution. Computing \(h\) is equivalent to computing the coefficients in its Dirichlet generating function \(H\). Since

\[ H(s) = F(s)G(s) = F(s)\prod_{p\in\mathbf P}G_p(s). \]

Here, \(G_p(s)\) is the factor in the Euler product decomposition of \(G(s)\), containing only coefficients at powers of \(p\):

\[ G_p(s) = \sum_{p^k\le n}\dfrac{f(p^k)}{p^{ks}} = 1 + \dfrac{f(p)}{p^s} + \dfrac{f(p^2)}{p^{2s}} + \cdots. \]

Starting from \(F(s)\), iterate over all primes \(p\) not exceeding \(n\), and multiply \(G_p(s)\) one by one to obtain the final result \(H(s)\). When multiplying \(G_p(s)\), directly apply the brute-force enumeration algorithm from the general case. The total enumeration count is

\[ \sum_{p\in\mathbf P,~p\le n}\sum_{k=1}^{\infty}\left\lfloor\dfrac{n}{p^k}\right\rfloor \le \sum_{p\in\mathbf P,~p\le n}\dfrac{n}{p-1} \le \sum_{p\in\mathbf P,~p\le n}\dfrac{2n}{p} \in O(n\log\log n). \]

The complexity estimate in the last step is similar to the proof for the complexity of the Sieve of Eratosthenes. Therefore, the time complexity of this algorithm is \(O(n\log\log n)\).

A reference implementation is as follows:

Reference Implementation
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
// Compute the Dirichlet convolution h = f * g.
// Assume that g is multiplicative.
auto dirichlet_convolute(const std::vector<int>& f, const std::vector<int>& g) {
  int n = f.size() - 1;
  std::vector<int> h(f);
  std::vector<bool> vis(n + 1);
  for (int i = 2; i <= n; ++i) {
    if (vis[i]) continue;
    // Reverse the order for in-place computation.
    for (int k = n / i * i; k; k -= i) {
      vis[k] = true;
      int d = k;
      while (true) {
        d /= i;
        h[k] += h[d] * g[k / d];
        if (d % i) break;
      }
    }
  }
  return h;
}

In particular, when the multiplicative function \(g\) is completely multiplicative or its Dirichlet inverse, e.g., when \(g = 1\) or \(g = \mu\), the algorithm can be further simplified. In this case, the computation of Dirichlet convolution \(h = f\ast g\) can use the more constant-efficient Dirichlet prefix sum/difference algorithm, but the time complexity remains \(O(n\log\log n)\).

Result is a Multiplicative Function

Finally, consider the case where \(h\) is multiplicative. In particular, when \(f,g\) are both multiplicative, \(h=f \ast g\) is multiplicative. To compute \(h\), it suffices to determine its values at prime powers, and then compute it in \(O(n)\) time via linear sieve. The values \(h(p^e)\) at prime powers \(p^e\) can be computed directly by brute force:

\[ h(p^e) = \sum_{i=0}^e f(p^i)g(p^{e-i}). \]

The enumeration required for these brute-force computations is

\[ \begin{aligned} \sum_{p\in\mathbf P,~p\le n}\sum_{e=1}^{\lfloor\log_p n\rfloor}(e+1) &\le \sum_{p\in\mathbf P,~p\le\sqrt{n}}\lfloor\log_p n\rfloor^2 + \sum_{p\in\mathbf P,~\sqrt{n} < p\le n}1 \\ &\le \sqrt{n}(\log_2 n)^2 + n \in O(n). \end{aligned} \]

Thus, the total time complexity of this algorithm is \(O(n)\).

A reference implementation is as follows:

Reference Implementation
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
// Compute the Dirichlet convolution h = f * g.
// Assume that h is multiplicative.
auto dirichlet_convolute(const std::vector<int>& f, const std::vector<int>& g) {
  int n = f.size() - 1;
  std::vector<int> h(n + 1), primes, rem(n + 1), lpf(n + 1);
  std::vector<bool> vis(n + 1);
  h[1] = 1;
  for (int x = 2; x <= n; ++x) {
    if (!vis[x]) {
      primes.push_back(x);
      rem[x] = 1;
      lpf[x] = x;
    }
    for (int p : primes) {
      if (x * p > n) break;
      vis[x * p] = true;
      rem[x * p] = x % p ? x : rem[x];
      lpf[x * p] = p;
      if (x % p == 0) break;
    }
    if (rem[x] == 1) {  // prime powers.
      for (int k = x; k; k /= lpf[x]) {
        h[x] += f[k] * g[x / k];
      }
    } else {  // other cases.
      h[x] = h[rem[x]] * h[x / rem[x]];
    }
  }
  return h;
}

References and Notes