..

数论学习笔记

基本定义

对于两个自然数 $a, b$，如果 $a = kb$ 且 $b \neq 0$，那么我们说 $b$ 整除 $a$ 或者说 $a$ 被 $b$ 整除，记作 $b \mid a$；否则 $b$ 不整除 $a$，记作 $b \nmid a$。如果 $b \mid a$，那么我们说 $b$ 是 $a$ 的因数，$a$ 是 $b$ 的倍数。如果一个数同时是两个数的因数，那么这个数是这两个数的公因数。如果两个数除了 $1$ 之外没有别的公因数，那么这两个数互质，记作 $a \ \bot \ b$。

欧几里得算法

欧几里得算法又称辗转相除法，用来求两个数的最大公因数。依据是 $\gcd(a, b) = \gcd(b, a\ \mathrm{mod} \ b)$。

证明：

$\gcd(a, b) = \gcd(b, a\ % \ b) \iff \gcd(a, b) = \gcd(b, a - b)$。使用反证法证明。设 $r = \gcd(a, b)$，则 $r$ 也是 $b, a - b$ 的因数。如果 $r \neq \gcd(b, a - b)$，那么一定有一个数 $s$ 满足 $s = \gcd(b, a - b)$。此时 $s \mid b$ 且 $s \mid a - b$，故 $s \mid a$，故 $s$ 也是 $a, b$ 的公因数。但 $s \gt r = \gcd(a, b)$，产生矛盾，故不可能有符合条件的 $s$，故 $r = \gcd(b, a - b)$。

证毕。

代码：

int gcd(int a, int b) {
    return b == 0 ? a : gcd(b, a % b);
}

时间复杂度为对数级别。

关于最大公约数还有一个性质是 $\gcd(a, b) \times \operatorname{lcm}(a, b) = a \times b$。

证明：

根据算术基本定理，$a, b$ 均可表示成如下形式：

$\begin{aligned} a = p_1^{a_1}p_2^{a_2} \ldots p_n^{a_n}\ b = p_1^{b_1}p_2^{b_2} \ldots p_m^{b_m}\ \end{aligned}$

其中 $p_i$ 为质数。我们将 $a$ 和 $b$ “对齐”，也就是将他们的质因数逐个对应，只在一个数里出现过的令其指数为 $0$。显然

$\begin{aligned} \gcd(a, b) = \prod_{i = 1}^{\max(n, m)}p_i^{\min(a_i, b_i)}\ \operatorname{lcm}(a, b) = \prod_{i = 1}^{\max(n, m)}p_i^{\max(a_i, b_i)} \end{aligned}$

故显然 $\gcd(a, b) \times \operatorname{lcm}(a, b) = a \times b$。

扩展欧几里得算法

扩展欧几里得算法用来求解形如 $ax + by = \gcd(a, b)$ 的不定方程。 $$\begin{aligned} ax + by &= \gcd(a, b) \\[2ex] &= \gcd(b, a \ \% \ b) \\[2ex] \end{aligned}$$ 因此我们可以列出新的方程

$$\begin{aligned} bx’ + (a \ \% \ b)y’ &= \gcd(b, a \ \% \ b) \\[2ex] &= bx’ + (a - \lfloor \frac{a}{b} \rfloor \times b) y’ \\[2ex] &= ay’ + b \cdot (x’ - \lfloor \frac{a}{b} \rfloor \times y’) \end{aligned}$$

故 $x = y’, y = x’ - \lfloor \frac{a}{b} \rfloor \times y’$。所以我们可以在进行欧几里得算法的同时递归得出每一组 $a, b$ 所对应的 $x, y$ 。递归终止条件是当 $b = 0$ 时 $x = 1, y = 0$。其实 $y$应该是可以随便取的，不过大家都取 $0$（可能是为了方便计算）。这样便可得到不定方程的一组解。已经求出一组 $x, y$ 以后，$ax + by = a(x - b) + b(y + a)$，所以我们可以由一组解计算出其他解。

代码：

int exgcd(int a, int b, int &x, int &y) {
    if (b == 0) {
        x = 1, y = 0;
        return a;
    }
    int res = exgcd(b, a % b, y, x);
    y -= x * a / b;
    return res;
}

扩展欧几里得算法还可以用来计算 $ax \equiv 1 \pmod b$ 的同余方程。因为这种方程可以化为 $ax + by = 1$，当 $\gcd(a, b) = 1$ 即 $a, b$ 互质的时候此方程有解，其他情况无解。

质数

质数的定义不再重复。

算术基本定理

任意一个正整数都可以表示成如下形式： $n = p_1^{a_1}p_2^{a_2} \ldots p_n^{a_n}$ ，其中 $p$ 是质数。

质因数分解

分解 $n$ 的质因数的方法是，枚举所有小于等于 $\sqrt{n}$ 的数，如果这个数是 $n$ 的因数，那么就把这个数从 $n$ 中除去并记录它的指数增加了一；如果除完后这个数仍能整除 $n$，那么继续除到无法整除为止。枚举完这些数以后如果 $n$ 被除后剩余的商仍不等于一，那么这个剩余的数就是 $n$ 的唯一一个大于 $\sqrt{n}$ 的质因数。这个算法的正确性十分显然，不再证明。

根据乘法原理，如果 $n = p_1^{a_1}p_2^{a_2} \ldots p_n^{a_n}$，那么 $n$ 的因数个数为 $\prod_{i=1}^n (a_i+1)$。

质数筛法

这里只记录线性筛相关知识，暴力筛和埃氏筛不再记录。

线性筛的基本思想是，每个合数只会被它的最小质因数筛去。

代码：

bool isPrime[MAXM];
int primes[7300], cnt, sum[MAXM];

void sieve(int m) {
    std::fill(isPrime + 1, isPrime + m + 1, true);
    isPrime[0] = isPrime[1] = false;
    for (int i = 2; i <= m; ++i) {
        if (isPrime[i]) {
            primes[++cnt] = i;
        }
        for (int j = 1; j <= cnt && i * primes[j] <= m; ++j) {
            isPrime[i * primes[j]] = false;
            if (i % primes[j] == 0) { // 线性之处
                break;
            }
        }
    }
}

线性筛还可以用来预处理出一系列数论函数的值。

逆元

因为除法在取模意义下不封闭，因此如果我们想在取模意义下做“除法”，就要用到逆元。

如果 $a \cdot a^{-1} \equiv 1 \pmod p$，那么我们就将 $a^{-1}$ 记作 $a$ 在 模 $p$ 意义下的逆元。

求逆元的方法有若干种。

费马小定理

如果 $p$ 是质数，那么 $a^{p - 1} \equiv 1 \pmod p$。所以 $a \cdot a^{p-2} \equiv 1 \pmod p$，即 $a$ 在模 $p$ 意义下的逆元是 $a^{p - 2}$。可以用快速幂在 $O(\log n)$ 的时间复杂度内求出 $a$ 的逆元。

线性求逆元

考虑取模的定义式： $p \bmod a = p - \lfloor \frac{p}{a} \rfloor \times a$ 将这个式子放在模 $p$ 同余系下： $p \bmod a \equiv -\lfloor \frac{p}{a} \rfloor \times a \pmod p$ 两边同时乘以 $(p \bmod a)^{-1} \cdot a^{-1}$： $a^{-1} \equiv -\lfloor \frac{p}{a} \rfloor \times (p \bmod a)^{-1} \pmod p$ 因为 $(p \bmod a)^{-1} \lt a$，因此逆元可以线性递推得出。初始条件为 $1^{-1} = 1$。写代码的时候需注意负数取模。

代码：

int inv[MAXN];

void get_inv(int n, int p) {
    inv[1] = 1;
    for (int i = 2; i <= n; ++i) {
        inv[i] = 1LL * (p - p / i) * inv[p % i] % p;
    }
}

积性函数与线性筛

积性函数

设 $f$ 是一个定义域为 $\mathrm{N}_+$ 的函数，如果对于任意两个互质的正整数 $a, b$ 都满足 $f(a \cdot b) = f(a) \cdot f(b)$，那么称 $f$ 为积性函数。如果对于任意两个正整数（不要求互质）都满足 $f(a \cdot b) = f(a) \cdot f(b)$，那么称 $f$ 为完全积性函数。

积性函数的性质

  1. 对于任意积性函数 $f$ 都有 $f(1) = 1$。证明略。

  2. 对于一个大于 $1$ 的正整数 $n$，设$n = \prod p_i^{a_i}$，其中 $p_i$ 为互不相同的素数，那么对于积性函数 $f$ 有

    $f(n) = f(\prod p_i^{a_i}) = \prod f(p_i^{a_i})$

    对于完全积性函数 $f$ 还有

    $f(n) = \prod f(p_i)^{a_i}$

    证明略。

常见的积性函数

欧拉函数

定义欧拉函数 $\varphi(n)$ 为小于 $n$ 的正整数中与 $n$ 互质的数的个数。

根据容斥原理可以推出：

若

$n = \prod_{i = 1}^{m} p_i^{a_i}$

，其中 $p_i$ 为质数，那么

$\varphi(n) = n \prod_{i = 1}^{m} (1 - \frac{1}{p_i})$

。证明略。（是因为我不会证）

欧拉函数的性质

  1. 欧拉函数是积性函数，但不是完全积性函数。

    证明：设 $a, b$ 为两个互质的正整数，则

    $\begin{aligned} \varphi(a) &= a \prod (1 - \frac{1}{p_{a_i}}) \ \varphi(b) &= b \prod (1 - \frac{1}{p_{b_i}}) \ \varphi(a) \varphi(b) &= a \prod p_{a_i} \cdot b \prod p_{b_i} \ &= ab \prod (1 - \frac{1}{p_{a_i}})(1 - \frac{1}{p_{b_i}}) \end{aligned}$

    因为 $a, b$ 互质，所以 $p_{a_1}, p_{b_2}$ 互不相同，所以 $\varphi(a)\varphi(b) = \varphi(ab)$。显然若 $a, b$ 不互质则 $\varphi(a)\varphi(b) \neq \varphi(ab)$，因此欧拉函数不是完全积性函数。

  2. 设 $p$ 为质数， $k$ 为正整数，那么 $\varphi(p^k) = p^k - p^{k - 1}$。

    证明：我们考虑枚举小于 $p^k$ 的正整数中与 $p^k$ 不互质的数。显然这样的数都是 $p$ 的倍数。它们是：

    $1 \times p, 2 \times p, \cdots, (p^{k-1} - 1) \times p$

    所以 $\varphi (p^k)$ 的值为小于 $p^k$ 的正整数的个数减去这些数中与 $p^k$ 互质的数的个数。即

    $(p^k - 1) - (p ^{k - 1} - 1) = p^k - p^{k - 1}$

    也可以从定义式直接证明。

  3. 欧拉定理：设 $a, n$ 为两个互质的正整数，那么 $a^{\varphi(n)} \equiv 1 \pmod n$。

    证明：设 $x_1, x_2, \cdots, x_{\varphi(n)}$ 是小于 $n$ 的数中与 $n$ 互质的数，显然这些数有 $\varphi(n)$ 个。令 $m_i = a \cdot x_i$。

    引理 1 ：$m_i \bmod n$ 的值互不相同。证明：设 $m_s \gt m_r$。若 $m_s \equiv m_r \pmod n$，则有 $m_s - m_r = a(x_s - x_r) = kn$，即 $n \mid a(x_s - x_r)$。但 $a \ \bot \ n$，而 $x_s - x_r \le n$，故 $n \nmid a(x_s - x_r)$，故不可能有两个 $m_i$ 模 $n$ 同余。所以 $\varphi(n)$ 个数有 $\varphi(n)$ 个不同的余数。

    引理 2：$(m_i \bmod n) \ \bot \ n$。证明：设 $r$ 为 $m_i \bmod n$ 与 $n$ 的公因数，则 $ax_i = m_i = pn + qr$，$ax_i$ 与 $n$ 不互质。但 $a \ \bot \ n$， $x_i \ \bot \ n$，所以 $ax_i \ \bot \ n$，故 $m_i \bmod n$ 与 $n$ 不可能有大于 $1$ 的公因数。

    由此可得 $m_i$ 模 $n$ 的余数是小于 $n$ 的数中与 $n$ 互质的数，即 ${m_i \bmod n} = {x_i}$。

    故

    $\begin{aligned} &\hspace{6ex}\prod m_i \equiv \prod x_i &\pmod n \ &\implies a^{\varphi(n)} \prod x_i \equiv \prod x_i &\pmod n \ &\implies (a^{\varphi(n)} - 1) \prod x_i \equiv 0 &\pmod n \ \end{aligned}$

    因为 $x_i$ 与 $n$ 互质，所以

    $a^{\varphi(n)} - 1 \equiv 0 \pmod n \implies a^{\varphi(n)} \equiv 1 \pmod n$

    证毕。

    费马小定理是当 $n$ 是质数的时候的欧拉定理的一种特殊情况。

  4. 对于正整数 $n$，$\sum_{d \mid n} \varphi(d) = n$。

    证明：当 $n = 1$ 时，原式显然成立。

    当 $n = p^k$，其中 $p$ 是质数时，

    $\begin{aligned} \sum_{d \mid n} \varphi(d) &= \sum_{i = 0}^{k} \varphi(p^i) \ &= 1 + \sum_{i = 1}^k \varphi(p^i) \ &= 1 + \sum_{i = 1}^k (-p^{i - 1} + p^i) \ &= p^i \end{aligned}$

    当 $n = p_1^{a_1}p_2^{a_2}…p_k^{a_k}$ 时，利用积性函数的性质可以证明。

    证毕。

莫比乌斯函数

莫比乌斯函数的定义如下：

$\mu(n)= \begin{cases} 1, & n = 1 \ (-1)^m, & n = p_1p_2 \cdots p_m \ 0, & \text{otherwise} \end{cases}$

莫比乌斯函数的性质

  1. 莫比乌斯函数有重要性质：

    $\sum_{d \mid n} \mu(d) = [n == 1]$

    ，其中 $\text{[condition}] = 1$ 当且仅当 $\text{condition}$ 为真。

    证明：

    当 $n = 1$ 时，显然成立。

    当 $n \neq 1$ 时，设 $n = p_1^{a_1}p_2^{a_2} \cdots p_m^{a_m}$，则由定义易得

    $\forall \ d = p_1^{b_1}p_2^{b_2} \cdots p_t^{b_t} \ \text{and}\ d \mid n, \mu(d) = 0$。故

    $\sum_{d \mid n} \mu(d) = \sum_{i=1}^m \binom{m}{i}(-1)^i = (1-1)^m = 0$

    。莫比乌斯函数的该性质使其在莫比乌斯反演中有重要应用。

  2. 莫比乌斯函数是积性函数。证明显然，此略。

其他常见积性函数

  • 约数个数 $d(n) = \sum_{i \mid n} 1$；

  • 约数和 $\sigma(n) = \sum_{i \mid n} i$；

  • 约数的 $m$ 次幂之和 $s(n) = \sum_{i \mid n} i ^m$，可以发现前面两个函数是这个函数的特殊情况；

  • 下面会提到的狄利克雷卷积单位元 $\epsilon(n) = [n == 1]$。

线性筛

利用积性函数的性质，我们可以在线性筛质数的同时筛出积性函数的值。通用过程是考虑将所有数分为三类：

  1. 质数
  2. 最小质因子的指数是 $1$
  3. 最小质因子的指数大于 $1$

据说还有不需要分三种情况而直接按照最小质因子的指数计算的方法，但是我没学过也还不会这里就不讲（逃

欧拉函数

计算 $\varphi(n)$ 的值，按照上面的三类考虑：

  1. $n$ 是质数时，$\varphi(n) = n - 1$；
  2. $n$ 的最小质因子 $q$ 的指数是 $1$ 时，$\varphi(n) = \varphi(m \cdot q) = \varphi(m) \cdot \varphi(q)$；
  3. $n$ 的最小质因子 $q$ 的指数大于 $1$ 时，$\varphi(n) = \varphi(m \cdot q) = q \cdot \varphi(m)$。证明：$\varphi(m \cdot q) = m \cdot q \prod (1- \frac{1}{p_i}) = q \cdot m \cdot \prod (1- \frac{1}{p_i}) = q \cdot \varphi(m)$。

莫比乌斯函数

计算 $\mu(n)$ 的值非常简单，按照上面的三类考虑：

  1. $n$ 是质数时，$\mu(n) = -1$；

  2. $n$ 的最小质因子 $q$ 的指数是 $1$ 时，$\mu(n) = \mu(m \cdot q) = \mu(m) \cdot \mu(q) = -\mu(m)$；

  3. $n$ 的最小质因子 $q$ 的指数大于 $1$ 时，$\mu(n) = 0$。

代码：

bool isPrime[MAXM];
int primes[7300], cnt, phi[MAXM];
void sieve() {
    std::fill(isPrime + 1, isPrime + m + 1, true);
    isPrime[0] = isPrime[1] = false;
    phi[1] = 1;
    mu[1] = 1;
    for (int i = 2; i <= m; ++i) {
        if (isPrime[i]) {
            primes[++cnt] = i;
            phi[i] = i - 1;
            mu[i] = -1;
        }
        for (int j = 1; j <= cnt && i * primes[j] <= m; ++j) {
            isPrime[i * primes[j]] = false;
            if (i % primes[j] == 0) {
                phi[i * primes[j]] = primes[j] * phi[i];
                mu[i * primes[j]] = 0;
                break;
            }
             phi[i * primes[j]] = phi[i] * phi[primes[j]];
             mu[i * primes[j]] = -mu[i];
        }
    }
}

莫比乌斯反演

莫比乌斯反演在 OI 中有许多应用，可以将许多本来难以求解的问题进行简化或加速。

数论分块（整除分块）

参考我的这篇博客。

考虑如何快速求如下函数 $f(n)$ 的值：

$f(n) = \sum_{i = 1}^{n} \left \lfloor \frac{n}{i} \right \rfloor$

写出前几项，可发现 $f(n)$ 的值呈块状分布，对于某个区间内的 $i$，$\left\lfloor\frac{n}{i}\right\rfloor$ 的值相同，这个区间的右端点为 $\frac{n}{n/i}$。这样的区间最多有 $2\sqrt{n}$ 个。整除分块总复杂度 $O(n\sqrt{n})$。

详细证明参见文末参考资料[1]。

狄利克雷卷积

对于两个数论函数 $f, g$，定义狄利克雷卷积运算 $()$ 为 $(fg)(n) = \sum_{d \mid n} f(d)g(\frac{n}{d})$。

卷积运算有如下性质：

  1. 满足交换律：$f*g = g * f$；

  2. 满足结合律：$(f * g) * h = f * (g * h)$；

  3. 存在单位元：

    定义函数

    $\epsilon(n) = \begin{cases} 1, &n = 1 \ 0, &n \neq 1 \end{cases}$

    显然该函数为狄利克雷卷积运算的单位元。

定义函数 $1(n) = 1$。可以发现，由于前文提到的莫比乌斯函数的重要性质，莫比乌斯函数即为函数 $1(n)$ 在狄利克雷卷积下的逆元。即

$(1 * \mu)(n) = \sum_{d \mid n}\mu(d) = [n == 1] = \epsilon(n)$

莫比乌斯反演

莫比乌斯反演的内容如下：

$f(n) = \sum_{d \mid n} g(n) \Leftrightarrow g(n) = \sum_{d \mid n} f(n)\mu(\frac{n}{d})$

证明：

$f = g * 1 \Leftrightarrow \mu * f = g * 1 * \mu \Leftrightarrow g = f * \mu$

从其他方向进行的证明可参见参考资料[2]。

莫比乌斯反演有枚举倍数的第二种形式[3]：

$f(n) = \sum_{n \mid d} g(d) \Rightarrow g(n) = \sum_{n \mid d} f(d)\mu(\frac{d}{n})$

证明：

设 $k = \frac{d}{n}$。若

$f(n) = \sum_{n \mid d} g(d)$

则

$\begin{aligned} \sum_{n \mid d} f(d)\mu(\frac{d}{n}) &= \sum_{k = 1}^{+\infty}\mu(k)f(nk) \ &= \sum_{k = 1}^{+\infty}\mu(k)\sum_{nk \mid t} g(t) \ &= \sum_{n \mid t} g(t) \sum_{k \mid \frac{t}{n}} \mu(k) \end{aligned}$

又由莫比乌斯函数的性质，当且仅当 $\frac{t}{n} = 1$ 即 $t = n$ 时，$\sum_{k \mid \frac{t}{n}} \mu(k) = 1$，否则等于 $0$。

故

$g(n) = \sum_{n \mid t} g(t) \sum_{k \mid \frac{t}{n}} \mu(k) = \sum_{n \mid d} f(d)\mu(\frac{d}{n})$

例题

  1. [HAOI2011]Problem b

    题意：对于给出的 $n$ 个询问，每次求有多少个数对 $(x,y)$，满足 $a \le x \le b$，$c \le y \le d$，且 $\gcd(x,y) = k$ 。

    思路：

    不妨设 $n \leq m$。

    首先我们发现，我们只要能发现一种可以快速求

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

    的方法，就可以根据容斥原理快速求解每个询问。考虑

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

    （考虑每个能对答案作出贡献的 $(i, j)$ 因为都满足 $\gcd(i, j) = k$，所以从这样的 $i$ 和 $j$ 中除去一个 $k$ 后二者一定互质，且其取值一定落在 $\left \lfloor \frac{n}{k} \right \rfloor$ 和 $\left \lfloor \frac{m}{k} \right \rfloor$ 内）

    $= \sum_{i = 1}^{\left\lfloor n / k \right \rfloor} \sum_{j = 1}^{\left\lfloor m / k \right \rfloor} \sum_{d \mid \gcd(i, j)} \mu(d)$

    （上文提到的莫比乌斯函数的性质）

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

    （我们所枚举的 $d$ 一定落在该范围内，考虑对答案有贡献的 $d$ 均满足 $[d \mid \gcd(i, j)]$，因此可以简单地交换求和顺序先枚举 $d$）

    $= \sum_{d = 1}^{\left\lfloor n / k \right \rfloor} \mu(d) \left\lfloor \frac{\left\lfloor n / k \right \rfloor}{d} \right \rfloor \left\lfloor \frac{\left\lfloor m / k \right \rfloor}{d} \right \rfloor$

    （满足 $[d \mid \gcd(i, j)]$ 的 $d$ 一定既是 $i$ 的约数也是 $j$ 的约数，所以能为答案产生 $1$ 的贡献的 $d$ 有 $i$ 与 $j$ 各自的约数个数之积即 $\left\lfloor \frac{\left\lfloor n / k \right \rfloor}{d} \right \rfloor \left\lfloor \frac{\left\lfloor m / k \right \rfloor}{d} \right \rfloor$ 个）

    观察到最后这个式子可以在预处理出 $\mu(n)$ 的前缀和后利用整除分块在 $O(\sqrt{m})$ 的时间内求得。

    代码：

    #include <algorithm>
#include <cstdio>
#include <cstring>
#include <numeric>

const int MAXN = 50000 + 5;

int N, a, b, c, d, k;
bool isNotPrime[MAXN];
int primes[6000], mu[MAXN], cnt;
int summu[MAXN], ans;

void sieve(int m) {
    isNotPrime[0] = isNotPrime[1] = true;
    mu[1] = 1;
    for (int i = 2; i <= m; ++i) {
        if (!isNotPrime[i]) {
            primes[++cnt] = i;
            mu[i] = -1;
        }
        for (int j = 1; j <= cnt && i * primes[j] <= m; ++j) {
            isNotPrime[i * primes[j]] = true;
            if (i % primes[j] == 0) {
                mu[i * primes[j]] = 0;
                break;
            }
            mu[i * primes[j]] = -mu[i];
        }
    }
}

int query(int n, int m) {
    int res = 0;
    for (int l = 1, r; l <= std::min(n, m); l = r + 1) {
        r = std::min(n / (n / l), m / (m / l));
        res += (summu[r] - summu[l - 1]) * (n / l) * (m / l);
    }
    return res;
}

int main() {
    sieve(50000);
    std::partial_sum(mu + 1, mu + 50000 + 1, summu + 1);
    scanf("%d", &N);
    while (N--) {
        scanf("%d %d %d %d %d", &a, &b, &c, &d, &k);
        ans = query(b / k, d / k) - query(b / k, (c - 1) / k) - query((a - 1) / k, d / k) + query((a - 1) / k, (c - 1) / k);
        printf("%d\n", ans);
    }
    return 0;
}

  2. YY 的 GCD

    题意：求

    $\sum_{i = 1}^n \sum_{j = 1}^m [\gcd(i, j) \text{ is prime}]$

    思路：

    设 $\mathbb{P}$ 为质数集。不妨令 $n \le m$。

    $\sum_{i = 1}^n \sum_{j = 1}^m [\gcd(i, j) \in \mathbb{P}] = \sum_{p \in \mathbb{P}}\sum_{i = 1}^n \sum_{j = 1}^m [\gcd(i, j) = p] \ = \sum_{p \in \mathbb{P}} \sum_{d = 1}^{\lfloor n / p \rfloor} \mu(d) \left \lfloor \frac{n}{pd} \right \rfloor \left \lfloor \frac{m}{pd} \right \rfloor$

    直接枚举显然复杂度不对。设 $t = pd$，则

    $\sum_{p \in \mathbb{P}} \sum_{d = 1}^{\lfloor n / p \rfloor} \mu(d) \left \lfloor \frac{n}{pd} \right \rfloor \left \lfloor \frac{m}{pd} \right \rfloor = \sum_{t = 1}^n \sum_{p \mid t, p \in \mathbb{P}}\mu(\frac{t}{p}) \left \lfloor \frac{n}{t} \right \rfloor \left \lfloor \frac{m}{t} \right \rfloor \ = \sum_{t = 1}^n \left \lfloor \frac{n}{t} \right \rfloor \left \lfloor \frac{m}{t} \right \rfloor \sum_{p \mid t, p \in \mathbb{P}}\mu(\frac{t}{p})$

    其中 $\left \lfloor \frac{n}{t} \right \rfloor \left \lfloor \frac{m}{t} \right \rfloor$ 可以整除分块，设 $f(t) = \sum_{p \mid t, p \in \mathbb{P}}\mu(\frac{t}{p})$ ，这个东西可以对于每个质数枚举其倍数 $t$，给 $f(t)$ 的值累加上 $\mu(\frac{t}{k})$ 就可以统计出 $f(t)$ 的值，再求个前缀和就能和整除分块一起计算了。复杂度 $O(能过)$，具体我不会算。

    代码：

    #include <algorithm>
#include <cstdio>
#include <numeric>

const int MAXN = 10000000 + 5;

int T, n, m;
int primes[700000], mu[MAXN], g[MAXN], cnt;
long long sumg[MAXN];
bool isNotPrime[MAXN];

void sieve(int s) {
    isNotPrime[1] = true;
    mu[1] = 1;
    for (int i = 2; i <= s; ++i) {
        if (!isNotPrime[i]) {
            primes[++cnt] = i;
            mu[i] = -1;
        }
        for (int j = 1; j <= cnt && i * primes[j] <= s; ++j) {
            isNotPrime[i * primes[j]] = true;
            if (i % primes[j] == 0) {
                mu[i * primes[j]] = 0;
                break;
            }
            mu[i * primes[j]] = -mu[i];
        }
    }
}

void calcg(int s) {
    for (int i = 1; i <= cnt; ++i) {
        for (int j = 1; primes[i] * j <= s; ++j) {
            g[primes[i] * j] += mu[j];
        }
    }
}

long long solve(int n, int m) {
    long long res = 0;
    for (int l = 1, r, s = std::min(n, m); l <= s; l = r + 1) {
        r = std::min(n / (n / l), m / (m / l));
        res += 1ll * (n / l) * (m / l) * (sumg[r] - sumg[l - 1]);
    }
    return res;
}

int main() {
    sieve(MAXN - 5);
    calcg(MAXN - 5);
    std::partial_sum(g + 1, g + MAXN - 4, sumg + 1);
    scanf("%d", &T);
    while (T--) {
        scanf("%d %d", &n, &m);
        printf("%lld\n", solve(n, m));
    }
    return 0;
}

  3. Crash 的数字表格

    题意：求

    $\sum_{i = 1}^{n} \sum_{j = 1}^{m} \text{lcm}(i, j)$

    思路：

    $\sum_{i = 1}^{n} \sum_{j = 1}^{m} \text{lcm}(i, j) = \sum_{i = 1}^{n} \sum_{j = 1}^{m} \frac{ij}{\gcd(i, j)} \ = \sum_{i = 1}^{n} \sum_{j = 1}^{m} \sum_{d \mid i, d \mid j, \gcd(\frac{i}{d}, \frac{j}{d}) = 1} \frac{ij}{d}$

    （枚举这个最大公因数 $d$，显然同除 $d$ 之后二者互质）

    $= \sum_{d = 1}^n d \sum_{i = 1}^{\left\lfloor n/d \right\rfloor} \sum_{j = 1}^{\left\lfloor m/d \right\rfloor} [\gcd(i, j) == 1] \cdot ij$

    （考虑令原来的 $i = i’ d$，$j$ 同理，然后把 $d$ 提出来，就等价于枚举 $d$ 再在新的上界里枚举互质数对之积之和）

    观察到后面这个东西是互质数对之积之和，我们令

    $s(n, m) = \sum_{i = 1}^{n} \sum_{j = 1}^{m} [\gcd(i, j) == 1] \cdot ij \ = \sum_{i = 1}^{n} \sum_{j = 1}^{m} \sum_{d \mid i, d \mid j} \mu(d) ij \ = \sum_{d = 1}^{n} \sum_{d \mid i}^{n} \sum_{d \mid j}^{m} \mu(d) ij$

    枚举的 $d$ 必是 $i$ 和 $j$ 的因数，所以我们改换顺序先枚举 $d$ 再枚举其倍数。再考虑令原来的 $i = i’ d$，$j$ 同理，然后把 $d$ 和 $\mu(d)$ 提出来：

    $= \sum_{d = 1}^{n} d^2 \mu(d) \sum_{i = 1}^{\left\lfloor n/d \right\rfloor} \sum_{j = 1}^{\left\lfloor n/d \right\rfloor} ij$

    根据简单的数学知识可以看出来上面这个式子的后面那两个和式的结果可以 $O(1)$ 地算出来，于是 $s(n, m)$ 的值可以用整除分块快速算出。然后我们回到一开始的式子：

    $\sum_{d = 1}^n d \sum_{i = 1}^{\left\lfloor n/d \right\rfloor} \sum_{j = 1}^{\left\lfloor m/d \right\rfloor} [\gcd(i, j) == 1] \cdot ij \ = \sum_{d = 1}^n d \cdot s(\left\lfloor n/d \right\rfloor, \left\lfloor m/d \right\rfloor)$

    我们发现如果 $s(n, m)$ 易求的话这个式子又可以整除分块。所以应用两次整除分块处理求解，整体复杂度呈线性。

    代码（取模部分很丑）：

    #include <algorithm>
#include <cstdio>
#include <numeric>

#define POSI(x) (((x) + MOD) % MOD)

const int MAXN = 10000000 + 5, MOD = 20101009;

int n, m;
int primes[700000], mu[MAXN], mud[MAXN], summud[MAXN], sum[MAXN], cnt;
bool isNotPrime[MAXN];

void sieve(int s) {
    isNotPrime[1] = true;
    mud[1] = mu[1] = 1;
    for (int i = 2; i <= s; ++i) {
        if (!isNotPrime[i]) {
            primes[++cnt] = i;
            mu[i] = -1;
            mud[i] = POSI(1ll * i * i % MOD * mu[i]);
        }
        for (int j = 1; j <= cnt && i * primes[j] <= s; ++j) {
            isNotPrime[i * primes[j]] = true;
            if (i % primes[j] == 0) {
                mud[i * primes[j]] = mu[i * primes[j]] = 0;
                break;
            }
            mu[i * primes[j]] = -mu[i];
            mud[i * primes[j]] = POSI(1ll * i * primes[j] % MOD * i % MOD *
                                 primes[j] % MOD * mu[i * primes[j]] % MOD);
        }
    }
}

inline int t(int a, int b) {
    return 1ll * (1ll * a * (a + 1) / 2 % MOD) * (1ll * b * (b + 1) / 2 % MOD) % MOD;
}

int s(int a, int b) {
    int res = 0;
    for (int l = 1, r; l <= std::min(a, b); l = r + 1) {
        r = std::min(a / (a / l), b / (b / l));
        res = (1ll * (summud[r] - summud[l - 1]) % MOD * t(a / l, b / l) % MOD + res) % MOD;
    }
    return res;
}

int f(int a, int b) {
    int res = 0;
    for (int l = 1, r; l <= std::min(a, b); l = r + 1) {
        r = std::min(a / (a / l), b / (b / l));
        res = (1ll * (sum[r] - sum[l - 1]) * s(a / l, b / l) % MOD + res) % MOD;
    }
    return res;
}

int main() {
    scanf("%d %d", &n, &m);
    if (n > m) std::swap(n, m);
    sieve(m);
    for (int i = 1; i <= m; ++i) {
        summud[i] = (1ll * summud[i - 1] + mud[i]) % MOD;
        sum[i] = (1ll * sum[i - 1] + i) % MOD;
    }
    printf("%d", POSI(f(n, m)));
    return 0;
}

参考资料