题解 P1829 [国家集训队]Crash的数字表格 / JZPTAB

题目背景

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

题目描述

今天的数学课上，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只想知道表格里所有数的和mod20101009的值。

输入输出格式

输入格式：

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

输出格式：

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

思路讲解

这是一道莫比乌斯反演和数论分块的练手题。

推一波公式：
$$\text{ans}=\sum _{i=1}^n\sum _{j=1}^m[i,j]=\sum _{i=1}^n\sum _{j=1}^m\frac{i\cdot j}{(i,j)}$$
令$d=(i,j)$，则必有：$(\frac{i}{d},\frac{j}{d})=1$，于是：
$$=\sum _{i=1}^n\sum _{j=1}^m\sum _{d|i\wedge d|j\wedge (\frac{i}{d},\frac{j}{d})=1}\frac{i\cdot j}{d}$$
考虑用$d$枚举$i,j$，则我们可以把和式提前：
$$=\sum _{d=1}^{\min (n,m)}d\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{d}\rfloor }[(i,j)=1]\cdot i\cdot j$$
我们令：
$$\text{Sum}(n,m)=\sum _{i=1}^n\sum _{j=1}^m[(i,j=1)]\cdot i\cdot j$$
则：
$$\text{ans}=\sum _{d=1}^{\min (n,m)}d\cdot \text{Sum}(\lfloor \frac{n}{d}\rfloor ,\lfloor \frac{m}{d}\rfloor )$$
考虑求$\text{Sum}(n,m)$:

根据套路，套莫比乌斯反演:
$$=\sum _{d=1}^{\min (n,m)}\sum _{d|i}^n\sum _{d|j}^m\mu (d)\cdot i\cdot j$$
可以将$\mu (d)$提前：
$$=\sum _{d=1}^{\min (n,m)}\mu (d)d^2\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{d}\rfloor }i\cdot j$$
由拉格朗日插值公式，我们可以知道：
$$L(n,m)=\sum _{i=1}^n\sum _{j=1}^{m}ij=\sum _{i=1}^{n}i\sum _{j=1}^{m}j=\frac{n(n+1)}{2}\times \frac{m(m+1)}{2}$$
于是：
$$\text{Sum}(n,m)=\sum _{d=1}^{\min (n,m)}\mu (d)d^{2}L(\lfloor \frac{n}{d}\rfloor ,\lfloor \frac{m}{d}\rfloor )$$
于是可以数论分块求。

定一个前缀和，令：
$$S_d=\sum _{i=1}^d\mu (i)i^2$$
就可以了。

时间复杂度：

我估计是$\mathcal{O}(\min (n,m))$的吧！有一定的常数。

#include <cstdio>
#include <iostream>

using namespace std;
const int N = 1E+7 + 1;
const int M = 20101009;
int S[N], prime[N/10], Mu[N], n, m, tot;
bool flag[N];

void sieve() {
	Mu[1] = 1;
	for(int i = 2; i <= min(n, m); i++) {
		if(!flag[i]) prime[++tot] = i, Mu[i] = -1;
		for(int j = 1; j <= tot && i * prime[j] <= min(n, m); j++) {
			flag[i * prime[j]] = 1;
			if(i % prime[j] == 0) { Mu[i * prime[j]] = 0; break; }
			Mu[i * prime[j]] = -Mu[i];
		}
	}
	for(int d = 1; d <= min(n, m); d++) 
		S[d] = (S[d - 1] + 1LL * d * d % M * (Mu[d] + M)) % M;
}
int G(int n, int m) {
	return (1LL * n * (n + 1) / 2 % M) * (1LL * m * (m + 1) / 2 % M) % M;
}
int Sum(int n, int m) {
	int res = 0;
	for(int l = 1, r; l <= min(n, m); l = r + 1) {
		r = min(n / (n / l), m / (m / l));
		res = (res + 1LL * (S[r] - S[l - 1] + M) * G(n / l, m / l) % M) % M;
	}
	return res;
}
int cal(int n, int m) {
	int res = 0;
	for(int l = 1, r; l <= min(n, m); l = r + 1) {
		r = min(n / (n / l), m / (m / l));
		res = (res + 1LL * (l + r) * (r - l + 1) / 2 % M * Sum(n / l, m / l) % M) % M;
	}
	return res;
}
int main()
{
	scanf("%d %d", &n, &m);
	sieve();
	printf("%d", cal(n, m));
	return 0;
}