BZOJ 2154: Crash的数字表格

Description

今天的数学课上，Crash小朋友学习了最小公倍数(Least Common Multiple)。对于两个正整数a和b，LCM(a, b)表示能同时被a和b整除的最小正整数。例如，LCM(6, 8) = 24。回到家后，Crash还在想着课上学的东西，为了研究最小公倍数，他画了一张NM的表格。每个格子里写了一个数字，其中第i行第j列的那个格子里写着数为LCM(i, j)。一个45的表格如下： 1 2 3 4 5 2 2 6 4 10 3 6 3 12 15 4 4 12 4 20 看着这个表格，Crash想到了很多可以思考的问题。不过他最想解决的问题却是一个十分简单的问题：这个表格中所有数的和是多少。当N和M很大时，Crash就束手无策了，因此他找到了聪明的你用程序帮他解决这个问题。由于最终结果可能会很大，Crash只想知道表格里所有数的和mod 20101009的值。

Input

输入的第一行包含两个正整数，分别表示N和M。

Output

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

Sample Input

1
4 5

Sample Output

1
122

HINT

100%的数据满足N, M ≤ 107。

Solution

先祭出宝典：莫比乌斯反演.pdf （by PoPoQQQ）

貌似pdf里没说怎么推出来的，其实不难：

\begin{eqnarray}

& & \sum _{1 \le i \le x} \sum _{1 \le j \le y} [(i,j)=1]*ij \

& = & \sum _{1 \le i \le x} \sum _{1 \le j \le y} \sum _{d | (i,j)} \mu (d) * ij \

& = & \sum _{1 \le d \le \min (x, y)} \mu (d) \sum _{1 \le i \le x , d|i } \sum _{1 \le j \le y , d|j} ij \

& = & \sum _{1 \le d \le \min (x, y)} \mu (d) d^2 \sum _{1 \le i \le \lfloor \frac{x}{d} \rfloor } \sum _{1 \le j \le \lfloor \frac{y}{d} \rfloor} ij \

& = & \sum _{1 \le d \le \min (x, y)} \mu (d) d^2 * {Sum}( \lfloor \frac{x}{d} \rfloor, \lfloor \frac{y}{d} \rfloor)

\end{eqnarray}

然后对于ans和F都\( \Theta ( \sqrt n) \)求，得到\( \Theta (n) \) 的总复杂度。

Code

 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
#include <cstdio>
#include <algorithm>
using namespace std;
typedef long long LL;
const LL mod = 20101009LL;
const int MAXN = 10000009;
LL N, M, pre[MAXN], Ans;
int mu[MAXN], cnt = 0, P[MAXN], vis[MAXN];
inline LL Sum(LL x, LL y) { return x*(x+1LL)%mod*y%mod*(y+1LL)%mod*15075757LL%mod; }
void Linear_Shaker(int N) {
    mu[1] = 1;
    for (int i = 2; i <= N; ++i) {
        if (!vis[i]) {
            P[++cnt] = i;
            mu[i] = -1;
        }
        for (int j = 1; j <= cnt && (LL)P[j] * i <= LL(N); ++j) {
            vis[i * P[j]] = 1;
            if (i % P[j] == 0) {
                mu[i * P[j]] = 0;
                break;
            } else mu[i * P[j]] = -mu[i];
        }
    }
    for (int i = 1; i <= N; ++i)
        pre[i] = (pre[i-1] + ((LL)mu[i] * i * i)) % mod;
}
LL F(LL x, LL y) {
    LL ret = 0LL; if (x > y) swap(x, y);
    for (LL i = 1LL, last; i <= x; i = last + 1LL) {
        last = min(x / (x / i), y / (y / i));
        ret += (pre[last] - pre[i - 1]) % mod * Sum(x / i, y / i) % mod;
        ret %= mod;
    }
    return ret;
}
LL Calc(LL N, LL M) {
    LL ret = 0LL;
    for (LL i = 1LL, last; i <= N; i = last + 1LL) {
        last = min(N / (N / i), M / (M / i));
        ret += ((i + last) * (last - i + 1LL) % mod * 10050505LL) % mod * F(N / i, M / i) % mod;
        ret %= mod;
    }
    return ret;
}
int main() {
    scanf("%lld%lld", &N, &M);
    if (N > M) swap(N, M);
    Linear_Shaker(M);
    Ans = (Calc(N, M) + mod) % mod;
    printf("%lld", Ans);
}