[BZOJ2154]Crash的数字表格

题目大意：
　　给定$n,m(n,m\leq10^7)$，求$\displaystyle\sum_{x=1}^n\sum_{y=1}^m{\rm lcm}(x,y)$。

思路：
　　令$d=\gcd(x,y),x=ad,y=bd$，则：
$$
\begin{align*}
原式&=\sum_{d=1}^{\min(n,m)}d\sum_{a=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{b=1}^{\lfloor\frac{m}{d}\rfloor}ab[\gcd(ab)=1]\\
&=\sum_{d=1}^{\min(n,m)}d\sum_{a=1}^{\lfloor\frac{n}{d}\rfloor}\sum_{b=1}^{\lfloor\frac{m}{d}\rfloor}ab\sum_{p|\gcd(a,b)}\mu(p)
\end{align*}
$$
　　令$p=\gcd(a,b),a=jp,b=kp$，则：
$$
\begin{align*}
原式&=\sum_{d=1}^{\min(n,m)}d\sum_{p=1}^{\min(\lfloor\frac{n}{d}\rfloor,\lfloor\frac{m}{d}\rfloor)}\mu(p)p^2\sum_{j=1}^{\lfloor\frac{n}{dp}\rfloor}\sum_{k=1}^{\lfloor\frac{m}{dp}\rfloor}jk\\
&=\sum_{d=1}^{\min(n,m)}d\sum_{p=1}^{\min(\lfloor\frac{n}{d}\rfloor,\lfloor\frac{m}{d}\rfloor)}\mu(p)p^2\frac{(\lfloor\frac{n}{dp}\rfloor+1)\lfloor\frac{n}{dp}\rfloor}{2}\cdot\frac{(\lfloor\frac{m}{dp}\rfloor+1)\lfloor\frac{m}{dp}\rfloor}{2}\\
\end{align*}
$$
　　预处理$\mu(p)p^2$的前缀和，暴力枚举$d$。对于$\lfloor\frac{n}{dp}\rfloor$和$\lfloor\frac{m}{dp}\rfloor$分别相等的$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\n",ans);
    return 0;
}