51nod1238. 最小公倍數之和 V3(莫比烏斯反演)
題目連結
https://www.51nod.com/Challenge/Problem.html#!#problemId=1238
題解
本來想做個杜教篩板子題結果用另一種方法過了......
所謂的“另一種方法”用到的技巧還是挺不錯的,因此這裡簡單介紹一下。
首先還是基本的推式子:
\[\begin{aligned}\sum_{i = 1}^n \sum_{j = 1}^n {\rm lcm}(i, j) &= \sum_{i = 1}^n \sum_{j = 1}^n \frac{ij}{{\rm gcd}(i, j)} \\ &= \sum_{d = 1}^{n} \sum_{i = 1}^n \sum_{j = 1}^n[{\rm gcd}(i, j) = d]\frac{ij}{d} \\ &= \sum_{d = 1}^n d\sum_{i = 1}^{\left\lfloor\frac{n}{d}\right\rfloor}\sum_{j = 1}^{\left\lfloor\frac{n}{d}\right\rfloor}[{\rm gcd}(i, j) = 1]ij\end{aligned}\]
設 \(f(x) = \sum_\limits{i = 1}^x \sum_\limits{j = 1}^x [{\rm gcd}(i, j) = 1]ij\),那麼答案即為 \(\sum_\limits{d = 1}^n d\times f(\left\lfloor\frac{n}{d}\right\rfloor)\)。顯然答案可以數論分塊求,因此我們的任務就是求 \(f\) 函式。
首先考慮當 \(x\) 較小時,我們能否直接預處理出 \(f(x)\)。考慮差分:當 \(x > 1\) 時,\(f(x)\) 較 \(f(x - 1)\) 而言,多了的部分為:\[\left(\sum_\limits{i = 1}^x\sum_\limits{j = 1}^x [{\rm gcd}(i, j) = 1]ij\right) - \left(\sum_\limits{i = 1}^{x - 1}\sum_\limits{j = 1}^{x - 1} [{\rm gcd}(i, j) = 1]ij\right) = 2x \sum_{i = 1}^x [{\rm gcd}(i, x) = 1]i\]
而由於小於等於 \(x(x > 1)\) 且與 \(x\) 互質的數的和為 \(\frac{\varphi(x)x}{2}\)(證明提示:當 \(n \geq 2\) 時,若 \({\rm gcd}(d, n) = 1\) 必然有 \({\rm gcd}(n - d, n) = 1\),與 \(n\) 互質的 \(d\) 共有 \(\varphi(n)\) 個),因此我們就得到了 \(f(x)\) 的遞推式:\(f(x) = f(x - 1) + 2x \times \frac{\varphi(x)x}{2}\),即 \(f(x) = f(x - 1) + \varphi(x)x^2\)
不過這隻能處理 \(x\) 較小的情況。當 \(x\) 較大時,我們仍然得另謀他路。在 \(f(x) = \sum_\limits{i = 1}^x \sum_\limits{j = 1}^x [{\rm gcd}(i, j) = 1]ij\) 當中,對 \(f(x)\) 有貢獻的 \(i, j\) 滿足 \(i\) 與 \(j\) 是互質的,我們考慮補集轉化,用總和減去不互質的 \(i, j\) 的貢獻:
\[\begin{aligned} f(x) &= \sum_{i = 1}^x\sum_{j = 1}^x ij - \sum_{d = 2}^x \sum_{i = 1}^x \sum_{j = 1}^x [{\rm gcd}(i, j) = d]ij \\ &= \left(\frac{x(x+1)}{2}\right)^2 - \sum_{d = 2}^x d^2 \sum_{i = 1}^{\left\lfloor\frac{x}{d}\right\rfloor} \sum_{j = 1}^{\left\lfloor\frac{x}{d}\right\rfloor} [{\rm gcd}(i, j) = 1]ij \\ &=\left(\frac{x(x+1)}{2}\right)^2 - \sum_{d = 2}^x d^2 f(\left\lfloor\frac{x}{d}\right\rfloor) \end{aligned}\]
這樣,我們就能夠遞迴地去求解當 \(x\) 較大時 \(f(x)\) 的值了。不難發現,該求解方法的時間複雜度和杜教篩是一樣的,為 \(O(n^{\frac{2}{3}})\),且非常好寫。
程式碼
#include<bits/stdc++.h>
using namespace std;
const int mod = 1000000007, inv2 = 500000004, inv6 = 166666668, up = 10000001;
int main() {
function<int (int, int)> mul = [&] (int x, int y) {
return (long long) x * y % mod;
};
function<void (int&, int)> add = [&] (int& x, int y) {
x += y;
if (x >= mod) {
x -= mod;
}
};
function<void (int&, int)> sub = [&] (int& x, int y) {
x -= y;
if (x < 0) {
x += mod;
}
};
vector<bool> is_prime(up, true);
vector<int> phi(up), primes;
phi[1] = 1;
for (int i = 2; i < up; ++i) {
if (is_prime[i]) {
primes.push_back(i);
phi[i] = i - 1;
}
for (auto v : primes) {
int d = v * i;
if (d >= up) {
break;
}
is_prime[d] = false;
if (i % v == 0) {
phi[d] = mul(phi[i], v);
break;
} else {
phi[d] = mul(phi[i], phi[v]);
}
}
}
for (int i = 2; i < up; ++i) {
phi[i] = mul(mul(phi[i], i), i);
add(phi[i], phi[i - 1]);
}
function<int (long long)> sum_pow2 = [&] (long long n) {
n %= mod;
return mul(mul(mul(n, n + 1), (n * 2 + 1)), inv6);
};
map<long long, int> value;
function<int (long long)> f = [&] (long long n) {
if (value.count(n)) {
return value[n];
} else {
int result = 0;
if (n < up) {
result = phi[n];
} else {
int x = n % mod;
x = mul(mul(x, x + 1), inv2);
result = mul(x, x);
for (long long i = 2, last; i <= n; i = last + 1) {
last = n / (n / i);
sub(result, mul((sum_pow2(last) - sum_pow2(i - 1) + mod) % mod, f(n / i)));
}
}
return value[n] = result;
}
};
long long n;
scanf("%lld", &n);
int answer = 0;
for (long long i = 1, last; i <= n; i = last + 1) {
last = n / (n / i);
add(answer, mul(mul(mul((i + last) % mod, (last - i + 1) % mod), inv2), f(n / i)));
}
printf("%d\n", answer);
return 0;
}