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

gate

(sumlimits_{i = 1}^{n}sumlimits_{j = 1}^{m} lcm(i,j))

(= sumlimits_{i = 1}^{n}sumlimits_{j = 1}^{m} dfrac{icdot j}{gcd(i,j)})

(= sumlimits_{i = 1}^{n}sumlimits_{j = 1}^{m} sumlimits_{d|i,d|j}dfrac{icdot j}{d} [gcd(dfrac{i}{d},dfrac{j}{d})=1])

(= sumlimits_{d=1}^{n}dsumlimits_{i = 1}^{frac{n}{d}}sumlimits_{j = 1}^{frac{m}{d}} icdot j [gcd(i,j)=1])

设(f(n,m)=sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}icdot j[gcd(i,j)=1])

(=sumlimits_{i=1}^{n}sumlimits_{j=1}^{m}sumlimits_{d|i,d|j}mu(d)cdot icdot j)

(=sumlimits_{d=1}^{n}mu(d)cdot d^2sumlimits_{i=1}^{frac{n}{d}}sumlimits_{j=1}^{frac{m}{d}} icdot j)

前半部分(sumlimits_{d=1}^{n}mu(d)cdot d^2)可以用前缀和预处理；

设(g(n,m) = sumlimits_{i=1}^{n}sumlimits_{j=1}^{m} icdot j)

(= sumlimits_{i=1}^{n}isumlimits_{j=1}^{m}j)

(= dfrac{n(n+1)}{2}cdot dfrac{m(m+1)}{2})

枚举(d)为(O(sqrt n))，计算(f(x,y))为(O(sqrt n))，总体为(O(n))

code

#include<cstdio>
#include<iostream>
#include<cmath>
#include<cstring>
#define MogeKo qwq
using namespace std;

#define int long long

const int maxn = 1e7+10;
const int N = 1e7;
const int mod = 20101009;
int n,m,cnt;
int prime[maxn],sum[maxn],mu[maxn];
bool vis[maxn];

void Prime(){
    sum[1] = mu[1] = 1;
    for(int i = 2;i <= N;i++){
        if(!vis[i]){
            prime[++cnt] = i;
            mu[i] = -1;
        }
        for(int j = 1;j <= cnt && i*prime[j] <= N;j++){
            vis[i*prime[j]] = true;
            if(i % prime[j] == 0) break;
            mu[i*prime[j]] = -mu[i];
        }
        sum[i] = sum[i-1] + 1ll * i*i % mod * (mu[i]+mod)%mod;
        sum[i] %= mod;
    }
}

int g(int a,int b){
    return (a*(a+1)/2 % mod) * (b*(b+1)/2 % mod) % mod;
}

int f(int a,int b){
    int ans = 0;
    if(a > b) swap(a,b);
    for(int i = 1,r;i <= a;i = r+1){
        r = min(a/(a/i),b/(b/i));
        ans += (sum[r]-sum[i-1] + mod) % mod * g(a/i,b/i) % mod;
        ans %= mod;
    }
    return ans;
}

int solve(int a,int b){
    int ans = 0;
    if(a > b) swap(a,b);
    for(int i = 1,r;i <= a;i = r+1){
        r = min(a/(a/i),b/(b/i));
        ans += f(a/i,b/i) * ((i+r)*(r-i+1)/2 % mod) % mod;
        ans %= mod;
    }
    return ans;
}

signed main(){   
    scanf("%lld%lld",&n,&m);
    Prime();
    printf("%lld
",solve(n,m));
    return 0;
}