BZOJ 2154 Crash的數字表格：雙重求和 LCM 莫比烏斯反演解法

問題重述

計算雙重求和：S = Σᵢ₌₁ⁿ Σⱼ₌₁ᵐ lcm(i, j)，對大數取模 (通常模 20101009 或類似質數)。n, m 可達 10⁷ 量級，需 O(√n) 或 O(n^(2/3)) 算法。

數學轉換核心

利用 lcm(i,j) = i*j / gcd(i,j)，令 d = gcd(i,j)，則 i = d*i', j = d*j' 且 gcd(i',j')=1：

S = Σ_d d * Σ_{i'≤n/d, j'≤m/d, gcd=1} i' * j'

內層求和可用莫比烏斯函數 μ 通過容斥展開：

Σ_{gcd=1} i'j' = Σ_k μ(k) * k² * (Σ_{i''} i'') * (Σ_{j''} j'')

其中 i'' ≤ n/(d*k), j'' ≤ m/(d*k)，前綴和為等差數列求和公式。

算法優化策略

1. 線性篩預計算 μ(k) 及其前綴和 M(x) = Σ_{i≤x} μ(i)，O(N)。
2. 整除分塊：對外層 d 和內層 k 合併枚舉 t = d*k，利用 n/t 和 m/t 只會有 O(√n) 種取值的性質，將雙重循環優化為單循環分塊處理。
3. 預計算前綴和：S₁(x) = Σ i, S₂(x) = Σ i² 等，O(1) 獲取區間和。

關鍵代碼框架 (C++)

#include <bits/stdc++.h>
using namespace std;
typedef long long ll;
const int MOD = 20101009;
const int MAXN = 10000000;

int mu[MAXN+5], prime[MAXN+5], cnt;
bool vis[MAXN+5];
ll sum_mu[MAXN+5]; // μ 前綴和

void linear_sieve(int n) {
    mu[1] = 1;
    for (int i=2; i<=n; ++i) {
        if (!vis[i]) prime[++cnt]=i, mu[i]=-1;
        for (int j=1; j<=cnt && i*prime[j]<=n; ++j) {
            vis[i*prime[j]] = true;
            if (i % prime[j] == 0) { mu[i*prime[j]] = 0; break; }
            else mu[i*prime[j]] = -mu[i];
        }
    }
    for (int i=1; i<=n; ++i) sum_mu[i] = (sum_mu[i-1] + mu[i]) % MOD;
}

ll S1(ll n) { return n*(n+1)/2 % MOD; } // Σ i
ll S2(ll n) { return n*(n+1)%MOD*(2*n+1)%MOD*16666667%MOD; } // Σ i², 16666667 = 6⁻¹ mod MOD

int main() {
    int n, m; cin >> n >> m;
    if (n > m) swap(n, m);
    linear_sieve(n);
    
    ll ans = 0;
    for (int l=1, r; l<=n; l=r+1) {
        int nd = n/l, md = m/l;
        r = min(n/nd, m/md);
        // 這裡需要結合 μ 的前綴和與 S1/S2 計算塊貢獻
        // 具體推導較複雜，核心是利用 sum_mu[r] - sum_mu[l-1]
        // 乘以對應的 S1(nd)*S1(md) 等項
    }
    cout << (ans%MOD+MOD)%MOD << endl;
    return 0;
}

複雜度分析

• 線性篩：O(n) 時間，O(n) 空間。
• 主循環分塊：O(√n) 次迭代。
• 總體：O(n) 預處理 + O(√n) 查詢，可處理 n,m ≤ 10⁷。

關鍵難點提示

1. 模數 20101009 是質數，可用費馬小定理求模逆元 (如 6⁻¹)。
2. 雙重求和交換順序後，需仔細處理邊界條件與取模運算。
3. μ 函數前綴和分塊是優化關鍵，類似杜教篩思想但更簡單。
4. 注意 long long 溢出，及時取模。

此題是數論、莫比烏斯反演、整除分塊綜合應用的經典範例，掌握其推導過程對解決類似雙重求和 GCD/LCM 問題極有價值。