【莫比乌斯反演】[HYSBZ\BZOJ2154]Crash的数字表格

题目
题目大意就是求$\sum_{i=1}^{n}\sum^{m}_{j=1}lcm(i,j)$.
分析：

$$ans=\sum_{i=1}^{n}\sum^{m}_{j=1}lcm(i,j)=\sum_{i=1}^{n}\sum^{m}_{j=1}\frac{i*j}{gcd(i,j)}$$
枚举d=gcd(i,j), $$令f(n,m,k)=\sum_{1<=i<=n}_{1<=j<=m}_{gcd(i,j)=k}i*j$$
显然可以得到
$$ans=\sum_{d=1}^{min(n,,m)}\frac{d^2*f(\lfloor\frac nd\rfloor,\lfloor\frac md\rfloor,1)}{d}=\sum_{d=1}^{min(n,,m)}d*f(\lfloor\frac nd\rfloor,\lfloor\frac md\rfloor,1)$$
即算出除以d后互质的对数，两个数都乘d的乘积就得到了两个数的乘积，在除以d就是他们的最小公倍数。
那怎么求f呢
令 $$sum(x,y)=\sum_{1<=i<=x}_{1<=j<=y}i*j=\frac{x*(x+1)}{2}*\frac{y*(y+1)}{2}$$
$$F(x,y,k)=\sum_{1<=i<=x}_{1<=j<=y}_{k|gcd(i,j)}i*j=k^2*sum(\lfloor\frac nk\rfloor\lfloor\frac mk\rfloor)=\sum_{k|d}f(i,j,d)$$
$$f(x,y,k)=\sum_{k|d}\mu(\frac dk)F(x,y,d)$$
所以 $$f(x,y,1)=\sum_{k=1}^{min(m,n)}\mu(k)F(x,y,k)=\sum_{k=1}^{min(x,y)}\mu(k)*k^2*sum(\lfloor\frac xk\rfloor\lfloor\frac yk\rfloor)\\=\sum_{k=1}^{min(x,y)}\mu(k)*k^2*\frac{\lfloor\frac xk\rfloor*(\lfloor\frac xk\rfloor+1)}{2}*\frac{\lfloor\frac yk\rfloor*(\lfloor\frac yk\rfloor+1)}{2}$$
然后运用分块优化求出ans即可（可参考 [HYSBZ/BZOJ2301]Problem b中的分块优化）

#include<cstdio>
#include<algorithm>
using namespace std;
#define MAXN 10000000
#define MOD 20101009
#define Sum(x,y) (1ll*x*(x+1)/2%MOD*(1ll*y*(y+1)/2%MOD)%MOD)
int mu[MAXN+10],p[MAXN+10],pcnt,sum[MAXN+10],n,m,ans;
bool f[MAXN+10];
void prepare(){
    int i,j,t=min(n,m);
    mu[1]=sum[1]=1;
    for(i=2;i<=t;i++){
        if(!f[i])
            p[++pcnt]=i,mu[i]=-1;
        for(j=1;p[j]*i<=t;j++){
            f[i*p[j]]=1;
            if(i%p[j]==0){
                mu[i*p[j]]=0;
                break;
            }   
            mu[i*p[j]]=-mu[i];
        }
        sum[i]=(sum[i-1]+1ll*i*i*mu[i])%MOD;
    }
}
void Read(int &x){
    char c;
    while(c=getchar(),~c)
        if(c>='0'&&c<='9'){
            x=c-'0';
            while(c=getchar(),c>='0'&&c<='9')
                x=x*10+c-'0';
            ungetc(c,stdin);
            return;
        }
}
int calf(int n,int m){
    int t=min(n,m),i,last,ret=0;
    for(i=1;i<=t;i=last+1){
        last=min(n/(n/i),m/(m/i));
        ret=(ret+((sum[last]-sum[i-1])%MOD*Sum(n/i,m/i))%MOD+MOD)%MOD;
    }
    return ret;
}
void solve(){
    int t=min(n,m),i,last;
    for(i=1;i<=t;i=last+1){
        last=min(n/(n/i),m/(m/i));
        ans=(ans+1ll*(i+last)*(last-i+1)/2%MOD*calf(n/i,m/i))%MOD;
    }
}
int main()
{
    Read(n),Read(m);
    prepare();
    solve();
    printf("%d\n",ans);
}