洛谷 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)}dfrac{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} dsum_{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} dsum_{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
",ans);
	return 0;
}