【题解】[国家集训队]Crash的数字表格 / JZPTAB

　　求解(sum_{i = 1}^{n}sum_{j = 1}^{m}lcmleft ( i,j ight ))。

有(lcmleft ( i,j ight )=frac{ij}{gcdleft ( i,j ight )})，

所以原本的式子转化为：(sum_{i = 1}^{n}sum_{j = 1}^{m}frac{ij}{gcdleft ( i,j ight )})。

注意到(i, j) 均为 (gcdleft ( i,j ight )) 的倍数，且原式中有除法不好处理，

所以我们改为枚举(gcdleft ( i,j ight )) 的倍数。

有：(sum_{d = 1}^{n}  d sum_{i = 1}^{frac{n}{d}}sum_{j = 1}^{frac{m}{d}}ijleft [ gcdleft ( i,j ight ) = 1 ight])。

后面的式子套路的来一发反演：

(sum_{d = 1}^{n}  d sum_{i = 1}^{frac{n}{d}}sum_{j = 1}^{frac{m}{d}}ijsum_{k|gcdleft ( i,j ight )}mu left ( k ight ))

注意这里面有一个乘积的项，可以理解为是任意数字的两两匹配，即：

(sum_{i = 1}^{n}sum_{j = 1}^{m}ij = left ( 1 + 2 + ... + n ight )left ( 1 + 2 + ... + m ight ))

所以转化为：

(sum_{d = 1}^{n}  d sum_{k = 1}^{frac{n}{d}} k^{2} * mu left ( k ight )sumleft ( left lfloor frac{n}{dk} ight floor ight )sumleft ( left lfloor frac{m}{dk} ight floor ight ))

依然是套路的改变枚举项为 (dk)

(sum_{T = 1}^{n}  sumleft ( left lfloor frac{n}{T} ight floor ight )sumleft ( left lfloor frac{m}{T} ight floor ight ) * T sum_{d|T}d*mu left ( d ight ))

　　到这里我们已经实现了第一步：前面的部分可以数论分块(Oleft ( sqrt{n} ight ))处理，只要我们能够通过线性筛处理出后面的一部分，这道题目就完成了。为了实现线性筛，我们对于后面部分进行观察。我们令(F[T] = T * sum_{d|T}d*mu left ( d ight )) 。

　　首先，(F[i])当 (i) 为质数时，(F[i]) 的值很容易确定为 (i - i^{2})。 注意到它实际上是积性函数。所以在线性筛中若 (i = x * y) ，（其中 (x) 为 (i) 的最小质因子），当 (y   mod x eq  0) 时说明二者互质，则 (F[i] = F[x] * F[y])。

　　然后考虑当(y   mod x =  0)的情况，这说明这两个部分中均含有最小的质因子。注意因为卷入了一个 (mu)，所以有平方因子时的值都不会造成贡献。也就是说取值范围和 (y) 仍然是相同的，只不过是系数改变了。所以此时 (F[i] = F[y] * x )。然后此题就圆满解决啦~~~

#include <bits/stdc++.h>
using namespace std;
#define maxn 10005000
#define int long long
#define mod 20101009
int n, m, N, maxx = maxn - 1e3;
int tot, pri[maxn], inv2;
int ans, f[maxn];
bitset <maxn> is_prime;

int read()
{
    int x = 0, k = 1;
    char c;
    c = getchar();
    while(c < '0' || c > '9') { if(c == '-') k = -1; c = getchar(); }
    while(c >= '0' && c <= '9') x = x * 10 + c - '0', c = getchar();
    return x * k;
}

int qpow(int x, int times) 
{
    int base = 1;
    for(; times; times >>= 1, x = (x * x) % mod)
        if(times & 1) base = (base * x) % mod;
    return base;
}
int Sum(int x) { x %= mod; return ((x * (x + 1)) % mod * inv2 % mod);}

void Get_F()
{
    f[1] = 1;
    for(int i = 2; i <= maxx; i ++)
    {
        if(!is_prime[i]) pri[++ tot] = i, f[i] = i * (1ll - i) % mod;
        for(int j = 1; j <= tot; j ++)
        {
            int K = i * pri[j]; if(K > maxx) break;
            is_prime[K] = 1;
            if(!(i % pri[j])) { f[K] = f[i] * pri[j] % mod; break; }
            else f[K] = f[i] * f[pri[j]] % mod;
        }
    }
    for(int i = 1; i <= maxx; i ++) f[i] = (f[i] + f[i - 1]) % mod;
}

signed main()
{
    n = read(), m = read(), N = min(n, m);
    maxx = min(n, m); inv2 = qpow(2, mod - 2);
    Get_F();
    for(int l = 1, r; l <= N; l = r + 1)
    {
        r = min((n / (n / l)), (m / (m / l)));
        int ret = Sum(n / l) * Sum(m / l) % mod;
        ans = (ans + (ret * (f[r] - f[l - 1]) % mod)) % mod; 
    }
    printf("%lld
", (ans + mod) % mod);
    return 0;
}