[Luogu P1829] [BZOJ 2154] [国家集训队]Crash的数字表格

洛谷传送门

BZOJ传送门

题目背景

提示：原 P1829 半数集问题 已经迁移至 P1028 数的计算

题目描述

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

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

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$的值。

输入输出格式

输入格式：

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

输出格式：

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

输入输出样例

输入样例#1：

4 5

输出样例#1：

122

说明

30%的数据满足$N, M≤ 10^3$。

70%的数据满足$N, M≤ 10^5$。

100%的数据满足$N, M≤ 10^7$。

解题分析

首先推一波式子，这里我们不妨设$n<=m$：

$$\sum_{i=1}^{n}\sum_{j=1}^{m}lcm(i,j) \\= \sum_{i=1}^{n}\sum_{j=1}^{m}\frac{ij}{gcd(i,j)} \\= \sum_{d=1}^{n}\sum_{i=1}^{n}\sum_{j=1}^{m}\frac{ij}{d}[gcd(i,j)=d] \\= \sum_{d=1}^{n}d\sum_{i=1}^{\lfloor \frac{n}{d}\rfloor}\sum_{j=1}^{\lfloor\frac{m}{d}\rfloor}ij[gcd(i,j)=1]$$
现在我们单独拿出右边的一堆进行处理。设 $F(x,y)=\sum_{i=1}^x\sum_{j=1}^yij[gcd(i,j)=1](x<=y)$，那么
$$F(x,y)=\sum_{i=1}^{x}\sum_{j=1}^{y}ij\sum_{k|gcd(i,j)}\mu(k) \\ =\sum_{k=1}^x\mu(k)\sum_{i=1}^x\sum_{j=1}^yij[k|gcd(i,j)]\ \\ =\sum_{k=1}^x\mu(k)\times k^2\sum_{i=1}^{\lfloor\frac{x}{k}\rfloor}\sum_{j=1}^{\lfloor \frac{y}{k}\rfloor}[1|gcd(i,j)] \\ =\sum_{k=1}^x\mu(k)\times k^2\times \frac{(1+\lfloor\frac{x}{k}\rfloor)\times\lfloor\frac{x}{k}\rfloor}{2}\times \frac{(1+\lfloor\frac{y}{k}\rfloor)\times\lfloor\frac{y}{k}\rfloor}{2}$$
将这个式子代回原式，得到：
$$\sum_{i=1}^{n}\sum_{j=1}^{m}lcm(i,j) \\= \sum_{d=1}^{n}d\sum_{k=1}^{\lfloor \frac{n}{d}\rfloor}\mu(k)\times k^2\times \frac{(1+\lfloor\frac{n}{kd}\rfloor)\times\lfloor\frac{n}{kd}\rfloor}{2}\times \frac{(1+\lfloor\frac{m}{kd}\rfloor)\times\lfloor\frac{m}{kd}\rfloor}{2}$$
看起来太凌乱， 而且已经没法化简了， 于是我们套路地更改枚举项，设 $T=kd$：
$$\sum_{d=1}^{n}d\sum_{k=1}^{\lfloor \frac{n}{d}\rfloor}\mu(k)\times k^2\times \frac{(1+\lfloor\frac{n}{kd}\rfloor)\times\lfloor\frac{n}{kd}\rfloor}{2}\times \frac{(1+\lfloor\frac{m}{kd}\rfloor)\times\lfloor\frac{m}{kd}\rfloor}{2} \\=\sum_{T=1}^{n}\frac{(1+\lfloor\frac{n}{T}\rfloor)\times\lfloor\frac{n}{T}\rfloor}{2}\times \frac{(1+\lfloor\frac{m}{T}\rfloor)\times\lfloor\frac{m}{T}\rfloor}{2} T\sum_{k|T}k\times \mu(k)$$
我们要做的最后一件事就是快速求出 $\sum_{k|T}k\times \mu(k)$了（因为前面可以下底分块快速求）。

考虑$f(x)=\sum_{k|x}k\times \mu(k)$。

• 当$x$为质数时， 显然$f(x)=1\times 1+x\times (-1)=1-x$。
• 当$x=pri\times m$时且$pri\nmid m$时， 除了原来的因数， 我们还引入了$pri\times a_1,pri\times a_2...(a_1,a_2...|m)$这些因数， 而$\mu(pri\times a_1)=-\mu(a_1)$， 所以$f(x)=-pri\times f(m)+f(m)$。
• 当$x=pri\times m$时且$pri|m$时， 没有新的约数， 而$\mu(pri^2)=\mu(pri^2\times a_1)...=0$， 所以$f(pri\times m)=f(m)$。

最后这个式子与$\mu$一点关系都没有， 筛都不需要筛出来…

我们对$x\times f(x)$做一个前缀和， 套上下底分块就可以$AC$了。

总复杂度$O(N+\sqrt N)$。

代码如下：

#include <cstdio>
#include <cstring>
#include <cstdlib>
#include <algorithm>
#include <cctype>
#include <cmath>
#define R register
#define IN inline
#define gc getchar()
#define W while
#define ll long long
#define MOD 20101009ll
#define MX 10000500
int pri[MX], pcnt, n, m;
ll F[MX], sum[MX];
bool npr[MX];
void get()
{
    F[1] = sum[1] = 1;
    R int i, j, tar;
    for (i = 2; i <= n; ++i)
    {
        if(!npr[i]) pri[++pcnt] = i, F[i] = 1 - i;
        for (j = 1; j <= pcnt; ++j)
        {
            tar = pri[j] * i;
            if(tar > n) break;
            npr[tar] = true;
            if(!(i % pri[j])) {F[tar] = F[i]; break;}
            F[tar] = F[i] * (1 - pri[j]);
        }
        sum[i] = (sum[i - 1] + F[i] * i % MOD + MOD) % MOD;
    }
}
int main(void)
{
    scanf("%d%d", &n, &m);
    if(n > m) std::swap(n, m);
    get(); ll ans = 0;
    for (R int lef = 1, rig, sum1, sum2; lef <= n; lef = rig + 1)
    {
        rig = std::min(n / (n / lef), m / (m / lef));
        sum1 = 1ll * (n / lef) * (n / lef + 1) / 2 % MOD;
        sum2 = 1ll * (m / lef) * (m / lef + 1) / 2 % MOD;
        ans  += 1ll * sum1 * sum2 % MOD * (sum[rig] - sum[lef - 1] + MOD) % MOD;
        ans %= MOD;
    }
    printf("%lld", (ans % MOD + MOD) % MOD);
}