[BZOJ2154] Crash的数字表格

题目的意思就是求出∑ni=1∑mj=1lcm(i,j).

我们可以把问题进行转换,有:

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

然后我们设

f(n,m,k)=∑i=1n∑j=1mi∗j∗[gcd(i,j)==k]

F(n,m,k)=∑i=1n∑j=1mi∗j∗[k|gcd(i,j)]

然后继续把F(n,m,k)化简一下,有:

F(n,m,k)=∑i=1n∑j=1mi∗j∗[k|gcd(i,j)]=k2∗∑i=1[nk]∑j=1[mk]i∗j=∑k|if(n,m,i)

通过莫比乌斯反演我们可以得到:

f(n,m,k)=∑k|dμ(dk)F(n,m,d)

k=1时,有：

f(n,m,1)=∑i=1min(n,m)μ(i)F(n,m,i)=∑k=1min(n,m)μ(k)∗k2∑i=1[nk]∑j=1[mk]i∗j

现在我们去枚举gcd(i,j),来求出最后的答案:

ans=∑k=1min(n,m)k2∗f([nk],[mk],1))k=∑k=1min(n,m)k∗f([nk],[mk],1)

然后通过在主函数中进行分块加速,以及f函数中的分块加速,使得复杂度为O(n√)∗O(n√)=O(n).

#include<cstdio>#include<cstring>#include<iostream>#include<algorithm>#define M 10000005#define P 20101009using namespace std;int f[M],mu[M],sum[M],sz=0,n,m;bool vis[M];void Init(){    mu[1]=1;    for(int i=2;i<=n;i++){        if(!vis[i])f[sz++]=i,mu[i]=-1;        for(int j=0,ll=f[j]*i;j<sz&&ll<=n;j++,ll=f[j]*i){            vis[ll]=1;            if(i%f[j]==0){                mu[ll]=0;                break;            }            mu[ll]=-mu[i];        }    }    for(int i=1;i<=n;i++)sum[i]=(sum[i-1]+1LL*i*i*mu[i])%P;}int Sum(int x,int y){    return (1LL*x*(x+1)/2%P)*(1LL*y*(y+1)/2%P)%P;}void Add(int &x,int y){    if(y<0)y+=P;    x+=y;    if(x>=P)x-=P;}int w(int x,int y){    int re=0;    for(int i=1;i<=x;i++){        int l=x/i,r=y/i;        l=x/l,r=y/r;        if(l>r)l=r;        Add(re,1LL*(sum[l]-sum[i-1])*Sum(x/i,y/i)%P);        i=l;    }    return re;}int solve(){    int re=0;    for(int i=1;i<=n;i++){        int l=n/i,r=m/i;        l=n/l,r=m/r;        if(l>r)l=r;        Add(re,1LL*w(n/i,m/i)*(Sum(l,1)-Sum(i-1,1))%P);        i=l;    }    return re;}int main(){    scanf("%d %d",&n,&m);    if(n>m)swap(n,m);    Init();    printf("%d\n",solve());    return 0;}