[国家集训队]Crash的数字表格
知识点: 莫比乌斯反演
原题面:Luogu
扯
取模,很有精神!
题意简述
给定 \(n,m\) , 求:
\[\sum_{i=1}^n\sum_{j=1}^{m} \operatorname{lcm}(i,j)\bmod 20101009 \]\(1\le n,m\le 10^7\)。
分析题意
易知原式等价于:
考虑枚举 \(\gcd(i,j)\),设枚举量为 \(d\)。
则 \(d=\gcd(i,j)\) 的充要条件是满足 \(d|i, d|j\) 且 \(\gcd(\dfrac{i}{d},\dfrac{j}{d}) = 1\),则原式等价于:
先枚举 \(d\),则原式等价于:
这个 \(d\) 很烦人,把 \(i,j\) 中的 \(d\) 提出来,变为枚举 \(\frac{i}{d}\),\(\frac{j}{d}\)。
消去 \(d\mid i\),\(d\mid j\) 的限定条件,则原式等价于:
单独考虑后半部分,设 \(f(x,y) = \sum\limits_{i=1}^{x} \sum\limits_{j=1}^{y}[\gcd(i,j)=1]ij\)。
发现 \(f(x,y)\) 的形式与 [HAOI2011] Problem b 中的式子类似,代入 \([\gcd(i,j) = 1] = \sum\limits_{d\mid \gcd(i,j)} {\mu (d)}\) 套路一波:
前半部分 \(\sum\limits_{d=1}\mu(d) d^2\),可以考虑筛出 \(\mu(d)\) 后求前缀和。
后半部分是等差数列乘等差数列的形式,设 \(g(p,q) = \sum\limits_{i=1}^{p} \sum\limits_{j=1}^{q}ij\),\(g_{p,q}\) 可以通过下式 \(O(1)\) 计算:
则对于 \(f(x,y)\),有:
数论分块求解即可。
再看回原式,原式等价于:
又是一个可以数论分块求解的形式。
线性筛预处理后 数论分块套数论分块,复杂度 \(O(n + m)\),瓶颈是线性筛。
爆零小技巧
处理时会出现求平方的运算,需要特别注意取模问题,ll 都会爆,被坑惨了。
在预处理前缀和的这个地方:
sum[i] = (sum[i - 1] + 1ll * i * i % kMod * (mu[i] + kMod)) % kMod; //仅令 mu + kMod
注意先对平方取模,在乘上 \(\mu\),否则就会爆掉。
以及可以仅令 \(mu + mod\)。
以及这个地方:
int g(int n_, int m_) {
return (1ll * n_ * (n_ + 1ll) / 2ll % kMod) * (1ll * m_ * (m_ + 1ll) / 2ll % kMod) % kMod;
}
平方计算,注意随时取模。
代码实现
//莫比乌斯反演
/*
By:Luckyblock
*/
#include <algorithm>
#include <cctype>
#include <cmath>
#include <cstdio>
#include <cstring>
#define ll long long
const ll kMod = 20101009;
const int kMaxn = 1e7 + 10;
//=============================================================
int pnum, p[kMaxn];
ll mu[kMaxn], sum[kMaxn];
bool vis[kMaxn];
//=============================================================
inline int read() {
int f = 1, w = 0;
char ch = getchar();
for (; !isdigit(ch); ch = getchar())
if (ch == '-') f = -1;
for (; isdigit(ch); ch = getchar()) w = (w << 3) + (w << 1) + (ch ^ '0');
return f * w;
}
void Getmax(int &fir_, int sec_) {
if (sec_ > fir_) fir_ = sec_;
}
void Getmin(int &fir_, int sec_) {
if (sec_ < fir_) fir_ = sec_;
}
void Euler(int lim_) {
vis[1] = true, mu[1] = 1ll;
for (int i = 2; i <= lim_; ++ i) {
if (! vis[i]) {
p[++ pnum] = i;
mu[i] = - 1;
}
for (int j = 1; j <= pnum && i * p[j] <= lim_; ++ j) {
vis[i * p[j]] = true;
if (i % p[j] == 0) { //勿忘平方因子
mu[i * p[j]] = 0;
break;
}
mu[i * p[j]] = - mu[i];
}
}
sum[1] = 1ll;
for (int i = 1; i <= lim_; ++ i) {
sum[i] = (sum[i - 1] + 1ll * i * i % kMod * (mu[i] + kMod)) % kMod; //仅令 mu + kMod
}
}
int g(int n_, int m_) {
return (1ll * n_ * (n_ + 1ll) / 2ll % kMod) * (1ll * m_ * (m_ + 1ll) / 2ll % kMod) % kMod;
}
int f(int n_, int m_) {
int lim = std :: min(n_, m_), ret = 0;
for (int l = 1, r; l <= lim; l = r + 1) {
r = std :: min(n_ / (n_ / l), m_ / (m_ / l));
ret = (ret + 1ll * (sum[r] - sum[l - 1] + kMod) * g(n_ / l, m_ / l) % kMod) % kMod;
}
return ret;
}
int Sum(ll l, ll r) {
return (1ll * (r - l + 1ll) * (l + r) / 2ll) % kMod;
}
//=============================================================
int main() {
int n = read(), m = read();
int lim = std :: min(n, m), ans = 0;
Euler(lim);
for (int l = 1, r; l <= lim; l = r + 1) {
r = std :: min(n / (n / l), m / (m / l));
ans = (ans + 1ll * Sum(l, r) * f(n / l, m / l) % kMod) % kMod;
}
printf("%d", ans);
return 0;
}
/*
7718820 8445880
*/