题目传送门 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$了
上代码：
1 #include <cstring> 2 #include <cstdio> 3 #include <algorithm> 4 #define register 5 #define inline 6 7 using namespace std; 8 typedef long long ll; 9 const int MOD = 20101009; 10 const int MAXN = 10000010; 11 12 int n, m; 13 14 bool vis[MAXN] = {0}; 15 int prime[MAXN], pidx = 0, miu[MAXN] = {0}, pf[MAXN] = {0}; 16 void prelude_miu() { 17 miu[1] = 1; 18 for( register int i = 2; i <= n; ++i ) { 19 if( !vis[i] ) { 20 prime[pidx++] = i; 21 miu[i] = -1; 22 } 23 for( register int j = 0; j < pidx; ++j ) { 24 register int k = i * prime[j]; 25 if( k > n ) break; 26 vis[k] = true; 27 if( i % prime[j] ) miu[k] = -miu[i]; 28 else break; 29 } 30 } 31 for( register int i = 1; i <= n; ++i ) // 计算miu[i] * i^2的前缀和，注意取模 32 pf[i] = (pf[i-1] + ((ll)miu[i] * i * i % MOD + MOD) % MOD) % MOD; 33 } 34 35 inline int g( register int n, register int m ) { // 函数g 36 return (ll(n)*(n+1)/2 % MOD) * (ll(m)*(m+1)/2 % MOD) % MOD; 37 } 38 int f( register int n, register int m ) { // 函数f 39 register int ans = 0; 40 for( register int i = 1, j; i <= n; i = j+1 ) { 41 j = min( n/(n/i), m/(m/i) ); // 分块优化 42 // printf( "g(%d,%d) = %lld ", n/j, m/j, g(n/j,m/j) ); 43 ans = (ans + g(n/j,m/j) * ll((pf[j] - pf[i-1] + MOD) % MOD) % MOD) % MOD; 44 } 45 return ans; 46 } 47 void solve() { 48 register int ans = 0; 49 for( register int i = 1, j; i <= n; i = j+1 ) { // 枚举最大公因数i 50 j = min( n/(n/i), m/(m/i) ); // 分块优化 51 // printf( "f(%d,%d) = %lld ", n/j, m/j, f(n/j,m/j) ); 52 ans = (ans + f(n/j,m/j) * (ll(i+j)*(j-i+1)/2 % MOD) % MOD) % MOD; 53 } 54 printf( "%d ", ans ); 55 } 56 57 int main() { 58 scanf( "%d%d", &n, &m ); 59 prelude_miu(); 60 if( n > m ) swap(n,m); // 保证n < m，不用每次取min 61 solve(); 62 return 0; 63 }