[BZOJ2154]Crash的数字表格

题目大意：
　　给定$n,m(n,mleq10^7)$，求$displaystylesum_{x=1}^nsum_{y=1}^m{ m lcm}(x,y)$。

思路：
　　令$d=gcd(x,y),x=ad,y=bd$，则：
$$
egin{align*}
原式&=sum_{d=1}^{min(n,m)}dsum_{a=1}^{lfloorfrac{n}{d} floor}sum_{b=1}^{lfloorfrac{m}{d} floor}ab[gcd(ab)=1]\
&=sum_{d=1}^{min(n,m)}dsum_{a=1}^{lfloorfrac{n}{d} floor}sum_{b=1}^{lfloorfrac{m}{d} floor}absum_{p|gcd(a,b)}mu(p)
end{align*}
$$
　　令$p=gcd(a,b),a=jp,b=kp$，则：
$$
egin{align*}
原式&=sum_{d=1}^{min(n,m)}dsum_{p=1}^{min(lfloorfrac{n}{d} floor,lfloorfrac{m}{d} floor)}mu(p)p^2sum_{j=1}^{lfloorfrac{n}{dp} floor}sum_{k=1}^{lfloorfrac{m}{dp} floor}jk\
&=sum_{d=1}^{min(n,m)}dsum_{p=1}^{min(lfloorfrac{n}{d} floor,lfloorfrac{m}{d} floor)}mu(p)p^2frac{(lfloorfrac{n}{dp} floor+1)lfloorfrac{n}{dp} floor}{2}cdotfrac{(lfloorfrac{m}{dp} floor+1)lfloorfrac{m}{dp} floor}{2}\
end{align*}
$$
　　预处理$mu(p)p^2$的前缀和，暴力枚举$d$。对于$lfloorfrac{n}{dp} floor$和$lfloorfrac{m}{dp} floor$分别相等的$p$可以分块计算。

#include<cstdio>
#include<cctype>
#include<algorithm>
typedef long long int64;
inline int getint() {
    register char ch;
    while(!isdigit(ch=getchar()));
    register int x=ch^'0';
    while(isdigit(ch=getchar())) x=(((x<<2)+x)<<1)+(ch^'0');
    return x;
}
const int N=10000001,M=664580,mod=20101009;
bool vis[N];
int mu[N],p[M],sum[N];
inline void sieve(const int lim) {
    mu[1]=1;
    for(register int i=2;i<=lim;i++) {
        if(!vis[i]) {
            p[++p[0]]=i;
            mu[i]=-1;
        }
        for(register int j=1;j<=p[0]&&i*p[j]<=lim;j++) {
            vis[i*p[j]]=true;
            if(i%p[j]==0) {
                mu[i*p[j]]=0;
                break;
            }
            mu[i*p[j]]=-mu[i];
        }
    }
    for(register int i=1;i<=lim;i++) {
        sum[i]=(sum[i-1]+(int64)mu[i]*i%mod*i%mod)%mod;
    }
}
int main() {
    const int n=getint(),m=getint(),lim=std::min(n,m);
    sieve(lim);
    int ans=0;
    for(register int d=1;d<=lim;d++) {
        int tmp=0;
        const int x=n/d,y=m/d,lim=std::min(x,y);
        for(register int i=1,j;i<=lim;i=j+1) {
            j=std::min(x/(x/i),y/(y/i));
            (tmp+=(sum[j]-sum[i-1]+mod)%mod*((int64)(x/i+1)*(x/i)/2%mod)%mod*((int64)(y/i+1)*(y/i)/2%mod)%mod)%=mod;
        }
        (ans+=((int64)tmp*d)%mod)%=mod;
    }
    printf("%d
",ans);
    return 0;
}