P1829 / SP8099 Crash的数字表格 题解

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\) 的前缀和。

\[f(T) = \sum\limits_{t|T} \mu(t)t^2 \times \dfrac 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;
}