#莫比乌斯反演，整除分块#bzoj 2154 bzoj 2693 jzoj 1938 洛谷 1829 Crash的数字表格 or 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 _{d=1}^{min(n,m)}d\times \sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{m}{d}\rfloor }[gcd(i,j)=1]\times i\times j$$
根据莫比乌斯反演可以得到原式=$$\sum _{d=1}^{min(n,m)}d\sum _{t=1}^{min(\lfloor \frac{n}{d}\rfloor ,\lfloor \frac{m}{d}\rfloor )}{t}^2\mu (t)\sum _{i=1}^{\lfloor \frac{n}{dt}\rfloor }i\sum _{j=1}^{\lfloor \frac{m}{dt}\rfloor }j$$
就可以预处理$t^2\mu (t)$,用整除分块求解，时间复杂度大概是$O(n)$

---

代码

#include <cstdio>
#define rr register
#define min(a,b) ((a)<(b)?(a):(b))
#define res(n,m) ((1ll*n*(n+1)>>1)%mod*((1ll*m*(m+1)>>1)%mod)%mod)
using namespace std;
const int mod=20101009,M=10000010;
int n,m,v[M],prime[M],s[M],cnt;
inline void prepa(int N){
    s[1]=1;
    for (rr int i=2;i<=N;++i){
        if (!v[i]) s[i]=-1,v[i]=prime[++cnt]=i;
        for (rr int j=1;j<=cnt&&prime[j]*i<=N;++j){
            v[i*prime[j]]=prime[j];
            if (i%prime[j]) s[i*prime[j]]=-s[i];
                else break;
        }
    }
    for (rr int i=1;i<=N;++i) s[i]=1ll*i*i%mod*(s[i]+mod)%mod;
    for (rr int i=2;i<=N;++i) (s[i]+=s[i-1])%=mod;
}
inline signed answ(int n,int m){
    rr int ans=0,t=min(n,m);
    for (rr int l=1,r;l<=t;l=r+1){
        r=min(n/(n/l),m/(m/l));
        (ans+=1ll*(s[r]-s[l-1]+mod)*res(n/l,m/l)%mod)%=mod;
    }
    return ans;
}
inline signed solve(int n,int m){
    rr int ans=0,t;
    prepa(t=min(n,m)); 
    for (rr int l=1,r;l<=t;l=r+1){
        r=min(n/(n/l),m/(m/l));
        (ans+=(1ll*(r-l+1)*(l+r)>>1)%mod*answ(n/l,m/l)%mod)%=mod;
    }
    return ans;
}
signed main(){
    scanf("%d%d",&n,&m);
    printf("%d",solve(n,m));
    return 0;
}