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

题目链接

两次数论分块。
式子：

\[\sum_{i=1}^{n}\sum_{j=1}^{m} lcm(i,j)\\ = \sum_{i=1}^{n}\sum_{j=1}^{m} \frac{ij}{gcd(i,j)}\\ = \sum_{i=1}^{n}\sum_{j=1}^{m} \sum_{d=gcd(i,j)}\frac{ij}{d}\\ = \sum_{i=1}^{n}\sum_{j=1}^{m} \sum_{d=gcd(i,j)}d\frac{i}{d}\frac{j}{d}\\ = \sum_{d=1}^{n} d \sum_{i=1}^{n}\sum_{j=1}^{m} \frac id \frac jd[gcd(i,j) = d]\\ = \sum_{d=1}^{n} d\sum_{i=1}^{\frac nd} \sum_{j=1}^{\frac md} ij[gcd(i,j) =1]\\ *\\ *\\ f(n,m) = \sum_{i=1}^{n} \sum_{j=1}^{m}ij[gcd(i,j) = 1]\\ = \sum_{i=1}^{n} \sum_{j=1}^{m} ij \sum_{p | gcd(i,j) } \mu(p)\\ = \sum_{i=1}^{n} \sum_{j=1}^{m} ij \sum_{p | i, p | j } \mu(p)\\ =\sum_{p=1}^{n} \mu(p) p^2 \sum_{i=1}^{\frac np}\sum_{j=1}^{\frac mp} ij\\ =\sum_{p=1}^{n} \mu(p) p^2 g(\frac np) g(\frac mp)\\ g(n) = \frac{n(n + 1)}{2}\\ *\\ *\\ \sum_{d=1}^{n} d\sum_{i=1}^{\frac nd} \sum_{j=1}^{\frac md} ij[gcd(i,j) =1]\\ = \sum_{d=1}^{n} d f(\frac nd , \frac md)\\ \]

#include<bits/stdc++.h>
using namespace std;

#define ll long long 
const int N = 1e7;
const int mod = 20101009;

int p[N + 5], prime[N / 10], cnt;
int mu[N + 5], sum[N + 5];

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

ll f(ll n, ll m){
    ll ret = 0;
    for(int l = 1,r; l <= n; l = r + 1){
        r = min(min(n / (n/l), m / (m/l)), n);
        ll v = 1ll * g(n/l) * g(m/l) % mod;
        ret = (ret + 1ll * ((sum[r] - sum[l - 1] + mod) % mod) * v % mod) % mod;  
    }
    return ret;
}

int n,m;

int main(){
    p[0] = p[1] = 1; mu[1] = 1;
    for(int i = 2; i <= N; ++ i){
        if(!p[i]) { prime[++ cnt] = i; mu[i] = -1; }
        for(int j = 1; j <= cnt && 1ll * prime[j] * i <= N; ++ j){
            p[prime[j] * i] = 1;
            if(i % prime[j] == 0) { mu[prime[j] * i] = 0; break; }
            else mu[prime[j] * i] = - mu[i];
        }
    }
    for(int i = 1; i <= N; ++ i) sum[i] = (sum[i - 1] + (1ll * mu[i] * (1ll * i * i % mod) % mod) % mod + mod) % mod;
    
    scanf("%d%d",&n,&m);
    if(n > m) swap(n,m);
    ll ans = 0;
    for(int l = 1, r; l <= n; l = r + 1){
        r = min(min(n / (n/l), m / (m/l)), n);
        ll v = f(n/l, m/l);
        ans = (ans + 1ll * v * ((1ll * (r - l + 1) * (l + r) / 2) % mod) % mod ) % mod; 
    }
    printf("%lld\n",ans);
	return 0;
}