BZOJ 2154: Crash的数字表格

2154: Crash的数字表格
思路：
莫比乌斯反演+整除分块
代码：

#pragma GCC optimize(2)
#pragma GCC optimize(3)
#pragma GCC optimize(4)
#include<bits/stdc++.h>
using namespace std;
#define y1 y11
#define fi first
#define se second
#define pi acos(-1.0)
#define LL long long
//#define mp make_pair
#define pb push_back
#define ls rt<<1, l, m
#define rs rt<<1|1, m+1, r
#define ULL unsigned LL
#define pll pair<LL, LL>
#define pli pair<LL, int>
#define pii pair<int, int>
#define piii pair<pii, int>
#define pdd pair<double, double>
#define mem(a, b) memset(a, b, sizeof(a))
#define debug(x) cerr << #x << " = " << x << "\n";
#define fio ios::sync_with_stdio(false);cin.tie(0);cout.tie(0);
//head
 
const int MOD = 20101009;
const int N = 1e7 + 10;
int prime[N/10], mu[N], cnt;
LL c[N], sum[N];
bool not_p[N];
int n, m;
inline void seive(int n) {
    mu[1] = 1;
    for (int i = 2; i <= n; ++i) {
        if(!not_p[i]) prime[++cnt] = i, mu[i] = -1;
        for (int j = 1; j <= cnt && i*prime[j] <= n; ++j) {
            not_p[i*prime[j]] = true;
            if(i%prime[j] == 0) {
                mu[i*prime[j]] = 0;
                break;
            }
            mu[i*prime[j]] = -mu[i];
        }
    }
    for (int i = 1; i <= n; ++i) sum[i] = (sum[i-1]+mu[i]*i*1LL*i)%MOD, c[i] = (c[i-1]+i)%MOD;
}
inline LL C(int x) {
    return (x*1LL*(x+1)/2) % MOD;
}
inline LL solve(int n, int m) {
    int up = min(n, m);
    LL ans = 0;
    for (int l = 1, r; l <= up; l = r+1) {
        r = min(n/(n/l), m/(m/l));
        ans += (sum[r]-sum[l-1])*C(n/l)%MOD*C(m/l)%MOD;
        ans %= MOD;
    }
    return ans;
}
int main() {
    scanf("%d %d", &n, &m);
    int up = min(n, m);
    seive(up);
    LL ans = 0;
    for (int l = 1, r; l <= up; l = r+1) {
        r = min(n/(n/l), m/(m/l));
        ans += (c[r]-c[l-1])*solve(n/l, m/l);
        ans %= MOD;
    }
    ans = (ans + MOD)%MOD;
    printf("%lld\n", ans);
    return 0;
}