【题解】Crash的数字表格 BZOJ 2154 莫比乌斯反演

题目传送门 http://www.lydsy.com/JudgeOnline/problem.php?id=2154

人生中第一道自己做出来的莫比乌斯反演

人生中第一篇用LaTeX写数学公式的博客

大家别看公式多就害怕了啊，这里面的公式大多是很显然的

首先，题目要我们求

$\Large\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}lcm(i,j)$

这个$lcm$很不好办，我们想办法转化成$gcd$，然后尝试搞莫比乌斯反演的套路

那么因为

$\Large lcm(i,j)=\frac{i \cdot j}{gcd(i,j)}$

所以我们要求的就是

$\Large\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\frac{i \cdot j}{gcd(i,j)}$

于是我们可以枚举每个最大公约数$g$，对于每个$g$，我们求

$\Large\sum\limits_{g=gcd(i,j)}\frac{i \cdot j}{g}$

最后把这些都加起来就好了

继续变形，得到

$\Large\sum\limits_{g=gcd(i,j)}\frac{i}{g} \cdot \frac{j}{g} \cdot {g}$

然后我们设函数

$\Large f(a,b,x)=\sum\limits_{x=gcd(i,j),i \le a,j \le b}\frac{i}{x} \cdot \frac{j}{x}$

于是我们要最终的答案就是

$\Large \sum\limits_{g=1}^{min(n,m)}f(n,m,g) \cdot g$

注意到

$\Large f(a,b,x)=f( \lfloor \frac{a}{x} \rfloor, \lfloor \frac{b}{x} \rfloor,1)$

这一步的原因请读者自行思考

所以最终的答案是

$\Large\sum\limits_{g=1}^{min(n,m)}f( \lfloor \frac{n}{g} \rfloor, \lfloor \frac{m}{g} \rfloor,1) \cdot g$

因为$\lfloor \frac{n}{g} \rfloor$和$\lfloor \frac{m}{g} \rfloor$都是整除，于是这一步我们就可以分块优化，不会分块优化的请先学习 http://sxysxy.org/blogs/77

既然有了分块优化，这一步就可以在$O(\sqrt{n})$的时间内解决，我们只要在$O(\sqrt{n})$的时间内求出来$f(a,b,1)$就好了

我们发现$f(a,b,1)$是和$gcd$相关的函数，于是可以很自然的想到莫比乌斯反演的套路，不会套路的请先学习 http://sxysxy.org/blogs/77

按照套路，我们设

$\Large g(a,b,x)=\sum\limits_{x|gcd(i,j), i \le a, j \le b} \frac{i}{x} \cdot \frac{j}{x}$

于是就有

$\Large g(a,b,x)=\sum\limits_{i=1}^{ \lfloor \frac{min(a,b)}{x} \rfloor }f(a,b,i \cdot x)$

于是按照套路，我们有莫比乌斯反演

$\Large f(a,b,x)=\sum\limits_{i=1}^{ \lfloor \frac{min(a,b)}{x} \rfloor }\mu(i)g(a,b,i \cdot x)$

不会这个套路的请学习 http://sxysxy.org/blogs/77，或者学习WC2016课件，这个也可以从容斥的角度理解

同样对于$g(a,b,x)$，我们有

$\Large g(a,b,x)=g( \lfloor \frac{a}{x} \rfloor, \lfloor \frac{b}{x} \rfloor,1) \cdot x^2$

请读者思考这一步

那么就有

$\Large f(a,b,x)=\sum\limits_{i=1}^{ \lfloor \frac{min(a,b)}{x} \rfloor }\mu(i)g( \lfloor \frac{a}{i} \rfloor, \lfloor \frac{b}{i} \rfloor,x) \cdot i^2$

于是$\lfloor\frac{a}{i}\rfloor$和$\lfloor\frac{b}{i}\rfloor$也都是整除，又可以 前缀和+分块 优化，我们就能在$O(\sqrt{n})$的时间内求出每个$f(a,b,x)$啦

于是我们只要在$O(1)$的时间内求出$g(a,b,x)$，总时间复杂度就是$O(n)$啦！

显然

$\Large g(a,b,x)=g( \lfloor \frac{a}{x} \rfloor, \lfloor \frac{b}{x} \rfloor,1) \cdot x^2=\frac{ \lfloor \frac{a}{x} \rfloor \cdot ( \lfloor \frac{a}{x} \rfloor +1)}{2} \cdot \frac{ \lfloor \frac{b}{x} \rfloor \cdot ( \lfloor \frac{b}{x} \rfloor +1)}{2} \cdot x^2$

为什么显然请读者思考

于是任务大功告成！总时间复杂度为$O(n)$

实际上对于函数$g$和$f$，最后的那个参数$x$是没必要传的，因为在我们的分析中，$x$都被化作$1$了

上代码：

#include <cstring>
#include <cstdio>
#include <algorithm>
#define register
#define inline

using namespace std;
typedef long long ll;
const int MOD = 20101009;
const int MAXN = 10000010;

int n, m;

bool vis[MAXN] = {0};
int prime[MAXN], pidx = 0, miu[MAXN] = {0}, pf[MAXN] = {0};
void prelude_miu() {
    miu[1] = 1;
    for( register int i = 2; i <= n; ++i ) {
        if( !vis[i] ) {
            prime[pidx++] = i;
            miu[i] = -1;
        }
        for( register int j = 0; j < pidx; ++j ) {
            register int k = i * prime[j];
            if( k > n ) break;
            vis[k] = true;
            if( i % prime[j] ) miu[k] = -miu[i];
            else break;
        }
    }
    for( register int i = 1; i <= n; ++i ) // 计算miu[i] * i^2的前缀和，注意取模
        pf[i] = (pf[i-1] + ((ll)miu[i] * i * i % MOD + MOD) % MOD) % MOD;
}

inline int g( register int n, register int m ) { // 函数g
    return (ll(n)*(n+1)/2 % MOD) * (ll(m)*(m+1)/2 % MOD) % MOD;
}
int f( register int n, register int m ) { // 函数f
    register int ans = 0;
    for( register int i = 1, j; i <= n; i = j+1 ) {
        j = min( n/(n/i), m/(m/i) ); // 分块优化
        // printf( "g(%d,%d) = %lld\n", n/j, m/j, g(n/j,m/j) );
        ans = (ans + g(n/j,m/j) * ll((pf[j] - pf[i-1] + MOD) % MOD) % MOD) % MOD;
    }
    return ans;
}
void solve() {
    register int ans = 0;
    for( register int i = 1, j; i <= n; i = j+1 ) { // 枚举最大公因数i
        j = min( n/(n/i), m/(m/i) ); // 分块优化
        // printf( "f(%d,%d) = %lld\n", n/j, m/j, f(n/j,m/j) );
        ans = (ans + f(n/j,m/j) * (ll(i+j)*(j-i+1)/2 % MOD) % MOD) % MOD;
    }
    printf( "%d\n", ans );
}

int main() {
    scanf( "%d%d", &n, &m );
    prelude_miu();
    if( n > m ) swap(n,m); // 保证n < m，不用每次取min
    solve();
    return 0;
}