P1829 JZPTAB(莫反)

P1829 JZPTAB(莫反)

∑ i = 1 n ∑ j = 1 m l c m ( i , j ) \large \sum\limits_{i=1}^n\sum\limits_{j=1}^m lcm(i,j) i=1∑nj=1∑mlcm(i,j)

= ∑ d = 1 n d ∑ k = 1 n d μ ( k ) k 2 g ( n k d , m k d ) \large=\sum\limits_{d=1}^n d\sum\limits_{k=1}^{\small\dfrac{n}{d}}\mu(k)k^2g(\dfrac{n}{kd},\dfrac{m}{kd}) =d=1∑ndk=1∑dnμ(k)k2g(kdn,kdm)

g ( n , m ) = ∑ i = 1 n i ∑ j = 1 m j = n ( n + 1 ) 2 × m ( m + 1 ) 2 \large g(n,m)=\sum\limits_{i=1}^n i\sum\limits_{j=1}^mj=\dfrac{n(n+1)}{2}\times \dfrac{m(m+1)}{2} g(n,m)=i=1∑nij=1∑mj=2n(n+1)×2m(m+1)

预处理 ∑ μ ( i ) i 2 \sum \mu(i)i^2 ∑μ(i)i2 即可，然后两次数论分块即可。

#include<bits/stdc++.h>
using namespace std;
typedef long long ll;
const int N=1e7+5,M=2e4+5,inf=0x3f3f3f3f,mod=20101009;
#define mst(a,b) memset(a,b,sizeof a)
#define PII pair<int,int>
#define fi first
#define se second
#define pb push_back
int p[N],cnt,mu[N],pre[N];
bitset<N>vis;
void getmu(int n){
    mu[1]=vis[1]=1;
    int cnt=0;
    for(int i=2;i<=n;i++){
        if(!vis[i]) p[++cnt]=i,mu[i]=-1;
        for(int j=1;j<=cnt&&i*p[j]<=n;j++){
            vis[i*p[j]]=1;
            if(i%p[j]==0){
                mu[i*p[j]]=0;
                break;
            }
            mu[i*p[j]]=-mu[i];
        }
    }
    for(int i=1;i<=n;i++) pre[i]=(pre[i-1]+1LL*mu[i]*i%mod*i%mod)%mod;
}
ll g(ll n,ll m){
    return ((n*(n+1)>>1)%mod)*((m*(m+1)>>1)%mod)%mod;
}
/*ll fun(ll n){
    return (n*(n+1)%mod*(2*n+1)%mod)/6%mod;
}
ll ff(ll l,ll r){
    return (fun(r)-fun(l-1)+mod)%mod;
}*/
ll solve(int n,int m){
     if(n>m) swap(n,m);
    ll s=0;
    for(int l=1,r;l<=n;l=r+1){
       r=min(n/(n/l),m/(m/l));
       s+=1LL*(pre[r]-pre[l-1])%mod*g(n/l,m/l)%mod;
       s%=mod;
    }
    return s;
 }
ll gg(ll n){
    return n*(n+1)/2%mod;
}
ll gfun(ll l,ll r){
    return (gg(r)-gg(l-1));
}
ll fun(ll n,ll m){
    if(n>m) swap(n,m);
    ll s=0;
    for(int l=1,r;l<=n;l=r+1){
        r=min(n/(n/l),m/(m/l));
        s+=1LL*gfun(l,r)*solve(n/l,m/l)%mod;
    }
    return (s%mod+mod)%mod;
}
int main(){
    getmu(N-5);
    int n,m;
    scanf("%d%d",&n,&m);
    printf("%lld\n",fun(n,m));
    return 0;
}