Möbius Inversion

Prerequisites: Number Theoretic Block, Dirichlet Convolution

Möbius inversion is an important topic in number theory. For some function \(f(n)\), if it is difficult to directly compute its value, but easy to compute its sum over multiples or divisors \(g(n)\), then Möbius inversion can simplify the computation and find the value of \(f(n)\).

Möbius Function

The Möbius function is defined as

\[ \mu(n)= \begin{cases} 1,&n=1,\\ 0,&n\text{ is divisible by a square }>1,\\ (-1)^k,&n\text{ is the product of }k\text{ distinct primes}.\\ \end{cases} \]

Specifically, suppose the positive integer \(n\) has prime factorization \(n=\prod_{i=1}^kp_i^{e_i}\), where \(p_i\) are primes and \(e_i\) are positive integers. Then the three cases correspond to:

\(\mu(1) = 1\);
When there exists \(i\) such that \(e_i > 1\), i.e., any prime factor appears more than once, \(\mu(n)=0\);
Otherwise, for all \(i\), \(e_i = 1\), i.e., every prime factor appears exactly once, \(\mu(n)=(-1)^k\), where \(k\) is the number of distinct prime factors.

Properties

From the definition, it is easy to verify that the Möbius function \(\mu(n)\) is multiplicative but not completely multiplicative. Besides, the most important property is the following identity:

Property

For a positive integer \(n\), we have

\[ \sum_{d\mid n}\mu(d) = [n = 1] = \begin{cases} 1,&n=1,\\ 0,&n\neq 1.\\ \end{cases} \]

where \([\cdot]\) is the Iverson bracket.

Proof

Let \(n=\prod_{i=1}^kp_i^{e_i}\), and set \(n' = \prod_{i=1}^kp_i\). By the binomial theorem, we have

\[ \sum_{d\mid n}\mu(d) = \sum_{d\mid n'}\mu(d) = \sum_{i=0}^k\binom{k}{i}(-1)^i = (1 + (-1))^k = [k = 0] = [n = 1]. \]

Using Dirichlet convolution, this expression can be written as \(\varepsilon = 1 * \mu\). That is, the Möbius function is the Dirichlet inverse of the constant function \(1\).

This property has a common application:

\[ [i\perp j] = [\gcd(i,j) = 1] = \sum_{d\mid\gcd(i,j)} \mu(d) = \sum_{d}[d\mid i][d\mid j]\mu(d). \]

It transforms the coprime condition into a sum involving the Möbius function, facilitating further derivation.

Computation

If we need to compute the Möbius function \(\mu(n)\) for a single \(n\), we can use its prime factorization. For example, when \(n\) is not too large, we can compute \(\mu(n)\) in \(O(\sqrt{n})\) time.

Reference Implementation

C++Python

int mu(int n) {
  int res = 1;
  for (int i = 2; i * i <= n; ++i) {
    if (n % i == 0) {
      n /= i;
      // Check if square-free.
      if (n % i == 0) return 0;
      res = -res;
    }
  }
  // The remaining factor must be prime.
  if (n > 1) res = -res;
  return res;
}

def mu(n):
    res = 1
    i = 2
    while i * i <= n:
        if n % i == 0:
            n //= i
            # Check if square-free
            if n % i == 0:
                return 0
            res = -res
        i += 1
    # The remaining factor must be prime
    if n > 1:
        res = -res
    return res

If we need to precompute \(\mu(n)\) for the first \(n\) positive integers, we can use the fact that it is multiplicative and compute it via linear sieve in \(O(n)\) time.

Reference Implementation

C++Python

std::vector<int> get_mu(int n) {
  std::vector<int> mu(n + 1), primes;
  std::vector<bool> not_prime(n + 1);
  primes.reserve(n);
  mu[1] = 1;
  for (int x = 2; x <= n; ++x) {
    if (!not_prime[x]) {
      primes.push_back(x);
      mu[x] = -1;
    }
    for (int p : primes) {
      if (x * p > n) break;
      not_prime[x * p] = true;
      if (x % p == 0) {
        mu[x * p] = 0;
        break;
      } else {
        mu[x * p] = -mu[x];
      }
    }
  }
  return mu;
}

def get_mu(n):
    mu = [0] * (n + 1)
    primes = []
    not_prime = [False] * (n + 1)

    mu[1] = 1
    for x in range(2, n + 1):
        if not not_prime[x]:
            primes.append(x)
            mu[x] = -1
        for p in primes:
            if x * p > n:
                break
            not_prime[x * p] = True
            if x % p == 0:
                mu[x * p] = 0
                break
            else:
                mu[x * p] = -mu[x]
    return mu

Möbius Inversion

The most important application of the Möbius function is Möbius inversion.

Möbius Inversion

Let \(f(n),g(n)\) be two arithmetic functions. Then,

\[ f(n) = \sum_{d\mid n}g(d) \iff g(n) = \sum_{d\mid n}\mu\left(\dfrac{n}{d}\right)f(d). \]

Proof 1

Direct verification:

\[ \begin{aligned} \sum_{d\mid n}\mu\left(\dfrac{n}{d}\right)f(d) &= \sum_{d\mid n}\mu\left(\dfrac{n}{d}\right)\sum_{k\mid d}g(k)\\ &= \sum_{k\mid n}g(k)\sum_d[k\mid d\mid n]\mu\left(\dfrac{n}{d}\right)\\ &= \sum_{k\mid n}g(k)\sum_{d\mid n}\left[\frac{n}{d}\mid\frac{n}{k}\right]\mu\left(\dfrac{n}{d}\right)\\ &= \sum_{k\mid n}g(k)\left[\frac{n}{k} = 1\right] \\ &= g(n). \end{aligned} \]

The key step is swapping the order of summation and noting that \(k\mid d\mid n\) is equivalent to \(\dfrac{n}{d}\mid\dfrac{n}{k}\). The second last equality sums the Möbius function over divisors \(\dfrac{n}{d}\) of \(\dfrac{n}{k}\), which equals \(\left[\dfrac{n}{k} = 1\right]\). This expression is nonzero only when \(n=k\), yielding \(g(n)\).

Proof 2

Using Dirichlet convolution, the proposition is equivalent to

\[ f = 1 * g \iff g = \mu * f. \]

Using \(1 * \mu = \varepsilon\), convolving both sides of the left equation with \(\mu\) gives

\[ f * \mu = (1 * g) * \mu = (1 * \mu) * g = \varepsilon * g = g. \]

In summations involving divisibility relations of arithmetic functions, Möbius inversion is a powerful transformation tool.

Examples

The Euler totient function \(\varphi(n)\) satisfies \(n = \sum_{d\mid n}\varphi(d)\), i.e., \(\mathrm{id}=1*\varphi\). Inverting it gives \(\varphi = \mu * \mathrm{id}\), i.e.,

\[ \varphi(n) = \sum_{d\mid n}d\mu\left(\dfrac{n}{d}\right). \]
The divisor function \(\sigma_k(n) = \sum_{d\mid n}d^k\), i.e., \(\sigma_k = 1 * \mathrm{id}_k\). Inverting it gives \(\mathrm{id}_k = \mu * \sigma_k\), i.e.,

\[ n^k = \sum_{d\mid n}\mu\left(\dfrac{n}{d}\right)\sigma_k(d). \]
The distinct prime factor count function \(\omega(n)=\sum_{d\mid n}[d\in\mathbf P]\), i.e., \(\omega = 1* \mathbf{1}_{\mathbf P}\), where \(\mathbf{1}_{\mathbf P}\) is the indicator function of the prime set \(\mathbf P\). Inverting it gives \(\mathbf{1}_{\mathbf P} = \mu * \omega\), i.e.,

\[ [n\in\mathbf P] = \sum_{d\mid n}\mu\left(\dfrac{n}{d}\right)\omega(d). \]
Consider the arithmetic function \(\Lambda(n)\) satisfying \(\log n = \sum_{d\mid n}\Lambda(d)\). It is the Möbius inversion of the logarithm function, also known as the von Mangoldt function:

\[ \Lambda(n) = \sum_{d\mid n}\mu\left(\dfrac{n}{d}\right)\log d = \begin{cases} \log p, & n = p^e,~p\in\mathbf P,~e\in\mathbf N_+, \\ 0, &\text{otherwise}. \end{cases} \]

Appendix: Proof of \(\Lambda(n)\) expression

For prime powers \(n=p^e~(e\in\mathbf N_+)\), we have

\[ \Lambda(n) = \sum_{i=0}^e\mu(p^{e-i})\log p^i = \log p^{e} - \log p^{e-1} = \log p. \]

For \(n=1\), obviously \(\Lambda(n)=\log 1=0\). For other composite \(n\), we have

\[ \Lambda(n) = \sum_{d\mid n}\mu(d)(\log n-\log d) = \left(\sum_{d\mid n}\mu(d)\right)\log n-\sum_{d\mid n}\mu(d)\log d. \]

By the property of the Möbius function, the coefficient of \(\log n\) is \([n=1]=0\). For the latter sum, further factor \(d\) into prime factors. For any prime \(p\mid n\), consider the coefficient of \(\log p\):

\[ -\sum_{p\mid d\mid n}\mu(d) = \sum_{(d/p)\mid(n/p)}\mu\left(\dfrac{d}{p}\right) = \left[\dfrac{n}{p}=1\right]=0. \]

Thus, for composite \(n\) with more than one prime factor, \(\Lambda(n)=0\).

Extended Forms

Besides the basic form, Möbius inversion has some common extended forms. First, consider its sum-over-multiples form.

Extension 1

Let \(f(n),g(n)\) be two arithmetic functions. Then,

\[ f(n) = \sum_{n\mid d}g(d) \iff g(n) = \sum_{n\mid d}\mu\left(\dfrac{d}{n}\right)f(d). \]

Proof

Direct verification:

\[ \begin{aligned} \sum_{n\mid d}\mu\left(\dfrac{d}{n}\right)f(d) &= \sum_{n\mid d}\mu\left(\dfrac{d}{n}\right)\sum_{d\mid k}g(k)\\ &= \sum_{n\mid k}g(k)\sum_{d}[n\mid d\mid k]\mu\left(\dfrac{d}{n}\right)\\ &= \sum_{n\mid k}g(k)\sum_{n\mid d}\left[\dfrac{d}{n}\mid\dfrac{k}{n}\right]\mu\left(\dfrac{d}{n}\right)\\ &= \sum_{n\mid k}g(k)\left[\dfrac{k}{n}=1\right]\\ &= g(n). \end{aligned} \]

This derivation is dual to the basic form.

Second, Möbius inversion is not limited to addition; it holds for any operation in an Abelian group. For example, it has the following multiplicative form:

Extension 2

Let \(f(n),g(n)\) be two arithmetic functions. Then,

\[ f(n) = \prod_{d\mid n}g(d) \iff g(n) = \prod_{d\mid n}f(d)^{\mu(n/d)}. \]

Proof

Direct verification:

\[ \begin{aligned} \prod_{d\mid n}f(d)^{\mu(n/d)} &= \prod_{d\mid n}\left(\prod_{k\mid d}g(k)\right)^{\mu(n/d)}\\ &= \prod_{k\mid n}g(k)\uparrow\left(\sum_d[k\mid d\mid n]\mu\left(\dfrac{n}{d}\right)\right)\\ &= \prod_{k\mid n}g(k)\uparrow\left(\sum_{d\mid n}\left[\frac{n}{d}\mid\frac{n}{k}\right]\mu\left(\dfrac{n}{d}\right)\right)\\ &= \prod_{k\mid n}g(k)\uparrow\left[\frac{n}{k} = 1\right] \\ &= g(n). \end{aligned} \]

Here, \(a\uparrow b = a^b\) is the Knuth arrow. Compared to the proof of the basic form, the only difference is addition replaced by multiplication and multiplication replaced by exponentiation.

From the perspective of Dirichlet convolution, Möbius inversion merely uses the fact that "the Möbius function is the Dirichlet inverse of the constant function". It is easy to imagine that similar inversion relations hold for general Dirichlet inverses.

Extension 3

Let \(f(n),g(n),\alpha(n)\) be arithmetic functions, and let \(\alpha^{-1}(n)\) be the Dirichlet inverse of \(\alpha(n)\), i.e.,

\[ [n=1] = \sum_{d\mid n}\alpha\left(\dfrac{n}{d}\right)\alpha^{-1}(d). \]

Then,

\[ f(n) = \sum_{d\mid n}\alpha\left(\dfrac{n}{d}\right)g(d) \iff g(n) = \sum_{d\mid n}\alpha^{-1}\left(\dfrac{n}{d}\right)f(d). \]

Proof

Direct verification:

\[ \begin{aligned} \sum_{d\mid n}\alpha^{-1}\left(\dfrac{n}{d}\right)f(d) &= \sum_{d\mid n}\alpha^{-1}\left(\dfrac{n}{d}\right)\sum_{k\mid d}\alpha\left(\dfrac{d}{k}\right)g(k)\\ &= \sum_{k\mid n}g(k)\sum_d[k\mid d\mid n]\alpha\left(\dfrac{d}{k}\right)\alpha^{-1}\left(\dfrac{n}{d}\right)\\ &= \sum_{k\mid n}g(k)\sum_{d\mid n}\left[\frac{n}{d}\mid\frac{n}{k}\right]\alpha\left(\dfrac{d}{k}\right)\alpha^{-1}\left(\dfrac{n/k}{d/k}\right)\\ &= \sum_{k\mid n}g(k)\left[\frac{n}{k} = 1\right] \\ &= g(n). \end{aligned} \]

Compared to the proof of the basic form, only the second last equality is replaced by the definition of Dirichlet inverse.

Corollary

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

\[ f(n) = \sum_{d\mid n}t\left(\dfrac{n}{d}\right)g(d) \iff g(n) = \sum_{d\mid n}\mu\left(\dfrac{n}{d}\right)t\left(\dfrac{n}{d}\right)f(d). \]

Proof

By the properties of Dirichlet convolution, for a completely multiplicative function \(t(n)\), its Dirichlet inverse is \(\mu(n)t(n)\).

Finally, Möbius inversion can be generalized to complex-valued functions on \([1,+\infty)\), not just arithmetic functions. The basic form of Möbius inversion can be seen as a special case where the complex-valued function is zero at all non-integer points.

Extension 4

Let \(F(x)\) and \(G(x)\) be complex-valued functions on \([1,+\infty)\). Then,

\[ F(x) = \sum_{n = 1}^{\lfloor x\rfloor}G\left(\dfrac{x}{n}\right) \iff G(x) = \sum_{n = 1}^{\lfloor x\rfloor}\mu(n)F\left(\dfrac{x}{n}\right). \]

Proof

Extend the definitions of \(F\) and \(G\) by setting \(F(x)=G(x)=0\) for \(x < 1\). Then the proposition is equivalent to:

\[ F(x) = \sum_n G\left(\dfrac{x}{n}\right) \iff G(x) = \sum_n \mu(n)F\left(\dfrac{x}{n}\right). \]

The sums are over \(n\in\mathbf N_+\).

Direct verification:

\[ \begin{aligned} \sum_n \mu(n)F\left(\dfrac{x}{n}\right) &= \sum_n\mu(n)\sum_d G\left(\dfrac{x/n}{d}\right)\\ &= \sum_k G\left(\dfrac{x}{k}\right)\sum_{n\mid k}\mu(n)\\ &= \sum_k G\left(\dfrac{x}{k}\right)[k=1]\\ &= G(x). \end{aligned} \]

To obtain the second equality, set \(k = nd\).

Corollary

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

\[ f(n) = \sum_{k=1}^ng\left(\left\lfloor\dfrac{n}{k}\right\rfloor\right) \iff g(n)=\sum_{k=1}^n\mu(k)f\left(\left\lfloor\dfrac{n}{k}\right\rfloor\right). \]

Proof

Simply take \(F(x)=f(\lfloor x\rfloor)\) and \(G(x)=g(\lfloor x\rfloor)\).

These extended forms can be combined to yield more complex inversion relations.

Dirichlet Prefix Sum

Prerequisites: Prefix Sum and Difference

Consider the basic form of Möbius inversion:

\[ f(n) = \sum_{d\mid n}g(d) \iff g(n) = \sum_{d\mid n}\mu\left(\dfrac{n}{d}\right)f(d). \]

In the left equation, \(f(n)\) is the sum of \(g(n)\) over all divisors of \(n\). If we interpret \(a\mid b\) as \(a\) preceding \(b\), then \(f(n)\) can be seen as a kind of prefix sum of \(g(n)\). Therefore, in the Chinese competitive programming community, the process of computing \(\{f(k)\}_{k=1}^n\) from \(\{g(k)\}_{k=1}^n\) is also called Dirichlet prefix sum, and the inverse process is called Dirichlet difference. These methods often appear when preprocessing the values of an arithmetic function at the first \(N\) points.

Next, we discuss the computation of Dirichlet prefix sums. If we treat each prime as a dimension, this becomes a kind of high-dimensional prefix sum. Recall the dimension-by-dimension prefix sum algorithm: iterate over all dimensions one by one, and accumulate the value at each position to its successor in that dimension. For arithmetic functions, this means iterating over all primes \(p\) from small to large, and accumulating the function value at \(n\) to \(np\). This traversal order is the same as in the Eratosthenes sieve. Therefore, this algorithm can compute the Dirichlet prefix sum of a sequence of length \(n\) in \(O(n\log\log n)\) time. Similarly, using dimension-by-dimension difference, the Dirichlet difference can be computed in the same time complexity.

Reference Implementation

Dirichlet Prefix SumDirichlet Difference

std::vector<int> dirichlet_presum(const std::vector<int>& g) {
  int n = g.size() - 1;
  std::vector<int> f(g);
  std::vector<bool> vis(n + 1);
  for (int x = 2; x <= n; ++x) {
    if (vis[x]) continue;
    for (int y = 1, xy = x; xy <= n; ++y, xy += x) {
      vis[xy] = true;
      f[xy] += f[y];
    }
  }
  return f;
}

std::vector<int> dirichlet_diff(const std::vector<int>& f) {
  int n = f.size() - 1;
  std::vector<int> g(f);
  std::vector<int> vis(n + 1);
  for (int x = 2; x <= n; ++x) {
    if (vis[x]) continue;
    for (int y = n / x, xy = x * y; y; --y, xy -= x) {
      vis[xy] = true;
      g[xy] -= g[y];
    }
  }
  return g;
}

This computation method can be generalized to the sum-over-multiples form (Extension 1), multiplicative form (Extension 2), and using completely multiplicative functions instead of the constant function (corollary of Extension 3).

Examples

This section demonstrates the application of Möbius inversion and some common transformation techniques through examples. First, an example to familiarize with basic techniques for handling greatest common divisor conditions in summations.

Luogu P2522 [HAOI 2011] Problem b

\(T\) groups of data. For each group, compute:

\[ \sum_{i=x}^{n}\sum_{j=y}^{m}[\gcd(i,j)=k]. \]

Data range: \(1\le T,x,y,n,m,k\le 5\times 10^4\).

Solution

By inclusion-exclusion, the original expression can be split into \(4\) blocks, each of the form

\[ f(n,m,k)=\sum_{i=1}^{n}\sum_{j=1}^{m}[\gcd(i,j)=k]. \]

For such expressions, the following is a standard derivation process: extract common factors, apply Möbius function property, swap summation order.

First, since \(i,j\) can only be multiples of \(k\), we can extract this factor—this is equivalent to substituting \(i=ki'\) and \(j=kj'\), yielding:

\[ f(n,m,k)=\sum_{i=1}^{\lfloor n/k\rfloor}\sum_{j=1}^{\lfloor m/k\rfloor}[\gcd(i,j)=1]. \]

Then using the property of the Möbius function:

\[ [\gcd(i,j)=1] = \sum_{d\mid\gcd(i,j)}\mu(d) = \sum_d[d\mid i][d\mid j]\mu(d). \]

Substituting into the expression and swapping summation order gives:

\[ f(n,m,k)=\sum_d\mu(d)\left(\sum_{i=1}^{\lfloor n/k\rfloor}[d\mid i]\right)\left(\sum_{j=1}^{\lfloor m/k\rfloor}[d\mid j]\right). \]

The advantage of this manipulation is that for fixed \(d\), the terms involving \(i\) and \(j\) are separated and can be summed independently. Next, since

\[ \sum_{i=1}^{\lfloor n/k\rfloor}[d\mid i] = \left\lfloor\dfrac{\lfloor n/k\rfloor}{d}\right\rfloor,~\sum_{j=1}^{\lfloor m/k\rfloor}[d\mid j]=\left\lfloor\dfrac{\lfloor m/k\rfloor}{d}\right\rfloor, \]

we have

\[ f(n,m,k)=\sum_d\mu(d)\left\lfloor\dfrac{\lfloor n/k\rfloor}{d}\right\rfloor\left\lfloor\dfrac{\lfloor m/k\rfloor}{d}\right\rfloor. \]

After preprocessing \(\mu(d)\) with linear sieve and its prefix sum, we can solve via number theoretic blocking. Total time complexity is \(O(N + T\sqrt{N})\), where \(N\) is the upper bound of \(n,m\), and \(T\) is the number of data groups.

Reference Code

#include <algorithm>
#include <iostream>
using namespace std;
constexpr int N = 50000;
int mu[N + 5], p[N + 5];
bool flg[N + 5];

void init() {
  int tot = 0;
  mu[1] = 1;
  for (int i = 2; i <= N; ++i) {
    if (!flg[i]) {
      p[++tot] = i;
      mu[i] = -1;
    }
    for (int j = 1; j <= tot && i * p[j] <= N; ++j) {
      flg[i * p[j]] = true;
      if (i % p[j] == 0) {
        mu[i * p[j]] = 0;
        break;
      }
      mu[i * p[j]] = -mu[i];
    }
  }
  for (int i = 1; i <= N; ++i) mu[i] += mu[i - 1];
}

int solve(int n, int m) {
  int res = 0;
  for (int i = 1, j; i <= min(n, m); i = j + 1) {
    j = min(n / (n / i), m / (m / i));
    res += (mu[j] - mu[i - 1]) * (n / i) * (m / i);  // Formula derived from the derivation
  }
  return res;
}

int main() {
  cin.tie(nullptr)->sync_with_stdio(false);
  int T, a, b, c, d, k;
  init();  // Preprocess the mu array
  cin >> T;
  for (int i = 1; i <= T; i++) {
    cin >> a >> b >> c >> d >> k;
    // According to the inclusion-exclusion principle, the answer in the range 1<=x<=b && 1<=y<=d
    //   minus the answer in the range 1<=x<=b && 1<=y<=c-1 and
    //   minus the answer in the range 1<=x<=a-1 && 1<=y<=d
    //   plus the answer in the range 1<=x<=a-1 && 1<=y<=c-1
    //   gives the answer in the range a<=x<=b && c<=y<=d
    // If this step is unclear, you can draw a coordinate diagram to understand
    cout << solve(b / k, d / k) - solve(b / k, (c - 1) / k) -
                solve((a - 1) / k, d / k) + solve((a - 1) / k, (c - 1) / k)
         << '\n';
  }
  return 0;
}

The next two examples demonstrate handling by enumerating common divisors and using sieve methods to compute values of general multiplicative functions.

SPOJ LCMSUM

\(T\) groups of data. For each group, compute:

\[ \sum_{i=1}^n \operatorname{lcm}(i,n). \]

Data range: \(1\le T\le 3\times 10^5,~1\le n\le 10^6\).

Solution 1

The problem involves least common multiple, but greatest common divisor is often easier to handle. First, transform:

\[ f(n)=\sum_{i=1}^n \operatorname{lcm}(i,n) = \sum_{i=1}^n \frac{i\cdot n}{\gcd(i,n)}. \]

Factor out \(n\) and enumerate the greatest common divisor \(k\):

\[ f(n)=n\sum_{k\mid n}\sum_{i=1}^n\dfrac{i}{k}[\gcd(i,n)=k]. \]

For the inner sum, this is the most common scenario involving gcd. Repeating the standard process:

\[ \begin{aligned} f(n) &= n\sum_{k\mid n}\sum_{i=1}^{n/k}i\left[\gcd\left(i,\dfrac{n}{k}\right)=1\right]\\ &= n\sum_{k\mid n}\sum_{i=1}^{n/k}i\sum_d\mu(d)[d\mid i]\left[d\mid \dfrac{n}{k}\right]\\ &= n\sum_{k\mid n}\sum_d\mu(d)\left[d\mid \dfrac{n}{k}\right]\left(\sum_{i=1}^{n/k}i[d\mid i]\right). \end{aligned} \]

Again, the sum over \(i\) is separated and can be handled independently. The final sum is essentially an arithmetic series sum: (set \(i=di'\))

\[ \sum_{i=1}^{n/k}i[d\mid i] = d\frac{1}{2}\left(\dfrac{n}{kd}+1\right)\dfrac{n}{kd}=:dG\left(\dfrac{n}{kd}\right). \]

Thus, we obtain:

\[ f(n) = n\sum_{k\mid n}\sum_d\mu(d)\left[d\mid \dfrac{n}{k}\right]dG\left(\dfrac{n}{kd}\right). \]

After enumerating common divisors, such double sums are common. There is a fixed handling method: set the product as a new variable \(\ell=kd\), then swap summation order again. Since \(d\mid(n/k)\) is equivalent to \(d\mid\ell\mid n\), the expression becomes:

\[ f(n) = n\sum_{\ell\mid n}G\left(\dfrac{n}{\ell}\right)\sum_{d\mid\ell}\mu(d)d. \]

Let \(F(\ell)=\sum_{d\mid\ell}\mu(d)d\). Then:

\[ f(n) = n\sum_{\ell\mid n}G\left(\dfrac{n}{\ell}\right)F(\ell). \]

Since \(\mu(d)d\) is multiplicative, its convolution with the constant function \(1\), \(F(n)\), is also multiplicative. Although the sum appears in Dirichlet convolution form, \(G(n)\) is not multiplicative. However, \(G(n)\) is a polynomial, so it is a linear combination of completely multiplicative functions. Thus,

\[ f(n) = \dfrac{1}{2}n\left(\sum_{\ell}\left(\dfrac{n}{\ell}\right)^2F(\ell) + \sum_{\ell}\dfrac{n}{\ell}F(\ell)\right). \]

These two terms (excluding coefficients) are multiplicative functions and can be preprocessed via linear sieve (or by sieving the inner function and then using Dirichlet prefix sum in \(O(N\log\log N)\) time). Specifically, let

\[ H_s(n) = \sum_{\ell}\left(\dfrac{n}{\ell}\right)^sF(\ell),~s=1,2. \]

To derive their expressions, it suffices to determine their values at prime powers. For prime \(p\) and positive exponent \(e\),

\[ \begin{aligned} F(p^e) &= \mu(1) + \mu(p)p + \sum_{j=2}^e\mu(p^j)p^j = 1-p,\\ H_s(p^e) &= (p^e)^{s}F(1) + \sum_{j=1}^e(p^{e-j})^sF(p^j) = p^{es} + (1-p)\dfrac{1-p^{es}}{1-p^s},~s=1,2. \end{aligned} \]

In particular, \(H_1(p^e)\equiv 1\) is constant, and

\[ H_2(p^e) = p^{2e} + (1-p)\dfrac{1-p^{2e}}{1-p^2} = H_2(p^{e-1}) + p^{2e} - p^{2e-1}. \]

This is easily computed via linear sieve. After preprocessing \(H_2(n)\) with linear sieve, each query can be answered in \(O(1)\) time via \(f(n)=(n/2)(H_2(n)+1)\). Total time complexity is \(O(N+T)\), where \(N\) is the upper bound of \(n\), and \(T\) is the number of groups.

The reference implementation further optimizes the linear sieve part using specific properties of the expression, which is not mandatory. Using values at prime powers, preprocessing can still be done in \(O(N)\) time. These derivations are detailed in Solution 2.

Solution 2

For this problem, there is a more flexible approach. From Solution 1,

\[ f(n) = n\sum_{k\mid n}\sum_{i=1}^{n/k}i\left[\gcd\left(i,\dfrac{n}{k}\right)=1\right] = n\sum_{k\mid n}F\left(\dfrac{n}{k}\right). \]

If we do not proceed with Möbius inversion here but observe that the sum is actually the sum of integers not exceeding \(d=n/k\) and coprime to \(d\). For \(d>1\), since integers coprime to \(d\) come in pairs, i.e., \(i\) and \(d-i\) are both coprime to \(d\), we have

\[ F(d)=\sum_{i=1}^{n'}i[i\perp d] = \sum_{i=1}^{d}(d-i)[i\perp d] = \dfrac{1}{2}d\sum_{i=1}^{d}[i\perp d] = \dfrac{1}{2}d\varphi(d). \]

For \(d=1\),

\[ F(d)=1=\dfrac{1}{2}+\dfrac{1}{2}d\varphi(d). \]

Thus, the original expression becomes

\[ f(n) = \dfrac{1}{2}n\left(\sum_{d\mid n}d\varphi(d) + 1\right). \]

Since \(G(n)=\sum_{d\mid n}d\varphi(d)\) is the Dirichlet convolution of the multiplicative function \(n\varphi(n)\) and the constant function \(1\), it is also multiplicative and can be preprocessed via linear sieve. To do so, determine its values at prime powers. For prime \(p\) and positive exponent \(e\),

\[ G(p^e) = 1 + \sum_{i=1}^ep^e(p^e-1) = G(p^{e-1}) + p^{2e} - p^{2e-1}. \]

This matches the result derived in Solution 1. The total time complexity remains \(O(N+T)\).

Finally, using the multiplicative function expression for this problem, the linear sieve computation can be further optimized. For prime \(p\),

\[ G(p) = 1 - p + p^2. \]

The key to linear sieve is computing \(G(pn)\) for general \(n\). This splits into two cases. When \(p\perp n\), since \(G\) is multiplicative,

\[ G(pn) = G(p)G(n). \]

Otherwise, when \(p\mid n\), let \(n=p^em\) with \(p\perp m\). Then

\[ \begin{aligned} G(pn) &= G(p^{e+1})G(m)\\ &= G(p^e)G(m) + (p^{2e+2} - p^{2e+1})G(m)\\ &= G(n) + (p^{2e+2} - p^{2e+1})G(m). \end{aligned} \]

Direct verification shows this also holds for \(p\perp n\). Therefore,

\[ G(n) - G\left(\dfrac{n}{p}\right) = (p^{2e}-p^{2e-1})G(m). \]

Substituting,

\[ G(pn) = G(n) + p^2\left(G(n) - G\left(\dfrac{n}{p}\right)\right). \]

This simplifies the linear sieve computation. Of course, this derivation is not necessary; for multiplicative functions without special properties, directly using \(G(pn)=G(p^{e+1})G(m)\) suffices.

Reference Code

#include <iostream>
constexpr int N = 1000000;
int tot, p[N + 5];
long long g[N + 5];
bool flg[N + 5];  // Flag array

void solve() {
  g[1] = 1;
  for (int i = 2; i <= N; ++i) {
    if (!flg[i]) {
      p[++tot] = i;
      g[i] = (long long)1 * i * (i - 1) + 1;
    }
    for (int j = 1; j <= tot && i * p[j] <= N; ++j) {
      flg[i * p[j]] = true;
      if (i % p[j] == 0) {
        g[i * p[j]] =
            g[i] + (g[i] - g[i / p[j]]) * p[j] * p[j];  // Substitute the derived formula
        break;
      }
      g[i * p[j]] = g[i] * g[p[j]];
    }
  }
}

using std::cin;
using std::cout;

int main() {
  cin.tie(nullptr)->sync_with_stdio(false);
  int T, n;
  solve();  // Preprocess the g array
  cin >> T;
  for (int i = 1; i <= T; ++i) {
    cin >> n;
    cout << (g[n] + 1) * n / 2 << '\n';
  }
  return 0;
}

BZOJ 2154 [National Training Team] Crash's Number Table

Compute:

\[ \sum_{i=1}^n\sum_{j=1}^m\operatorname{lcm}(i,j)\mod{20101009}. \]

Data range: \(1\le n,m\le 10^7\).

Solution

Ignore modulus during derivation. Let

\[ f(n,m) = \sum_{i=1}^n\sum_{j=1}^m\operatorname{lcm}(i,j). \]

Convert lcm to gcd, enumerate common divisor, and apply standard process:

\[ \begin{aligned} f(n,m) &= \sum_k\sum_{i=1}^n\sum_{j=1}^m\dfrac{ij}{k}[\gcd(i,j)=k] \\ &= \sum_k\sum_{i=1}^{\lfloor n/k\rfloor}\sum_{j=1}^{\lfloor m/k\rfloor} kij[\gcd(i,j)=1]\\ &= \sum_k\sum_{i=1}^{\lfloor n/k\rfloor}\sum_{j=1}^{\lfloor m/k\rfloor} kij\sum_d\mu(d)[d\mid i][d\mid j]\\ &= \sum_kk\sum_d\mu(d)\left(\sum_{i=1}^{\lfloor n/k\rfloor}i[d\mid i]\right)\left(\sum_{j=1}^{\lfloor m/k\rfloor}j[d\mid j]\right). \end{aligned} \]

Again, the sums separate for \(i\) and \(j\). Compute these inner sums first by extracting factors (set \(i=di'\)):

\[ \sum_{i=1}^{\lfloor n/k\rfloor}i[d\mid i] = d\sum_{i=1}^{\lfloor\lfloor n/k\rfloor/d\rfloor}i = dG\left(\left\lfloor\dfrac{\lfloor n/k\rfloor}{d}\right\rfloor\right) = dG\left(\left\lfloor\dfrac{n}{kd}\right\rfloor\right). \]

where \(G(n)=\dfrac{1}{2}n(n+1)\) is the arithmetic series sum; the last equality uses properties of the floor function. Symmetrically for the other sum. Substituting back,

\[ f(n,m) = \sum_k k\sum_{d}\mu(d)d^2G\left(\left\lfloor\dfrac{n}{kd}\right\rfloor\right)G\left(\left\lfloor\dfrac{m}{kd}\right\rfloor\right). \]

As before, for such expressions enumerating common divisors, we often enumerate the product \(\ell = kd\) and swap summation order again:

\[ f(n,m) = \sum_{\ell}\left(\sum_{d\mid\ell}\mu(d)d\ell\right)G\left(\left\lfloor\dfrac{n}{\ell}\right\rfloor\right)G\left(\left\lfloor\dfrac{m}{\ell}\right\rfloor\right). \]

Let

\[ F(\ell) = \sum_{d\mid\ell}\mu(d)d\ell. \]

This is the product of the multiplicative function \(\ell\) and the multiplicative function \(\sum_{d\mid\ell}\mu(d)d\), so it is also multiplicative and can be preprocessed via linear sieve, along with its prefix sum. Then \(f(n,m)\) can be computed via number theoretic blocking. Total time complexity is \(O(\min\{n,m\})\).

Reference Code

#include <algorithm>
#include <iostream>
#include <vector>

int main() {
  constexpr int M = 20101009;
  int n, m;
  std::cin >> n >> m;
  if (n > m) std::swap(n, m);
  std::vector<int> f(n + 1), vis(n + 1), prime;
  prime.reserve(n);
  f[1] = 1;
  for (int x = 2; x <= n; ++x) {
    if (!vis[x]) {
      prime.emplace_back(x);
      f[x] = 1 - x;
    }
    for (int p : prime) {
      if (1LL * p * x > n) break;
      vis[x * p] = 1;
      f[x * p] = x % p ? (1 - p) * f[x] : f[x];
      if (x % p == 0) break;
    }
  }
  for (int x = 1; x <= n; ++x) {
    f[x] = 1LL * x * f[x] % M + M;
    f[x] = (f[x] + f[x - 1]) % M;
  }
  long long res = 0;
  for (int l = 1, r; l <= n; l = r + 1) {
    int nn = n / l, mm = m / l;
    r = std::min(n / nn, m / mm);
    int g_n = (1LL * nn * (nn + 1) / 2) % M;
    int g_m = (1LL * mm * (mm + 1) / 2) % M;
    res += 1LL * g_n * g_m % M * (f[r] - f[l - 1] + M) % M;
  }
  std::cout << (res % M) << '\n';
  return 0;
}

The next example is special, requiring transformation of the divisor count function of a product.

LOJ 2185. [SDOI2015] Sum of Divisors

\(T\) groups of data. For each group, compute:

\[ \sum_{i=1}^n\sum_{j=1}^m\sigma_0(ij). \]

where \(\sigma_0(n)=\sum_{d \mid n}1\) denotes the number of divisors of \(n\).

Data range: \(1\le n,m,T\le 5\times 10^4\).

Solution

The difficulty lies in transforming \(\sigma_0(ij)\) into an expression involving gcd. Since \(\sigma_0\) is multiplicative, first consider prime power cases. For prime \(p\) and nonnegative exponents \(e_1,e_2\), let \(i=p^{e_1},~j=p^{e_2}\). Then

\[ \sigma_0(ij) = 1 + e_1 + e_2 = \sum_{x\mid i}\sum_{y\mid j}[x\perp y]. \]

For the general case, let \(i=\prod_p i_p\) and \(j=\prod_p j_p\), where \(i_p,j_p\) are the powers of \(p\) in the prime factorizations of \(i,j\) respectively. Then

\[ \sigma_0(ij) = \prod_p\sigma_0(i_pj_p)= \prod_p\sum_{x_p\mid i_p}\sum_{y_p\mid j_p}[x_p\perp y_p]. \]

Note that for each prime power factor \(i_p\) of \(i\), enumerating its divisors \(x_p\) is equivalent to enumerating divisors \(x\) of \(i\) and then extracting all prime power factors \(x_p\); similarly for \(j\). Thus, using distributive law,

\[ \sigma_0(ij) = \sum_{x\mid i}\sum_{y\mid j}\prod_p[x_p\perp y_p] = \sum_{x\mid i}\sum_{y\mid j}[x\perp y]. \]

The last step uses the fact: \(x\perp y\) if and only if for every prime factor \(p\), \(x_p\perp y_p\).

With this expression, apply the standard process:

\[ \begin{aligned} \sigma_0(ij) &= \sum_{x\mid i}\sum_{y\mid j}[x\perp y]\\ &= \sum_{x\mid i}\sum_{y\mid j}\sum_d\mu(d)[d\mid x][d\mid y]\\ &= \sum_d\mu(d)\left(\sum_{x}[d\mid x\mid i]\right)\left(\sum_{y}[d\mid y\mid j]\right)\\ &= \sum_d\mu(d)[d\mid i][d\mid j]\sigma_0\left(\dfrac{i}{d}\right)\sigma_0\left(\dfrac{j}{d}\right). \end{aligned} \]

The last step means: the function is nonzero only when \(d\mid i\) and \(d\mid j\), and then enumerating \(x\) with \(d\mid x\mid i\) is equivalent to enumerating divisors \(\dfrac{x}{d}\) of \(\dfrac{i}{d}\); similarly for \(y\).

Substitute back into the original sum and swap summation order:

\[ \begin{aligned} f(n,m) &= \sum_{i=1}^n\sum_{j=1}^m\sigma_0(ij)\\ &= \sum_{i=1}^n\sum_{j=1}^m\sum_d\mu(d)[d\mid i][d\mid j]\sigma_0\left(\dfrac{i}{d}\right)\sigma_0\left(\dfrac{j}{d}\right)\\ &= \sum_d\mu(d)\left(\sum_{i=1}^n[d\mid i]\sigma_0\left(\dfrac{i}{d}\right)\right)\left(\sum_{j=1}^m[d\mid j]\sigma_0\left(\dfrac{j}{d}\right)\right)\\ &= \sum_d\mu(d)\left(\sum_{i=1}^{\lfloor n/d\rfloor}\sigma_0(i)\right)\left(\sum_{j=1}^{\lfloor m/d\rfloor}\sigma_0(j)\right). \end{aligned} \]

Let \(G(n)=\sum_{i=1}^n\sigma_0(i)\). Then

\[ f(n,m)=\sum_{d}\mu(d)G\left(\left\lfloor\dfrac{n}{d}\right\rfloor\right)G\left(\left\lfloor\dfrac{m}{d}\right\rfloor\right). \]

This can be solved via number theoretic blocking. Only need to precompute prefix sums of \(\mu(n)\) and \(\sigma_0(n)\). Total time complexity is \(O(N+T\sqrt{N})\), where \(N\) is the upper bound of \(n,m\), and \(T\) is the number of groups.

Reference Code

#include <algorithm>
#include <iostream>
using namespace std;
constexpr long long N = 5e4 + 5;
long long n, m, T, pr[N], mu[N], d[N], t[N],
    cnt;  // t represents the count of the smallest prime factor of i
bool bp[N];

void prime_work(long long k) {
  bp[0] = bp[1] = true, mu[1] = 1, d[1] = 1;
  for (long long i = 2; i <= k; i++) {  // Linear sieve
    if (!bp[i]) pr[++cnt] = i, mu[i] = -1, d[i] = 2, t[i] = 1;
    for (long long j = 1; j <= cnt && i * pr[j] <= k; j++) {
      bp[i * pr[j]] = true;
      if (i % pr[j] == 0) {
        mu[i * pr[j]] = 0;
        d[i * pr[j]] = d[i] / (t[i] + 1) * (t[i] + 2);
        t[i * pr[j]] = t[i] + 1;
        break;
      } else {
        mu[i * pr[j]] = -mu[i];
        d[i * pr[j]] = d[i] << 1;
        t[i * pr[j]] = 1;
      }
    }
  }
  for (long long i = 2; i <= k; i++)
    mu[i] += mu[i - 1], d[i] += d[i - 1];  // Compute prefix sums
}

long long solve() {
  long long res = 0, mxi = min(n, m);
  for (long long i = 1, j; i <= mxi; i = j + 1) {  // Integer division block
    j = min(n / (n / i), m / (m / i));
    res += d[n / i] * d[m / i] * (mu[j] - mu[i - 1]);
  }
  return res;
}

int main() {
  cin.tie(nullptr)->sync_with_stdio(false);
  cin >> T;
  prime_work(50000);  // Preprocessing
  while (T--) {
    cin >> n >> m;
    cout << solve() << '\n';
  }
  return 0;
}

The final example demonstrates applying the multiplicative version of Möbius inversion.

Luogu P5221 Product

Compute:

\[ \prod_{i=1}^n\prod_{j=1}^n\dfrac{\operatorname{lcm}(i,j)}{\gcd(i,j)}\pmod{104857601}. \]

Data range: \(1\le n\le 1\times 10^6\).

Solution 1

Ignore modulus during derivation. Let

\[ f(n) = \prod_{i=1}^n\prod_{j=1}^n\dfrac{\operatorname{lcm}(i,j)}{\gcd(i,j)}. \]

Convert lcm to gcd:

\[ f(n) = \prod_{i=1}^n\prod_{j=1}^n\dfrac{ij}{(\gcd(i,j))^2}. \]

Note that the product factors are independent and can be computed separately. Let

\[ g(n) = \prod_{i=1}^n\prod_{j=1}^n\gcd(i,j). \]

Then the original expression equals:

\[ f(n) = \dfrac{(n!)^{2n}}{g(n)^2}. \]

The key is computing \(g(n)\). The process is similar to before but adapted to multiplication. First, enumerate and extract common factors:

\[ \begin{aligned} g(n) &= \prod_k\prod_{i=1}^n\prod_{j=1}^nk\uparrow[\gcd(i,j)=k]\\ &= \prod_k\prod_{i=1}^{\lfloor n/k\rfloor}\prod_{j=1}^{\lfloor n/k\rfloor}k\uparrow[\gcd(i,j)=1]. \end{aligned} \]

Here, \(a\uparrow b=a^b\) is the Knuth arrow. Then substitute \([\gcd(i,j)=1]=\sum_d\mu(d)[d\mid i][d\mid j]\), and convert the sum in the exponent to a product of powers:

\[ g(n) = \prod_k\prod_d\prod_{i=1}^{\lfloor n/k\rfloor}\prod_{j=1}^{\lfloor n/k\rfloor}k\uparrow(\mu(d)[d\mid i][d\mid j]). \]

Further extract factors (set \(i=di'\), \(j=dj'\)) and use properties of the floor function:

\[ g(n) = \prod_k\prod_d\prod_{i=1}^{\lfloor n/(kd)\rfloor}\prod_{j=1}^{\lfloor n/(kd)\rfloor}k\uparrow\mu(d). \]

Separating the product over \(i,j\), we find they do not appear, so it's equivalent to raising the product to a power:

\[ g(n) = \prod_k\prod_d k\uparrow\left(\mu(d)\left\lfloor\dfrac{n}{kd}\right\rfloor^2\right). \]

Since we enumerated common divisors, we need to swap product order again. Let \(\ell = kd\):

\[ \begin{aligned} g(n) &= \prod_{\ell}\prod_{d\mid\ell}\left(\dfrac{\ell}{d}\right)\uparrow\left(\mu(d)\left\lfloor\dfrac{n}{\ell}\right\rfloor^2\right)\\ &= \prod_\ell\left(\prod_{d\mid\ell}\left(\dfrac{\ell}{d}\right)\uparrow\mu(d)\right)\uparrow\left\lfloor\dfrac{n}{\ell}\right\rfloor^2. \end{aligned} \]

Let

\[ F(n) = \prod_{d\mid n}\left(\dfrac{n}{d}\right)\uparrow\mu(d). \]

This is easily recognized as the multiplicative Möbius inversion of \(\tilde F(n)=n\). Even without knowing its expression, it can be preprocessed in \(O(n\log\log n)\) time using Dirichlet difference. Of course, since \(\tilde F(n)\) is simple, \(F(n)\) can be directly derived:

\[ F(n) = \begin{cases} p, & n = p^e,~p\in\mathbf P,~e\in\mathbf N_+, \\ 1, &\text{otherwise}. \end{cases} \]

The von Mangoldt function is its natural logarithm. After obtaining \(F(n)\), applying multiplicative number theoretic blocking yields \(g(n)\) in \(O(\sqrt{n})\) time, and thus \(f(n)\). Total time complexity is \(O(n)\).

Note: when dealing with products, Euler's theorem is often needed, so the modulus for exponentiation differs from the given modulus.

Solution 2

The difficulty with multiplicative derivation is handling products and exponents. Alternatively, take logarithms first. For this problem, consider only \(g(n)\). Taking logarithm:

\[ \log g(n) = \sum_{i=1}^n\sum_{j=1}^n\log\gcd(i,j). \]

For such expressions involving gcd, apply the standard process:

\[ \begin{aligned} \log g(n) &= \sum_k\log k\sum_{i=1}^n\sum_{j=1}^n[\gcd(i,j)=k]\\ &= \sum_k\log k\sum_{i=1}^{\lfloor n/k\rfloor}\sum_{j=1}^{\lfloor n/k\rfloor}[\gcd(i,j)=1]\\ &= \sum_k\log k\sum_d\mu(d)\left(\sum_{i=1}^{\lfloor n/k\rfloor}[i\mid d]\right)\left(\sum_{j=1}^{\lfloor n/k\rfloor}[j\mid d]\right)\\ &= \sum_k\log k\sum_d\mu(d)\left\lfloor\dfrac{n}{kd}\right\rfloor^2\\ &= \sum_{\ell}\left(\sum_d\mu(d)\log\dfrac{\ell}{d}\right)\left\lfloor\dfrac{n}{\ell}\right\rfloor^2\\ &= \sum_{\ell}\Lambda(\ell)\left\lfloor\dfrac{n}{\ell}\right\rfloor^2. \end{aligned} \]

where \(\Lambda(n)\) is the von Mangoldt function. Exponentiating this yields the result from Solution 1.

Reference Code

#include <algorithm>
#include <iostream>
#include <vector>

int pow(int a, int b, int m) {
  int res = 1;
  for (int po = a; b; b >>= 1) {
    if (b & 1) res = 1LL * res * po % m;
    po = 1LL * po * po % m;
  }
  return res;
}

int main() {
  constexpr int M = 104857601;
  int n;
  std::cin >> n;
  // Get F(n), i.e., exp(Lambda(n)).
  std::vector<int> vis(n + 1), prime, pf(n + 1), rem(n + 1);
  prime.reserve(n);
  for (int x = 2; x <= n; ++x) {
    if (!vis[x]) {
      prime.emplace_back(x);
      pf[x] = x;
      rem[x] = 1;
    }
    for (int p : prime) {
      if (1LL * p * x > n) break;
      vis[x * p] = 1;
      pf[x * p] = p;
      rem[x * p] = x % p ? x : rem[x];
      if (x % p == 0) break;
    }
  }
  pf[1] = 1;
  for (int x = 2; x <= n; ++x) {
    pf[x] = rem[x] == 1 ? pf[x] : 1;
  }
  // Prefix products and their inverses.
  std::vector<int> pm(n + 1), ip(n + 1);
  pm[0] = 1;
  for (int x = 1; x <= n; ++x) {
    pm[x] = 1LL * pm[x - 1] * pf[x] % M;
  }
  int inv = pow(pm[n], M - 2, M);
  ip[0] = 1;
  for (int x = n; x >= 1; --x) {
    ip[x] = inv;
    inv = 1LL * inv * pf[x] % M;
  }
  // Calculate g(n).
  int res = 1;
  for (int l = 1, r; l <= n; l = r + 1) {
    r = n / (n / l);
    int a = 1LL * pm[r] * ip[l - 1] % M;
    int b = 1LL * (n / l) * (n / l) % (M - 1);
    res = 1LL * res * pow(a, b, M) % M;
  }
  // Take square and then inverse.
  res = pow(res, M - 3, M);
  // Get factorials.
  int fac = 1;
  for (int x = 1; x <= n; ++x) {
    fac = 1LL * fac * x % M;
  }
  // Final answer.
  res = 1LL * res * pow(fac, 2 * n, M) % M;
  std::cout << res << std::endl;
  return 0;
}

Exercises

References

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