洛谷P1829 [国家集训队]Crash的数字表格 / JZPTAB

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

思路：

$\sum_{i = 1}^{n}\sum_{j = 1}^{m}lcm(i, j)$

$= \sum_{i = 1}^{n}\sum_{j = 1}^{m}\frac{i * j}{\gcd(i, j)}$

$= \sum_{d}\sum_{i = 1}^{n}\sum_{j = 1}^{m}\frac{i * j}{d}[\gcd(i, j)=d]$

$= \sum_{d}\sum_{i = 1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j = 1}^{\lfloor\frac{m}{d}\rfloor}i * j * d[\gcd(i, j)=1]$

$= \sum_{d}\sum_{i = 1}^{\lfloor\frac{n}{d}\rfloor}\sum_{j = 1}^{\lfloor\frac{m}{d}\rfloor}i * j * d \sum_{k | gcd(i, j)} mu(k)$

$= \sum_{d} d \sum_{k} mu(k) \sum_{i = 1}^{\lfloor\frac{n}{d}\rfloor}[k|i]*i \sum_{j = 1}^{\lfloor\frac{m}{d}\rfloor}[k|j] * j$

$= \sum_{d} d \sum_{k} mu(k) \sum_{i = 1}^{\lfloor\frac{n}{dk}\rfloor}k*i \sum_{j = 1}^{\lfloor\frac{m}{dk}\rfloor}k * j$

$= \sum_{d} d \sum_{k} k^2 mu(k) \sum_{i = 1}^{\lfloor\frac{n}{dk}\rfloor} i \sum_{j = 1}^{\lfloor\frac{m}{dk}\rfloor} j$

 

Code:

#pragma GCC optimize(3)
#pragma GCC optimize(2)
#include <map>
#include <set>
// #include <array>
#include <cmath>
#include <queue>
#include <stack>
#include <vector>
#include <cstdio>
#include <cstring>
#include <sstream>
#include <iostream>
#include <stdlib.h>
#include <algorithm>
// #include <unordered_map>

using namespace std;

typedef long long ll;
typedef pair<int, int> PII;

#define Time (double)clock() / CLOCKS_PER_SEC

#define sd(a) scanf("%d", &a)
#define sdd(a, b) scanf("%d%d", &a, &b)
#define slld(a) scanf("%lld", &a)
#define slldd(a, b) scanf("%lld%lld", &a, &b)

const int N = 1e7 + 20;
const int M = 1e5 + 20;
const int mod = 20101009;
const double eps = 1e-6;

ll cnt = 0, primes[N], mu[N], sum[N] = {0}, sd[N] = {0};
map<ll, ll> su;
bool st[N];

int n, m, p, k;

void get(ll n){
    mu[1] = 1;
    d[1] = 1;
    for(ll i = 2; i <= n; i ++){
        if(!st[i]){
            primes[cnt ++] = i;
            mu[i] = -1;
        }
        for(int j = 0; primes[j] <= n / i; j ++){
            st[i * primes[j]] = true;
            if(i % primes[j] == 0){
                mu[i * primes[j]] = 0;
                break;
            }
            mu[i * primes[j]] = - mu[i];
        }
    }

    for(int i = 1; i <= n; i ++){
        sum[i] = sum[i - 1] + (ll)i * i * mu[i] % mod;
        sum[i] %= mod;
    }

}

ll s(int n){
    return (ll)n * (n + 1) / 2 % mod;
}



ll cal(int a, int b){
    ll res = 0;
    if(a > b) swap(a, b);
    for(int l = 1, r; l <= a; l = r + 1){
        r = min(a / (a / l), b / (b / l));
        ll ans = 0;
        int ad = a / l, bd = b / l;
        for(int i = 1, j; i <= min(ad, bd); i = j + 1){
            j = min(ad / (ad / i), bd / (bd / i));
            ans += (ll)(sum[j] - sum[i - 1]) * s(ad / i) % mod * s(bd / i) % mod;
            ans %= mod;
        }
        res += (ll)(s(r) - s(l - 1)) * ans % mod;
        res %= mod;
    }
    return (res + mod) % mod;
}


void solve()
{
    int n, m;
    sdd(n, m);
    printf("%lld\n", cal(n, m));

}

int main()
{
#ifdef ONLINE_JUDGE
#else
    freopen("/home/jungu/code/in.txt", "r", stdin);
    // freopen("/home/jungu/code/out.txt", "w", stdout);
    // freopen("/home/jungu/code/practice/out.txt","w",stdout);
#endif
    // ios::sync_with_stdio(false);
    cin.tie(0), cout.tie(0);

    int T = 1;
    // sd(T);
    k = 10000000;
    get(k);
    // int cas = 1;
    while (T--)
    {
        // printf("Case #%d:", cas++);
        solve();
    }

    return 0;
}