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

gate

\(\sum\limits_{i = 1}^{n}\sum\limits_{j = 1}^{m} lcm(i,j)\)

\(= \sum\limits_{i = 1}^{n}\sum\limits_{j = 1}^{m} \dfrac{i\cdot j}{gcd(i,j)}\)

\(= \sum\limits_{i = 1}^{n}\sum\limits_{j = 1}^{m} \sum\limits_{d|i,d|j}\dfrac{i\cdot j}{d}\ [gcd(\dfrac{i}{d},\dfrac{j}{d})=1]\)

\(= \sum\limits_{d=1}^{n}d\sum\limits_{i = 1}^{\frac{n}{d}}\sum\limits_{j = 1}^{\frac{m}{d}} i\cdot j\ [gcd(i,j)=1]\)

设\(f(n,m)=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}i\cdot j[gcd(i,j)=1]\)

\(=\sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\sum\limits_{d|i,d|j}\mu(d)\cdot i\cdot j\)

\(=\sum\limits_{d=1}^{n}\mu(d)\cdot d^2\sum\limits_{i=1}^{\frac{n}{d}}\sum\limits_{j=1}^{\frac{m}{d}}\ i\cdot j\)

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

设\(g(n,m) = \sum\limits_{i=1}^{n}\sum\limits_{j=1}^{m}\ i\cdot j\)

\(= \sum\limits_{i=1}^{n}i\sum\limits_{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\n",solve(n,m));
    return 0;
}