Crash的数字表格
【题目描述】
今天的数学课上，Crash小朋友学习了最小公倍数(Least Common Multiple)。对于两个正整数\(a\)和\(b\)，\(\operatorname{lcm}(a,b)\)表示能同时被\(a\)和\(b\)整除的最小正整数。例如，\(\operatorname{lcm}(6,8)\) = \(24\)。回到家后，Crash还在想着课上学的东西，为了研究最小公倍数，他画了一张\(N\times M\)的表格。每个格子里写了一个数字，其中第i行第j列的那个格子里写着数为\(\operatorname{lcm}(i,j)\)。一个\(4\times 5\)的表格如下：
\(\begin{matrix}1&2&3&4&5\\2&2&6&4&10\\3&6&3&12&15\\4&4&12&4&20\\ \end{matrix}\)
看着这个表格，Crash想到了很多可以思考的问题。不过他最想解决的问题却是一个十分简单的问题：这个表格中所有数的和是多少。当N和M很大时，Crash就束手无策了，因此他找到了聪明的你用程序帮他解决这个问题。由于最终结果可能会很大，Crash只想知道表格里所有数的和 \(\operatorname{mod} 20101009\) 的值。
【输入格式】
输入的第一行包含两个正整数，分别表示\(N\)和\(M\)。
【输出格式】
输出一个正整数，表示表格中所有数的和 \(\operatorname{mod} 20101009\) 的值。
\(n \le 10^7\)
万幸20101009是个质数
题解
不妨设$n\le m$那就直接开始推式子吧。。。
懵逼警告
需要求的是 \(\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\operatorname{lcm}(i,j)\)
\(=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\frac{i*j}{\gcd(i,j)}\)
枚举\(\gcd(i,j)\)的值
\(=\sum\limits_{g=1}^{n}\frac{1}{g}*\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}[\gcd(i,j)=g]*i*j\)
把后半部分提出来一个\(g^2\)
\(=\sum\limits_{g=1}^{n}g*\sum\limits_{i=1}^{\lfloor\frac{n}{g}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{g}\rfloor}[\gcd(i,j)=1]*i*j\)
套用公式\([n=1]=\sum\limits_{d|n}\mu(d)\)
\(=\sum\limits_{g=1}^{n}g*\sum\limits_{i=1}^{\lfloor\frac{n}{g}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{g}\rfloor}i*j*\sum\limits_{d|\gcd(i,j)}\mu(d)\)
如法炮制 枚举\(d\)的值
\(=\sum\limits_{g=1}^{n}g*\sum\limits_{d=1}^{\lfloor\frac{n}{g}\rfloor}\mu(d)*\sum\limits_{i=1}^{\lfloor\frac{n}{g*d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{g*d}\rfloor}d^2*i*j\)
提出\(d^2\)
\(=\sum\limits_{g=1}^{n}g*\sum\limits_{d=1}^{\lfloor\frac{n}{g}\rfloor}\mu(d)*d^2*\sum\limits_{i=1}^{\lfloor\frac{n}{g*d}\rfloor}\sum\limits_{j=1}^{\lfloor\frac{m}{g*d}\rfloor}i*j\)
后面的两个\(\sum\)可以用等差数列求和的方式表示
\(=\sum\limits_{g=1}^{n}g*\sum\limits_{d=1}^{\lfloor\frac{n}{g}\rfloor}\mu(d)*d^2*\frac{\lfloor\frac{n}{g*d}\rfloor*(\lfloor\frac{n}{g*d}\rfloor+1)}{2}*\frac{\lfloor\frac{m}{g*d}\rfloor*(\lfloor\frac{m}{g*d}\rfloor+1)}{2}\)
好了 这个式子就是最终的形式了
观察一下 枚举\(g\) 我们只需要知道 \(\frac{n}{g}\) 及 \(\frac{m}{g}\) 是多少， \(\sum\limits_{d=1}^{\lfloor\frac{n}{g}\rfloor}\mu(d)*d^2*\frac{\lfloor\frac{n}{g*d}\rfloor*(\lfloor\frac{n}{g*d}\rfloor+1)}{2}*\frac{\lfloor\frac{m}{g*d}\rfloor*(\lfloor\frac{m}{g*d}\rfloor+1)}{2}\) 这一串就能算出来
所以用数论分块来枚举\(g\)就可以了 记\(\lfloor\frac{n}{g}\rfloor =n_0,\lfloor\frac{m}{g}\rfloor =m_0\)
同理 只需要知道\(\frac{n_0}{d}\)和\(\frac{m_0}{d}\)， 里面那串 \(\frac{\lfloor\frac{n}{g*d}\rfloor*(\lfloor\frac{n}{g*d}\rfloor+1)}{2}*\frac{\lfloor\frac{m}{g*d}\rfloor*(\lfloor\frac{m}{g*d}\rfloor+1)}{2}\) 也就能算出来了
所以\(d\)也可以用数论分块枚举 处理一下 \(\mu(d)*d^2\) 这个东西的前缀和就可以了
时间复杂度是\(O(\sqrt{n}^2)=O(n)\)，很完美
注意处理负数 \(\mu(d)*d^2\)的前缀和会有负数
推导很冗长 代码很简洁
#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
inline ll read() {
ll x = 0, f = 1; char ch = getchar();
for (; ch > '9' || ch < '0'; ch = getchar()) if (ch == '-') f = -1;
for (; ch <= '9' && ch >= '0'; ch = getchar()) x = (x << 1) + (x << 3) + (ch ^ '0');
return x * f;
}
const ll mod = 20101009;
const ll div2 = mod / 2 + 1;
ll n, m, ans;
bool np[10000005];
ll prime[1000005], mob[10000005], sum[10000005], tot;
inline void getprime() {
mob[1] = 1;
for (ll i = 2; i <= max(n, m); i++) {
if (!np[i]) prime[++tot] = i, mob[i] = -1;
for (ll j = 1; j <= tot && i * prime[j] <= max(n, m); j++) {
np[i * prime[j]] = 1;
if (i % prime[j] == 0) {
mob[i * prime[j]] = 0;
break;
} else mob[i * prime[j]] = -mob[i];
}
}
for (ll i = 1; i <= max(n, m); i++) {
sum[i] = (sum[i-1] + i * i % mod * mob[i] % mod) % mod;
}
}
inline ll calc(ll n0, ll m0) {
ll ret = 0;
for (ll l = 1, r; l <= min(n0, m0); l = r + 1) {
r = min(n0 / (n0 / l), m0 / (m0 / l));
ll s = (sum[r] - sum[l-1] + mod) % mod;
ll s1 = (n0 / l) * (n0 / l + 1) % mod * div2 % mod, s2 = (m0 / l) * (m0 / l + 1) % mod * div2 % mod;
ret = ((ret + s * s1 % mod * s2 % mod) % mod + mod) % mod;
}
return ret;
}
int main() {
n = read(); m = read();
getprime();
for (ll l = 1, r; l <= min(n, m); l = r + 1) {
r = min(n / (n / l), m / (m / l));
ans = (((r - l + 1) * (r + l) % mod * div2 % mod * calc(n / l, m / l) % mod + ans) % mod + mod) % mod;
}
printf("%lld\n", ans);
return 0;
}
