[题解] Crash的数字表格 / JZPTAB

Table of Contents

题目描述

今天的数学课上，Crash 小朋友学习了最小公倍数 (Least Common Multiple)。对于两个正整数 $a$ 和 $b$ , $l c m (a, b)$ 表示能同时整除 $a$ 和 $b$ 的最小正整数。例如， $l c m (6, 8) = 24$ 。

回到家后，Crash 还在想着课上学的东西，为了研究最小公倍数，他画了一张 $n \times m$ 的表格。每个格子里写了一个数字，其中第 $i$ 行第 $j$ 列的那个格子里写着数为 $l c m (i, j)$ 。一个 $4 \times 5$ 的表格如下：

	1	2	3	4	5
1	1	2	3	4	5
2	2	2	6	4	10
3	3	6	3	12	15
4	4	4	12	4	20

看着这个表格，Crash想到了很多可以思考的问题。不过他最想解决的问题却是一个十分简单的问题：这个表格中所有数的和是多少。当 $n$ 和 $m$ 很大时，Crash 就束手无策了，因此他找到了聪明的你用程序帮他解决这个问题。由于最终结果可能会很大，Crash 只想知道表格里所有数的和 $mod 20101009$ 的值。

输入格式

输入的第一行包含两个正整数，分别表示 $n$ 和 $m$

输出格式

输出一个正整数，表示表格中所有数的和 $mod 20101009$ 的值

样例输入

4 5

样例输出

数据范围及提示

$30 %$ 的数据满足 $n, m \leq 10^{3}$

$70 %$ 的数据满足 $n, m \leq 10^{5}$

$100 %$ 的数据满足 $n, m \leq 10^{7}$

解题思路

$\begin{aligned} \sum_{i = 1}^{n} \sum_{j = 1}^{m} \frac{i j}{gcd (i, j)} \\ = & \sum_{i = 1}^{n} \sum_{j = 1}^{m} \sum_{d = 1}^{min (n, m)} \frac{i j}{d} [g c d (i, j) = d] \\ = & \sum_{i = 1}^{n} \sum_{j = 1}^{m} \sum_{d = 1}^{min (n, m)} \frac{i j}{d} [g c d (i, j) = d] \\ = & \sum_{d = 1}^{min (n, m)} d^{2} \sum_{i = 1}^{⌊ \frac{n}{d} ⌋} \sum_{j = 1}^{⌊ \frac{m}{d} ⌋} \frac{i j}{d} [g c d (i, j) = 1] \\ = & \sum_{d = 1}^{min (n, m)} d \sum_{i = 1}^{⌊ \frac{n}{d} ⌋} \sum_{j = 1}^{⌊ \frac{m}{d} ⌋} [g c d (i, j) = 1] i j \\ = & \sum_{d = 1}^{min (n, m)} d \sum_{i = 1}^{⌊ \frac{n}{d} ⌋} \sum_{j = 1}^{⌊ \frac{m}{d} ⌋} \sum_{k ∣ gcd (i, j)} μ (k) i j \\ = & \sum_{d = 1}^{min (n, m)} d \sum_{i = 1}^{⌊ \frac{n}{d} ⌋} \sum_{j = 1}^{⌊ \frac{m}{d} ⌋} \sum_{k ∣ i} \sum_{k ∣ j} μ (k) i j \end{aligned}$

这个式子可以被分成两部分解决，我们令
$f (n, m) = \sum_{i = 1}^{n} \sum_{j = 1}^{m} \sum_{k ∣ i} \sum_{k ∣ j} μ (k) i j$
那么原式就等于
$\sum_{d = 1}^{min (n, m)} d \cdot f (⌊ \frac{n}{d} ⌋, ⌊ \frac{m}{d} ⌋)$
可以通过整除分块计算，而 $f (n, m)$ 的求解则是经典问题
$\begin{aligned} f (n, m) & = \sum_{i = 1}^{n} \sum_{j = 1}^{m} \sum_{k ∣ i} \sum_{k ∣ j} μ (k) i j \\ = \sum_{k = 1}^{min (n, m)} μ (k) \sum_{i = 1}^{n} \sum_{i = 1}^{m} i j [k ∣ i] [k ∣ j] \\ = \sum_{k = 1}^{min (n, m)} μ (k) k^{2} \sum_{i = 1}^{⌊ \frac{n}{k} ⌋} \sum_{j = 1}^{⌊ \frac{m}{k} ⌋} i j \end{aligned}$
令
$\begin{aligned} g (n, m) & = \sum_{i = 1}^{n} \sum_{j = 1}^{m} i j \\ = \frac{(n + 1) \times n}{2} \times \frac{(m + 1) \times m}{2} \end{aligned}$
那么式子可以写作
$f (n, m) = \sum_{k = 1}^{min (n, m)} μ (k) k^{2} g (⌊ \frac{n}{k} ⌋, ⌊ \frac{m}{k} ⌋)$
可以整除分块计算

代码

#include <iostream>
#include <cstdio>
#include <cstring>
#include <cmath>
#include <algorithm>

using namespace std;
const int MAXN = 10000000 + 10;
const int MOD = 20101009;

inline int read() {
    int x = 0;
    char ch = getchar();

    while (ch < '0' || ch > '9') {
        ch = getchar();
    }

    while ('0' <= ch && ch <= '9') {
        x = x * 10 + ch - '0';
        ch = getchar();
    }

    return x;
}

int prime[MAXN], mu[MAXN], cnt;
bool isNotPrime[MAXN];

long long sum[MAXN];

inline void init (int n) {
    mu[1] = 1;

    for (register int i=2; i<=n; i++) {
        if (!isNotPrime[i]) {
            prime[++cnt] = i;
            mu[i] = -1;
        }

        for (register int j=1; j<=cnt; j++) {
            int v = i * prime[j];
            if (v > n) break;

            isNotPrime[v] = true;

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

    for (register int i=1; i<=n; i++) {
        sum[i] = (sum[i-1] + ((1LL * i * i) % MOD * (mu[i] + MOD) % MOD)) % MOD;
    }
}

inline long long calc2(int n, int m) {
    return (1LL * n * (n + 1) / 2 % MOD) * (1LL * m * (m + 1) / 2 % MOD) % MOD;
}

inline long long calc(int n, int m) {
    long long ans = 0;
    int r = 0;

    for (register int l=1; l<=min(n, m); l = r + 1) {
        r = min(n / (n / l), m / (m / l));
        ans = (ans + (sum[r] - sum[l-1] + MOD) % MOD * calc2(n / l, m / l) % MOD) % MOD;
    }

    return ans;
}

inline long long solve(int n, int m) {
    long long ans = 0;
    int r = 0;

    for (register int l=1; l<=min(n, m); l = r + 1) {
        r = min(n/(n/l), m/(m/l));
        ans = (ans + ((1LL * (l + r) * (r - l + 1) / 2)) % MOD * calc(n / l, m / l)) % MOD;
    }

    return ans;
}

int main(){
    int n = read();
    int m = read();

    init(max(n, m));

    printf("%lld\n", solve(n, m));
    return 0;
}