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

题目
题目大意就是求∑ni=1∑mj=1lcm(i,j).
分析：

ans=∑i=1n∑j=1mlcm(i,j)=∑i=1n∑j=1mi∗jgcd(i,j)

枚举d=gcd(i,j),

令f(n,m,k)=∑1<=i<=n1<=j<=mgcd(i,j)=ki∗j

显然可以得到

ans=∑d=1min(n,,m)d2∗f(⌊nd⌋,⌊md⌋,1)d=∑d=1min(n,,m)d∗f(⌊nd⌋,⌊md⌋,1)

即算出除以d后互质的对数，两个数都乘d的乘积就得到了两个数的乘积，在除以d就是他们的最小公倍数。
那怎么求f呢
令

sum(x,y)=∑1<=i<=x1<=j<=yi∗j=x∗(x+1)2∗y∗(y+1)2

F(x,y,k)=∑1<=i<=x1<=j<=yk|gcd(i,j)i∗j=k2∗sum(⌊nk⌋⌊mk⌋)=∑k|df(i,j,d)

f(x,y,k)=∑k|dμ(dk)F(x,y,d)

所以

f(x,y,1)=∑k=1min(m,n)μ(k)F(x,y,k)=∑k=1min(x,y)μ(k)∗k2∗sum(⌊xk⌋⌊yk⌋)=∑k=1min(x,y)μ(k)∗k2∗⌊xk⌋∗(⌊xk⌋+1)2∗⌊yk⌋∗(⌊yk⌋+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);}