BZOJ 2154 Crash的数字表格（sigma(lcm(i,j))，莫比乌斯反演）

题目链接；BZOJ 2154 Crash的数字表格
题意：

$$\sum_{i=1}^{i=n} \sum_{j=1}^{j=m}lcm(i,j)\ \%\ 100000009\ \ (n,m\leq 10^{7})$$
分析：
$$ans=\sum_{d=1}^{d=n}\sum_{i=1}^{i=n}\sum_{j=1}^{j=m}\frac{i*j}{d}\ (gcd(i,j)=d)$$
$$=\sum_{d=1}^{d=n}]\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{m}{d}}\frac{i*j*d^2}{d}(gcd(i,j)=1)$$
$$=\sum_{d=1}^{n}d\sum_{i=1}^{\frac{n}{d}}\sum_{j=1}^{\frac{m}{d}}(i*j)(gcd(i,j)=1)$$
$定义f(d)=\sum_{i=1}^{a}\sum_{j=1}^{b}(i*j)\ (gcd(i,j)=d)$
$定义F(d)=\sum_{i=1}^{i=a}\sum_{j=1}^{j=b}(i*j)\ (gcd(i,j)\ \%\ d=0),$
$易得F(1)=sum(a,b)=\frac{a(a+1)}{2}*\frac{b(b+1)}{2}$
$反演可得:f(1)=\sum_{i=1}^{i=a}\sum_{j=1}^{j=b}(i*j)(gcd(i,j)=1,a\leq b)$
$=\sum_{x=1}^{a}\mu(x)*x^2\sum_{i=1}^{\frac{a}{x}}\sum_{j=1}^{\frac{b}{x}}(i*j)$
$=\sum_{x=1}^{a}\mu(x)*x^2\sum_{i=1}^{\frac{a}{x}}i*\sum_{j=1}^{\frac{b}{x}}j$
$=\sum_{x=1}^{a}\mu_(x)*x^2*sum(\frac{a}{x},\frac{b}{x})$
$所以ans = \sum_{d=1}^{n}d\sum_{x=1}^{n'}\mu(x)*x^2*sum(\frac{n'}{x},\frac{m'}{x})(n'=\frac{n}{d},m'=\frac{m}{d})$
此时如果预处理出 $\mu(x) * x^2$的前缀和求 $ans$的复杂度是 $O(\sqrt {n} *\sqrt{n})=O(n)$，对于本题而言是不够的。令 $T= d*x$可得： $ans= \sum_{T=1}^{n}sum(\frac{n}{T},\frac{m}{T})\sum_{x|T}\frac{T}{x}*x^2*\mu(x)$.令 $h[T] = \sum_{x|n}\frac{T}{x}*x^2*\mu(x)$，预处理出 $h[T]$【因为 $h(T)$是积性函数可以线性筛】，那么时间复杂度就变为 $O(\sqrt n)$总的时间复杂度为 $O(Tn)(T是测试组数)$

#include <cstdio>
#include <cstring>
#include <cmath>
#include <string>
#include <algorithm>
#include <bitset>
using namespace std;
typedef long long ll;
const int MAX_N =10000010;
const ll mod = 100000009;

bitset<MAX_N> bs;
int prime_cnt, prime[MAX_N / 100 * 7];
ll h[MAX_N], sum[MAX_N];

void GetMu()
{
    bs.set();
    prime_cnt = 0;
    h[1] = sum[1] = 1;
    for(int i = 2; i < MAX_N; ++i) {
        if(bs[i]) {
            prime[prime_cnt++] = i;
            h[i] = (ll)i * (1 - i) % mod;
        }
        for(int j = 0; j < prime_cnt && i * prime[j] < MAX_N; ++j ){
            bs[i * prime[j]] = 0;
            if(i % prime[j]) {
                h[i * prime[j]] = h[i] * h[prime[j]] % mod;
            }else {
                h[i * prime[j]] = h[i] * prime[j] % mod;
                break;
            }
        }
    }
    for(int i = 1; i < MAX_N; ++i ) {
        sum[i] = ((sum[i - 1] + h[i]) % mod + mod) % mod;
    }
}

inline ll  work(int n, int  m) 
{
    ll res1 = (ll) n * (n + 1) / 2 % mod;
    ll res2 = (ll) m * (m + 1) / 2 % mod;
    return res1 * res2 % mod;   
}

inline ll solve(int n, int m)
{
    int top = min(n, m), last;
    ll res = 0;
    for(int i = 1; i <= top; i = last + 1) {
        last = min( n / (n / i), m / (m / i));
        res = (res + (sum[last] - sum[i - 1] + mod) % mod * work(n / i, m / i) % mod) % mod;
    }
    return res;
}

int main()
{
    GetMu();
    int T,  n, m;
    scanf("%d", &T);
    while(T--) {
        scanf("%d%d", &n, &m);
        printf("%lld\n", solve(n, m));
    }
    return 0;
}