Möbius Inversion

For a positive integer $n$,

Here $[\cdot]$ is the Iverson bracket.

Let $n=\prod_{i=1}^kp_i^{e_i}$ and let $n' = \prod_{i=1}^kp_i$. By the binomial theorem,

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

Direct verification gives:

The key steps are exchanging the order of summation and observing that $k\mid d\mid n$ is equivalent to $\dfrac{n}{d}\mid\dfrac{n}{k}$. The penultimate equality sums the Möbius function over the divisors $\dfrac{n}{k}$ of $\dfrac{n}{d}$, so it equals $\left[\dfrac{n}{k} = 1\right]$. This expression is nonzero only when $n=k$; in that case, it is not $0$, and the result is $g(n)$.

Using Dirichlet convolution, the proposition is equivalent to

Using $1 * \mu = \varepsilon$, convolve both sides of the equality on the left with $\mu$; this gives

For a prime power $n=p^e~(e\in\mathbf N_+)$, there is

For $n=1$, clearly $\Lambda(n)=\log 1=0$. For other composite $n$,

By the properties of the Möbius function, the coefficient of the term $\log n$ is $[n=1]=0$. For the latter term, factor $d$ further into primes. For any prime $p\mid n$, consider the coefficient of $\log p$; it is:

Thus every composite $n$ with more than one prime factor satisfies $\Lambda(n)=0$.

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

Direct verification gives:

This derivation is completely dual to that of the basic form.

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

Direct verification gives:

Here $a\uparrow b = a^b$ is the Knuth arrow. Comparing this proof with the basic one, the only differences are replacing addition by multiplication and multiplication by exponentiation.

Let $f(n),g(n),\alpha(n)$ be arithmetic functions, and suppose $\alpha^{-1}(n)$ is the Dirichlet inverse of $\alpha(n)$, namely

Then

Direct verification gives:

Compared with the proof of the basic form, only the penultimate equality needs to be replaced by the definition of the Dirichlet inverse.

Let $f(n),g(n)$ be arithmetic functions, and let $t(n)$ be completely multiplicative. Then

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

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

Extend $F$ and $G$ to the region $x < 1$, where they satisfy $F(x)=G(x)=0$. Then the proposition is equivalent to:

All these sums range over $n\in\mathbf N_+$.

Direct verification gives:

To obtain the second equality, set $k = nd$.

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

It suffices to take $F(x)=f(\lfloor x\rfloor)$ and $G(x)=g(\lfloor x\rfloor)$.

By inclusion–exclusion, the original expression can be split into $4$ blocks, each having the form

For expressions of this kind, the standard derivation is: factor out common divisors, apply properties of the Möbius function, and exchange the order of summation.

Since $i,j$ can only be multiples of $k$, factor this out first—equivalently substitute $i=ki'$ and $j=kj'$, obtaining:

By the Möbius-function identity,

Substituting this and exchanging the order of summation gives:

The benefit is that, with $d$ fixed, the terms involving $i$ and $j$ separate and can be summed independently. Next, since

Therefore,

After precomputing $\mu(d)$ with a linear sieve and its prefix sums, number-theoretic block decomposition solves the problem. The total time complexity is $O(N + T\sqrt{N})$, where $N$ is the upper bound for $n,m$, and $T$ is the number of test cases.

#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 above
  }
  return res;
}

int main() {
  cin.tie(nullptr)->sync_with_stdio(false);
  int T, a, b, c, d, k;
  init();  // Precompute the mu array
  cin >> T;
  for (int i = 1; i <= T; i++) {
    cin >> a >> b >> c >> d >> k;
    // By inclusion–exclusion, subtract the count in 1<=x<=b&&1<=y<=c-1 from the count in 1<=x<=b&&1<=y<=d, and
    //   add the count in 1<=x<=a-1&&1<=y<=c-1 to the count in 1<=x<=a-1&&1<=y<=d,
    //   obtaining the count in a<=x<=b&&c<=y<=d.
    // If this step is unclear, draw a coordinate diagram to understand it
    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 problem gives a least common multiple, but a greatest common divisor is often easier to handle. First transform it:

Factor out $n$ and enumerate the greatest common divisor $k$:

For the inner sum, this is the standard greatest-common-divisor case; repeating the usual procedure gives:

Again, the sum over $i$ separates from the other terms and can be handled independently. The final sum is an arithmetic-progression sum (take $i=di'$):

This yields:

After enumerating common divisors, double sums of this form are common. Use the standard method: set the product to a new variable $\ell=kd$, then exchange the order of summation again. Since $d\mid(n/k)$ is equivalent to $d\mid\ell\mid n$, the expression becomes:

Let $F(\ell)=\sum_{d\mid\ell}\mu(d)d$; then the expression has the form:

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

Both terms (without their coefficients) are multiplicative and can be precomputed directly with a linear sieve (or the inner function can be sieved first and then preprocessed with a Dirichlet prefix sum in $O(N\log\log N)$ time). Specifically, let

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

In particular, $H_1(p^e)\equiv 1$ is the constant function, while

This is easy to solve with a linear sieve. After precomputing $H_2(n)$, a query can be answered using $f(n)=(n/2)(H_2(n)+1)$ in $O(1)$ time. The total time complexity is $O(N+T)$, where $N$ is the upper bound for $n$ and $T$ is the number of test cases.

The reference implementation derives the linear-sieve part further using the special form of this problem; this is not necessary. Using only the values at prime powers, preprocessing can still be completed in $O(N)$ time. The derivation is given in Solution 2.

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

Instead of continuing Möbius inversion here, observe that the following sum is actually the sum of integers not exceeding $d=n/k$ and coprime to it. For $d>1$, the integers not exceeding $d$ that are coprime to it occur in pairs: $i$ and $d-i$ are both coprime to $d$, so

For $d=1$,

Therefore, the original expression can be written as

Since $G(n)=\sum_{d\mid n}d\varphi(d)$ is the Dirichlet convolution of the multiplicative function $n\varphi(n)$ with the constant function $1$, it is multiplicative and can be precomputed with a linear sieve. It remains to determine its values at prime powers. For a prime $p$ and positive exponent $e$,

This agrees with the result derived in Solution 1. The total time complexity remains $O(N+T)$.

Finally, the expression of this multiplicative function permits a further optimization of the linear sieve. For a prime $p$,

For a general $n$, the key to the linear sieve is computing $G(pn)$. This splits into two cases. When $p\perp n$, multiplicativity of $G$ gives

Otherwise, when $p\mid n$, write $n=p^em$ with $p\perp m$; then

Direct verification shows that this expression also holds when $p\perp n$. Therefore,

Substituting into the equation above gives

This simplifies the linear-sieve computation. The derivation is not necessary: for a multiplicative function without special properties, $G(pn)=G(p^{e+1})G(m)$ directly completes the linear sieve.

#include <iostream>
constexpr int N = 1000000;
int tot, p[N + 5];
long long g[N + 5];
bool flg[N + 5];  // Marker 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];  // Formula derived above
        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();  // Precompute the g array
  cin >> T;
  for (int i = 1; i <= T; ++i) {
    cin >> n;
    cout << (g[n] + 1) * n / 2 << '\n';
  }
  return 0;
}

Ignore the modulus during the derivation. Let

Again convert the least common multiple to a greatest common divisor, enumerate the divisor, and apply the standard procedure to obtain

Again, the sums over $i$ and $j$ separate. First compute these inner sums by extracting a factor (take $i=di'$); then

Here $G(n)=\dfrac{1}{2}n(n+1)$ is the arithmetic-progression sum, and the final equality uses a property of the floor function. The other sum is symmetric. Substituting into the previous expression gives

As before, expressions enumerating common divisors usually require enumerating the product $\ell = kd$ and exchanging the order of summation again:

Let

This is the product of the multiplicative function $\ell$ and the multiplicative function $\sum_{d\mid\ell}\mu(d)d$, so it is multiplicative. Precompute it and its prefix sums with a linear sieve, then use number-theoretic block decomposition to compute $f(n,m)$. The total time complexity is $O(\min\{n,m\})$.

#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 difficulty is converting $\sigma_0(ij)$ into an expression involving greatest common divisors. Since $\sigma_0$ is multiplicative, first consider prime powers. For a prime $p$ and nonnegative exponents $e_1,e_2$, let $i=p^{e_1},~j=p^{e_2}$; then

In general, write $i=\prod_p i_p$ and $j=\prod_p j_p$, where $i_p,j_p$ are the powers in the prime factorizations of $i,j$, respectively, corresponding to the prime $p$. Then

Notice that, for each $i$, enumerating its prime-power factors $i_p$ and their divisors $x_p$ is equivalent to enumerating $i$'s divisors $x$ and then decomposing them into all prime-power factors $x_p$; the same holds for $j$. By distributivity, this gives

conclusion：$x\perp y$，if and only ifForeach $p$， $x_p\perp y_p$ ．

obtain ， canApplications ：

：function $d\mid i$ $d\mid j$ nonzerovalue， ， satisfy $d\mid x\mid i$ $x$ When $\dfrac{i}{d}$ $\dfrac{x}{d}$， satisfy $d\mid y\mid j$ $y$ ．

， findsum ：

Let $G(n)=\sum_{i=1}^n\sigma_0(i)$，

canbynumber theoryblock decompositionfind ． need to $\mu(n)$ sum $\sigma_0(n)$ prefix sum ． time complexity $O(N+T\sqrt{N})$，where，$N$ $n,m$ ，$T$ ．

#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 denotes the number of times the smallest prime factor of i occurs.
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) {  // Division into quotient blocks.
    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);  // Precompute.
  while (T--) {
    cin >> n >> m;
    cout << solve() << '\n';
  }
  return 0;
}