BZOJ 2154 莫比乌斯反演 Crash的数字表格

BZOJ 2154 莫比乌斯反演 Crash的数字表格

题意：${\sum }_{i=1}^n{\sum }_{j=1}^{m}lcm(i,j)$

#include <bits/stdc++.h>
#define FOR(i,s,t) for(int i=(s);i<=(t);i++)
#define ROF(i,s,t) for(int i=(s);i>=(t);i--)
#define pb push_back
#define mp make_pair
#define eb emplace_back
#define fi first
#define se second
#define endl '\n'
using namespace std;
typedef unsigned long long ull;
typedef long long ll;
const int maxn = 1e7 + 6;
const ll mod = 20101009;
const int inf = 0x3f3f3f3f;
const ll INF = 0x3f3f3f3f3f3f3f3f;
int readInt(){
    int x=0;
    bool sign=false;
    char c=getchar();
    while(!isdigit(c)){
        sign=c=='-';
        c=getchar();
    }
    while(isdigit(c)){
        x=x*10+c-'0';
        c=getchar();
    }
    return sign?-x:x;
}
ll readLong(){
    ll x=0;
    bool sign=false;
    char c=getchar();
    while(!isdigit(c)){
        sign=c=='-';
        c=getchar();
    }
    while(isdigit(c)){
        x=x*10+c-'0';
        c=getchar();
    }
    return sign?-x:x;
}
string readString(){
    string s;
    char c=getchar();
    while(isspace(c)){
        c=getchar();
    }
    while(!isspace(c)){
        s+=c;
        c=getchar();
    }
    return s;
}

vector <int> pri;
int vis[maxn];
int mu[maxn];
ll n, m;
ll sum[maxn];

void getPrime(int maxn){
    maxn++;
    mu[1] = 1 ;
    for(int i = 2; i < maxn; i++){
        if (!vis[i]){
            pri.pb(i);
            mu[i] = -1;
        }
        for (int j = 0; j < pri.size() && i * pri[j] <maxn; j++){
            vis[i *pri[j]] = 1;
            if (i % pri[j] == 0){
                mu[i * pri[j]] = 0;
                break;
            }
            mu[i *pri[j]] = - mu[i];
        }
    }
    for (int i = 1; i <maxn; i++ ){
        sum[i] = (sum[i-1] % mod + (((ll)mu[i] * (ll)i) % mod * i) % mod) % mod;
        //cout << sum[i] << endl;
    }
}

ll getsum(ll n, ll m){
    return ((((n + 1) * n / 2) % mod) * (((m + 1) * m / 2)% mod)) % mod;
}
ll solve (ll x, ll y){
    ll res = 0;
    int pos = 1;
    if (x > y) swap(x, y);
    for (int i = 1; i <= x; i = pos + 1){
        pos = min(x/(x/i), y/(y/i));
        res =(res + (sum[pos] - sum[i-1]) % mod * getsum(x/i, y/i) %mod) % mod;
    }
    return res % mod;
}

int main(){
    n = readLong();
    m = readLong();
    if (n > m ) swap(n, m);
    getPrime(n);
    int pos = 0;
    ll res = 0;
    for (int i= 1; i <= n; i= pos +1){
        pos = min(n/(n/i), m/(m/i));
        res =(res + ((ll)(pos + i) * (ll)(pos - i + 1) / 2ll % mod) * (solve(n/ i, m / i) % mod)) % mod;
        //cout << res << endl;
    }
    printf("%lld\n", (res + mod) % mod);
    return 0;
}