[国家集训队]Crash的数字表格 / JZPTAB 莫比乌斯反演

～～～题面～～～

题解：
　　$$ans = sum_{i = 1}^{n}sum_{j = 1}^{m}{frac{ij}{gcd(i, j)}}$$
　　改成枚举d(设n < m)
　　$$ans = sum_{d = 1}^{n}sum_{i = 1}^{n}sum_{j = 1}^{m}[gcd(i, j) == d]frac{ij}{d}$$
　　考虑枚举$id$
　　设$N = lfloor{frac{n}{d}} floor$,$M = lfloor{frac{m}{d}} floor$
　　$$ans = sum_{d = 1}^{n}{d}sum_{i = 1}^{N}sum_{j = 1}^{M}sum_{t | gcd(i, j)}{mu(t)ij}$$
　　把后面改成枚举$mu(t)$
　　$$ans = sum_{d = 1}^{n}dsum_{t = 1}^{N}{mu(t)} sum_{i = 1}^{lfloor{frac{N}{t}} floor}sum_{j = 1}^{lfloor{frac{M}{t}} floor}ijt^2$$
　　$$ans = sum_{d = 1}^{n}dsum_{t = 1}^{N}{mu(t)} t^2 sum_{i = 1}^{lfloor{frac{N}{t}} floor}sum_{j = 1}^{lfloor{frac{M}{t}} floor}ij$$
　　$$ans = sum_{d = 1}^{n}dsum_{t = 1}^{N}{mu(t)} t^2 sum_{i = 1}^{lfloor{frac{N}{t}} floor}isum_{j = 1}^{lfloor{frac{M}{t}} floor}j$$
　　因为$sum_{i = 1}^{lfloor{frac{N}{t}} floor}i$显然是可以$O(n)$预处理的，而$lfloor{frac{N}{t}} floor lfloor{frac{M}{t}} floor$可以整数分块，$t^2mu(t)$可以线性筛$O(n)$预处理，总复杂度$O(n + nsqrt{n})$

　　不过其实还有更优的方法，以后再说吧。。。

#include<bits/stdc++.h>
using namespace std;
#define R register int
#define AC 10000010
#define p 20101009
#define LL long long

int n, m, N, M, tot;
LL ans;
int prime[AC], mu[AC];
LL sum[AC], s[AC];//质数，mu函数，1~n求和，mu(t)*t^2求和
bool z[AC];

void pre()
{
    scanf("%d%d", &n, &m);
    mu[1] = 1;
    int b = max(n, m), x;
    s[1] = sum[1] = 1;//初始化1，，，，因为下面的i是从2开始的
    for(R i = 2; i <= b; i++)
    {
        if(!z[i]) prime[++tot] = i, mu[i] = -1;
        s[i] = (s[i - 1] + (LL)mu[i] * ((LL)i * (LL)i)%p) % p;//error !!, mu要LL
        sum[i] = (sum[i - 1] + i) % p;
        for(R j = 1; j <= tot; j ++)
        {
            x = prime[j];
            if(i * x > b) break;
            z[i * x] = true;
            if(!(i % x)) break;
            mu[i * x] = -mu[i];
        }
    }
}

void work()
{
    if(m < n) swap(n, m);
    //for(R i = 1; i <= m; i ++) printf("%lld ", sum[i]);
    //printf("
");
    LL summ;
    for(R i = 1; i <= n; i++)//枚举d
    {
        int pos;
        N = n / i, M = m / i;
        summ = 0;
        for(R j = 1; j <= N; j = pos + 1)//枚举t
        {
            pos = min(N / (N / j), M / (M / j));
            summ += (((s[pos] - s[j - 1]) % p * sum[N / j]) % p * sum[M / j]) % p;
            summ %= p;
        }
        summ = (summ * i + p) % p;//error!!!别把d给忘了！
        ans = (ans + summ + p) % p;
        //printf("%lld
", ans);
    }
    printf("%lld
", ans);
}

int main()
{
//    freopen("in.in", "r", stdin);
    pre();
    work();
//    fclose(stdin);
    return 0;
}