P1829 / SP8099 Crash的数字表格 题解
Published:
P1829 / SP8099 Crash的数字表格 题解
莫比乌斯反演、数论分块、杜教筛。
题意简述
求以下式子的值，对 $20101009$（一个质数）取模：
\[\sum\limits_{i=1}^n \sum\limits_{j=1}^m \operatorname{lcm}(i,j)\]$n,m \le 10^7$。加强：$n,m \le 10^{10}$。
解法
莫比乌斯反演
设 $s_1(n) = \sum\limits_{i=1}^n i = \dfrac {n (n+1)} 2$，$s_2(n) = \sum\limits_{i=1}^n i^2 = \dfrac {n (n+1) (2n+1)} 6$。
则答案可以通过莫比乌斯反演得到（这一段可以看其他题解），为：
\[\sum\limits_{T=1}^n s_1\left(\left\lfloor \dfrac n T \right\rfloor\right) \times s_1\left(\left\lfloor \dfrac m T \right\rfloor\right) \times T \times \sum\limits_{t|T} \mu(t)t\]直接做可以 $O(n \log \log n)$ 或者 $O(n)$，但我们显然是不满意的。
杜教筛
| 两个下取整除法可以数论分块（瓶颈显然不在这），所以现在问题变为求 $f(T) = T \times \sum\limits_{t | T} \mu(t)t$ 的前缀和。 |
设函数 $\text{Id}(n) = n$，$\text{Id}_2(n) = n^2$。
\[f = (\mu \cdot \text{Id}_2) * \text{Id}\]观察其贝尔级数：
\[f_p(x) = \dfrac{1 - p^2 x} {1 - px}\]不难发现：
\[\dfrac{1 - p^2 x} {1 - px} \times \dfrac 1 {1 - p^2 x} = \dfrac 1 {1 - p x}\] \[f * \text{Id}_2 = \text{Id}\]杜教筛即可。时间复杂度 $O(n^{\frac 2 3})$。
代码
代码是 Luogu P1829 次优解（截止至 2024/3/1），注意提交到 SP8099 需要将语言改为 C++98（当然 auto 之类的也要一起改掉）。
#include <bits/stdc++.h>
#define umap unordered_map
using namespace std;
typedef long long ll;
const ll MOD = 20101009;
const ll inv2 = 10050505;
const ll inv6 = 16750841;
const ll CBRT2 = 46415 + 100;
const ll MAXN2 = CBRT2 + 100;
ll s1(ll x) {
return x * (x + 1) % MOD * inv2 % MOD;
}
ll s2(ll x) {
return x * (x + 1) % MOD * (2 * x + 1) % MOD * inv6 % MOD;
}
ll n, m;
bool vis[MAXN2];
vector<ll> primes;
ll mu[MAXN2];
ll f[MAXN2], sum[MAXN2];
void sieve(ll n) {
vis[1] = 1;
mu[1] = 1;
f[1] = 1;
for (ll i = 2; i <= n; i++) {
if (!vis[i]) {
primes.push_back(i);
mu[i] = -1;
f[i] = i - i * i;
}
for (auto j : primes) {
ll k = i * j; if (k > n) break;
vis[k] = 1;
if (i % j == 0) {
mu[k] = 0;
f[k] = f[i] * j;
break;
}
mu[k] = mu[i] * mu[j];
f[k] = f[i] * f[j];
}
}
for (ll i = 1; i <= n; i++)
sum[i] = (sum[i-1] + f[i]) % MOD;
}
umap<ll, ll> s;
ll S(ll x) {
if (x <= CBRT2)
return sum[x];
if (s.count(x))
return s[x];
ll ans = s1(x);
for (ll l = 2, r; l <= x; l = r + 1) {
ll v = x / l;
r = x / v;
ans -= (s2(r) - s2(l-1)) * S(v) % MOD;
}
(ans += MOD) %= MOD;
return s[x] = ans;
}
ll solve(ll n, ll m) {
ll ans = 0;
for (ll l = 1, r; l <= n; l = r + 1) {
ll v1 = n / l, v2 = m / l;
r = min(n / v1, m / v2);
(ans += s1(v1) * s1(v2) % MOD * (S(r) - S(l-1) + MOD) % MOD) %= MOD;
}
return ans;
}
int main() {
cin >> n >> m;
if (n > m) swap(n, m);
sieve(min(n, CBRT2));
cout << solve(n, m) << '\n';
return 0;
}