P1829 [国家集训队]Crash的数字表格（莫比乌斯反演+数论分块）

题意：

$给定n,m,要求计算{\sum }_{i=1}^n{\sum }_{j=1}^{m}lcm(i,j)$

$答案对20101009取模$

数据范围：n,m<=1e7

解法：

$$\begin{gathered} ans=\sum _{i=1}^n\sum _{j=1}^{m}lcm(i,j) \\ =\sum _{i=1}^n\sum _{j=1}^m\frac{i\ast j}{gcd(i,j)} \\ =\sum _{d=1}^{min(n,m)}\sum _{i=1}^n\sum _{j=1}^m\frac{i\ast j}{d}[gcd(i,j)=d] \\ =\sum _{d=1}^{min(n,m)}\sum _{i=1}^n\sum _{j=1}^m\frac{i}{d}\ast \frac{j}{d}\ast d\ast [gcd(i,j)=d] \\ =\sum _{d=1}^{min(n,m)}d\sum _{i=1}^n\sum _{j=1}^m\frac{i}{d}\ast \frac{j}{d}\ast [gcd(i,j)=d] \\ =\sum _{d=1}^{min(n,m)}d\sum _{i=1}^{n/d}\sum _{j=1}^{m/d}i\ast j\ast [gcd(i,j)=1] \\ 设f(n,m)=\sum _{i=1}^n\sum _{j=1}^{m}i\ast j\ast [gcd(i,j)=1] \\ 存在反演结论:[gcd(i,j)=1]=\sum _{d|gcd(i,j)}\mu (d) \\ 因此f(n,m)=\sum _{i=1}^n\sum _{j=1}^{m}i\ast j\sum _{d|gcd(i,j)}\mu (d) \\ =\sum _{d=1}^{min(n,m)}\mu (d)\sum _{d|i,i<=n}\sum _{d|j,j<=m}i\ast j \\ =\sum _{d=1}^{min(n,m)}\mu (d)\sum _{i=1}^{n/d}\sum _{j=1}^{m/d}(i\ast d)\ast (j\ast d) \\ =\sum _{d=1}^{min(n,m)}\mu (d)d^2\sum _{i=1}^{n/d}\sum _{j=1}^{m/d}i\ast j \\ 设s(i)=\frac{i\ast (i+1)}{2} \\ 那么f(n,m)=\sum _{d=1}^{min(n,m)}\mu (d)d^{2}s(\frac{n}{d})s(\frac{m}{d}) \\ f(n,m)可以数论分块优化到O(\sqrt{n}) \\ 继续之前对ans的推导: \\ ans=\sum _{d=1}^{min(n,m)}d\sum _{i=1}^{n/d}\sum _{j=1}^{m/d}i\ast j\ast [gcd(i,j)=1] \\ =\sum _{d=1}^{min(n,m)}d\ast f(\frac{n}{d},\frac{m}{d}) \\ 可以用数论分块优化到O(\sqrt{n}) \\ 由于里面还套了一个f(n,m)的O(\sqrt{n}) \\ 所以总复杂度为O(\sqrt{n}\ast \sqrt{n})=O(n) \end{gathered}$$

code：

#include<bits/stdc++.h>
#define int long long
using namespace std;
const int maxm=1e7+5;
const int mod=20101009;
int p[maxm],c;
int np[maxm];
int mu[maxm];
int sum[maxm];
int n,m;
void init(){
    mu[1]=1;
    for(int i=2;i<maxm;i++){
        if(!np[i])p[c++]=i,mu[i]=-1;
        for(int j=0;j<c;j++){
            if(p[j]*i>=maxm)break;
            np[p[j]*i]=1;
            mu[p[j]*i]=(i%p[j]?-mu[i]:0);
            if(i%p[j]==0)break;
        }
    }
    for(int i=1;i<maxm;i++){
        sum[i]=(sum[i-1]+mu[i]*i*i%mod)%mod;
    }
}
int s(int x){
    return x*(x+1)/2%mod;
}
int f(int n,int m){
    int mi=min(n,m);
    int ans=0;
    for(int l=1,r;l<=mi;l=r+1){
        r=min(n/(n/l),m/(m/l));
        ans+=s(n/l)*s(m/l)%mod*(sum[r]-sum[l-1])%mod;
        ans%=mod;
    }
    return ans;
}
signed main(){
    ios::sync_with_stdio(0);
    init();
    cin>>n>>m;
    int mi=min(n,m);
    int ans=0;
    for(int l=1,r;l<=mi;l=r+1){
        r=min(n/(n/l),m/(m/l));
        ans+=(s(r)-s(l-1))%mod*f(n/l,m/l)%mod;
        ans%=mod;
    }
    ans=(ans%mod+mod)%mod;
    cout<<ans<<endl;
    return 0;
}

---