【2011集训队出题】Crash的数字表格

Description

求

$$\sum_{i=1}^{n}\sum_{j=1}^{m}lcm(i,j)$$
n,m<=10^7

Solution

（注意，以下内容默认$n<=m$）
看着lcm不爽，把它变一变：

$$\sum_{i=1}^{n}\sum_{j=1}^{m}{ij\over \gcd(i,j)}$$
莫比乌斯反演常用，枚举gcd
$$\sum_{d=1}^{n}{{1\over d}\sum_{i=1}^{n}\sum_{j=1,\gcd(i,j)=d}^{m}}ij$$
设 $$f(d)=\sum_{i=1}^{n}\sum_{j=1,\gcd(i,j)=d}^{m}ij$$
则 $$Ans=\sum_{d=1}^{n}{f(d)\over d}$$
考虑反演常用，设 $$g(d)=\sum_{i=1}^{n}\sum_{j=1,d|\gcd(i,j)}^{m}ij$$
即 $$g(d)=\sum_{x=1}^{\lfloor {n\over d}\rfloor}f(xd)$$
可以发现， $$g(d)=\sum_{i=1}^{\lfloor {n\over d}\rfloor}\sum_{j=1}^{\lfloor {m\over d}\rfloor}d^2ij$$
$$g(d)=d^2\sum_{i=1}^{\lfloor {n\over d}\rfloor}\sum_{j=1}^{\lfloor {m\over d}\rfloor}ij$$
$$g(d)=d^2\sum_{i=1}^{\lfloor {n\over d}\rfloor}i\sum_{j=1}^{\lfloor {m\over d}\rfloor}j$$
$$g(d)={d^2{\lfloor {n\over d}\rfloor}({\lfloor {n\over d}\rfloor+1}){\lfloor {m\over d}\rfloor}({\lfloor {n\over d}\rfloor}+1)\over 4}$$
根据反演得 $$f(d)=\sum_{i=1}^{\lfloor {n\over d}\rfloor}\mu(i)g(di)$$
$$f(d)=\sum_{i=1}^{\lfloor {n\over d}\rfloor}{\mu(i)d^2i^2{\lfloor {n\over di}\rfloor}({\lfloor {n\over di}\rfloor+1}){\lfloor {m\over di}\rfloor}({\lfloor {n\over di}\rfloor}+1)\over 4}$$
$$Ans={1\over 4}\sum_{d=1}^{n}d\sum_{i=1}^{\lfloor {n\over d}\rfloor}\mu(i)i^2{\lfloor {n\over di}\rfloor}({\lfloor {n\over di}\rfloor+1}){\lfloor {m\over di}\rfloor}({\lfloor {n\over di}\rfloor}+1)$$
这样发现后面的值只和 $\lfloor {n\over di}\rfloor$和 $\lfloor {m\over di}\rfloor$有关于是就可以分块了。（分块大法好！）
然后求和再用一次分块。（分块大法好！）
这样就可以 $O(n)$解决问题了。（分块大法好！）

Code

#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
#define fo(i,a,b) for(ll i=a;i<=b;i++)
#define ll long long
#define mo 20101009
#define N 10000005
using namespace std;
int p[N],tot,mi[N],mu[N];
bool f[N];
ll n,m,ans;
ll sum(ll l,ll r) {return (r+l)*(r-l+1)/2%mo;}
ll get(ll x,ll y) {return ((x*(x+1)/2)%mo)*((y*(y+1)/2)%mo)%mo;}
int main() {
    scanf("%lld%lld",&n,&m);
    if (n>m) swap(n,m);mu[1]=1;
    fo(i,2,n) {
        if (!f[i]) p[++tot]=i,mu[i]=-1;
        fo(j,1,tot) {
            int k=p[j]*i;if (k>n) break;f[k]=1;
            if (!(i%p[j])) break;mu[k]=-mu[i];
        }
    }
    fo(i,1,n) mi[i]=(mi[i-1]+mu[i]*i%mo*i%mo)%mo;
    ll l=1,r;
    while (l<=n) {
        ll nn=n / l,mm=m / l;r=min(n/nn,m / mm);
        ll le=1,ri,an=0;
        while (le<=nn) {
            ll x=nn/le,y=mm/le;ri=min(nn/x,mm/y);
            an=(an+(mi[ri]-mi[le-1])%mo*get(x,y)%mo)%mo;le=ri+1;
        }
        ans=(ans+sum(l,r)*an%mo)%mo;l=r+1;
    }
    printf("%lld",(ans+mo)%mo);
}