「SDIO2015」约数个数和 - 莫比乌斯反演

发表于 更新于 分类于 ， 阅读次数： Disqus：

题意描述

洛谷链接

LibreOJ 链接

设 \(d(x)\) 为 \(x\) 的约数个数，给定 \(n, m\)，求

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

对于所有的数据，\(1 \leq N, M \leq 50000,\ 1 \leq T \leq 50000\)。

解题思路

约数个数的变换

对这个式子，我们首先要解决的是 \(d(ij)\)。对于 \(d(ij)\)，我们有以下结论：

\[ d(ij) = \sum_{x | i}\sum_{y | j}[\gcd(x, y) = 1] \]

证明：我们可以对 \(i\)，\(j\) 分解质因数：\(i = \prod{ p_i^{a_i} }\)，\(j = \prod{ p_i^{b_i} }\)，故 \(ij = \prod{ p_i^{c_i = a_i + b_i} }\)。对于 \(ij\) 的每个约数中的 \(p_i\)，我们总可以找到一个 \(d_i (d_i \le c_i)\) 进行对应。而对于每个 \(d_i\)，我们分为两种情况：

  1. \(d_i <= a_i\)，则 \(p_i\) 全部在 \(i\) 中取；
  2. \(a_i < d_i \le a_i + b_i\)，则 \(p_i\) 在 \(i\) 中取完后再在 \(b_i\) 中补充。

我们可以发现，对于每一种 \(d_i\)，用上述方法我们总可以找到唯一一种用 \(a_i\)、\(b_i\) 对应的方法；所以对应每一种 \(ij\) 的约数，用上述方法我们总可以找到对应且为唯一对应。

于是我们可以利用 \(i\)、\(j\) 的约数进行统计。令 \(x | i\)，\(y | j\)；当 \(\gcd(i, j) = 1\) 时，对于 \(p_i\)，\(i\)、\(j\) 中仅有一个含有 \(p_i\)：

  • 当仅有 \(x\) 含有 \(p_i\) 时，对应上述方法 \(1\)；
  • 当仅有 \(y\) 含有 \(p_i\) 时，其中「\(x\) 不含有 \(p_i\)，\(y\) 含有 \(p_i\)」所含情况数和「\(i\) 中的 \(p_i\) 被取完，\(y\) 含有 \(p_i\)」相同，对于 \(d(ij)\) 只统计个数来说是等价的。所以对应上述方法 \(2\)。

综上，我们仅需统计满足 \(\gcd(i, j) = 1\) 的 \((i, j)\) 的个数即可。即 \(d(ij) = \sum_{x | i}\sum_{y | j}[\gcd(x, y) = 1]\)。

莫比乌斯反演

接下来原式就变成了下列形式：

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

我们设 \(A(n, m) = \sum_{i = 1}^n \sum_{j = 1}^m\sum_{x | i}\sum_{y | j}[\gcd(x, y) = 1]\)，显然 \(A(n, m)\) 即为答案。

对于 \([\gcd(x, y) = 1]\) 其中的 \(1\) 不好化简，于是我们可以设下列函数：

\[ \begin{align*} f(a) &= \sum_{i = 1}^n\sum_{j = 1}^m\sum_{x | i}\sum_{y | j}[\gcd(x, y) = a] \\ &= \sum_{x = 1}^n\sum_{y = 1}^m \lfloor \frac{n}{x} \rfloor \lfloor \frac{m}{y} \rfloor [\gcd(x, y) = a] \end{align*} \]

显然 \(A(n, m) = f(1)\)。接下来我们要化简 \(f(a)\)。我们设下列函数：

\[ \begin{align*} F(b) &= \sum_{b | a}f(a) \\ &= \sum_{b | a}\sum_{x = 1}^n\sum_{y = 1}^m \lfloor \frac{n}{x} \rfloor \lfloor \frac{m}{y} \rfloor [\gcd(x, y) = a] \\ &= \sum_{t = 1}^{ \min\{ \lfloor \frac{n}{a} \rfloor \lfloor \frac{m}{a} \rfloor \} }\sum_{x = 1}^n\sum_{y = 1}^m \lfloor \frac{n}{x} \rfloor \lfloor \frac{m}{y} \rfloor [\gcd(x, y) = a] \\ &= \sum_{x = 1}^n\sum_{y = 1}^m \lfloor \frac{n}{x} \rfloor \lfloor \frac{m}{y} \rfloor [b | \gcd(x, y)] \\ &= \sum_{x = 1}^n\sum_{y = 1}^m \lfloor \frac{n}{x} \rfloor \lfloor \frac{m}{y} \rfloor [b | x][b | y] \\ &= \sum_{x = 1}^n\sum_{y = 1}^m \lfloor \frac{n}{x} \rfloor \lfloor \frac{m}{y} \rfloor [b | x][b | y] \\ &= \sum_{i = 1}^{ \lfloor\frac{n}{b}\rfloor }\sum_{j = 1}^{ \lfloor\frac{m}{b}\rfloor } \lfloor\frac{n}{bi}\rfloor\lfloor\frac{m}{bj}\rfloor \\ \end{align*} \]

由于 \(F(b) = \sum_{b | a}f(a)\)，我们可以进行莫比乌斯反演：

\[ f(a) = \sum_{a | b}\mu(\frac{b}{a})F(b) \]

接下来可以带入 \(A(n, m)\)：

\[ \begin{align*} A(n, m) &= f(1) \\ &= \sum_{b}^{ \min\{n, m\} }\mu(b)F(b) \\ &= \sum_{b}^{ \min\{n, m\} }\mu(b)\sum_{i = 1}^{ \lfloor\frac{n}{b}\rfloor }\sum_{j = 1}^{ \lfloor\frac{m}{b}\rfloor } \lfloor\frac{n}{bi}\rfloor\lfloor\frac{m}{bj}\rfloor \\ &= \sum_{b}^{ \min\{n, m\} }\mu(b)(\sum_{i = 1}^{ \lfloor\frac{n}{b}\rfloor }\lfloor\frac{n}{bi}\rfloor)(\sum_{j = 1}^{ \lfloor\frac{m}{b}\rfloor } \lfloor\frac{m}{bj}\rfloor) \\ &= \sum_{b}^{ \min\{n, m\} }\mu(b)(\sum_{i = 1}^{ \lfloor\frac{n}{b}\rfloor }\lfloor\frac{\lfloor\frac{n}{b}\rfloor}{i}\rfloor)(\sum_{j = 1}^{ \lfloor\frac{m}{b}\rfloor }\lfloor\frac{\lfloor\frac{m}{b}\rfloor}{j}\rfloor) \end{align*} \]

接下来预处理一下 \(\mu(n)\) 的前缀和与 \(\sum_{i = 1}^n\lfloor\frac{n}{i}\rfloor\) 的前缀和，在询问时用数论分块求 \(A(n, m)\) 即可。预处理复杂度 \(O(n\sqrt{n})\)，单次询问复杂度 \(O(\sqrt{n})\)。

代码演示

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
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
#include <cstdio>
#include <iostream>

const int MAXN = 5e4;

bool isNotPrime[MAXN + 1];
int mu[MAXN + 1], primes[MAXN + 1], cnt;
long long sumMu[MAXN + 1], sum[MAXN + 1];

inline void getPrimes() {
    isNotPrime[0] = isNotPrime[1] = true;
    mu[1] = 1;
    for (int i = 2; i <= MAXN; i++) {
        if (!isNotPrime[i]) {
            primes[++cnt] = i;
            mu[i] = -1;
        }

        for (int j = 1; j <= cnt; j++) {
            int t = i * primes[j];
            if (t > MAXN) break;

            isNotPrime[t] = true;

            if (i % primes[j] == 0) {
                mu[t] = 0;
                break;
            } else mu[t] = -mu[i];
        }
    }

    sumMu[0] = 0;
    for (int i = 1; i <= MAXN; i++) sumMu[i] = sumMu[i - 1] + mu[i];
    for (int i = 1; i <= MAXN; i++) {
        for (int l = 1, r = 0; l <= i; l = r + 1) {
            r = i / (i / l);
            sum[i] += (long long)(r - l + 1) * (i / l);
        }
    }
}

inline long long f(int n, int m) {
    long long ans = 0;
    for (int l = 1, r = 0; l <= std::min(n, m); l = r + 1) {
        r = std::min(n / (n / l), m / (m / l));
        ans += (sumMu[r] - sumMu[l - 1]) * sum[n / l] * sum[m / l];
    }
    return ans;
}

int main() {
    getPrimes();

    int t;

    scanf("%d", &t);

    while (t--) {
        int n, m;

        scanf("%d %d", &n, &m);

        printf("%lld\n", f(n, m));
    }

    return 0;
}