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

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

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

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

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

首先，题目要我们求

$Largesumlimits_{i=1}^{n}sumlimits_{j=1}^{m}lcm(i,j)$

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

那么因为

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

所以我们要求的就是

$Largesumlimits_{i=1}^{n}sumlimits_{j=1}^{m}frac{i cdot j}{gcd(i,j)}$

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

$Largesumlimits_{g=gcd(i,j)}frac{i cdot j}{g}$

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

继续变形，得到

$Largesumlimits_{g=gcd(i,j)}frac{i}{g} cdot frac{j}{g} cdot {g}$

然后我们设函数

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

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

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

注意到

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

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

所以最终的答案是

$Largesumlimits_{g=1}^{min(n,m)}f( lfloor frac{n}{g} floor, lfloor frac{m}{g} floor,1) cdot g$

因为$lfloor frac{n}{g} floor$和$lfloor frac{m}{g} floor$都是整除，于是这一步我们就可以分块优化，不会分块优化的请先学习 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)=sumlimits_{x|gcd(i,j), i le a, j le b} frac{i}{x} cdot frac{j}{x}$

于是就有

$Large g(a,b,x)=sumlimits_{i=1}^{ lfloor frac{min(a,b)}{x} floor }f(a,b,i cdot x)$

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

$Large f(a,b,x)=sumlimits_{i=1}^{ lfloor frac{min(a,b)}{x} floor }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} floor, lfloor frac{b}{x} floor,1) cdot x^2$

请读者思考这一步

那么就有

$Large f(a,b,x)=sumlimits_{i=1}^{ lfloor frac{min(a,b)}{x} floor }mu(i)g( lfloor frac{a}{i} floor, lfloor frac{b}{i} floor,x) cdot i^2$

于是$lfloorfrac{a}{i} floor$和$lfloorfrac{b}{i} floor$也都是整除，又可以 前缀和+分块 优化，我们就能在$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} floor, lfloor frac{b}{x} floor,1) cdot x^2=frac{ lfloor frac{a}{x} floor cdot ( lfloor frac{a}{x} floor +1)}{2} cdot frac{ lfloor frac{b}{x} floor cdot ( lfloor frac{b}{x} floor +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/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/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
", ans );
}

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