bzoj2154 Crash的数字表格

题目链接：bzoj2154
题目大意：
给定n、m，求∑ni=1∑mj=1lcm(i∗j)。
1≤n，m≤10000000

题解：
莫比乌斯反演
把原式化开

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

枚举约数(省略前面的了)

∑i=1n∑j=1m∑d|i,d|ji∗jd[gcd(i,j)=d]

交换枚举倍数和约数

∑d=1min(n,m)∑i=1⌊nd⌋∑j=1⌊md⌋id∗jdd[gcd(i,j)=1]

[gcd(i,j)=1]就是ϵ((i,j))，套1∗μ=ϵ再约下分什么的就变成了

∑d=1min(n,m)d∑i=1⌊nd⌋∑j=1⌊md⌋i∗j∑t|i,t|jμ(t)

再次交换枚举倍数和约数

∑d=1min(n,m)d∑t=1min(⌊nd⌋,⌊md⌋)∑i=1⌊ndt⌋∑j=1⌊mdt⌋it∗jt∗μ(t)

∑d=1min(n,m)d∑t=1min(⌊nd⌋,⌊md⌋)t2∗μ(t)∑i=1⌊ndt⌋∑j=1⌊mdt⌋i∗j

设sum(n,m)=∑ni=1∑mj=1i∗j=∑ni=1i∑mj=1j这个直接上求和公式
则上式变成

∑d=1min(n,m)d∑t=1min(⌊nd⌋,⌊md⌋)t2∗μ(t)∗sum(⌊ndt⌋,⌊mdt⌋)

设f(n,m)=∑min(n,m)t=1t2∗μ(t)∗sum(⌊nt⌋,⌊mt⌋)这个分块求。
于是上式又变成了

∑d=1min(n,m)d×f(⌊nd⌋,⌊md⌋)

这个也可以分块了。
所以时间复杂度大概是O(n)..吧
LL、mod什么的真的要注意啊不然WA到死。

#include<cstdio>#include<cstdlib>#include<cstring>#include<iostream>#include<algorithm>using namespace std;typedef long long LL;#define N 10000100const LL ny=10050505;const LL mod=20101009;int cnt,pri[N/10];bool ispri[N];LL mu[N],pre[N];int mymin(int x,int y){return (x<y)?x:y;}void mobius(int lim){    cnt=0;mu[1]=1;    for (int i=2;i<=lim;i++)    {        if (!ispri[i]) {pri[++cnt]=i;mu[i]=-1;}        for (int j=1;j<=cnt && i*pri[j]<=lim;j++)        {            ispri[i*pri[j]]=true;            if (i%pri[j]==0)            {                mu[i*pri[j]]=0;                break;            }            mu[i*pri[j]]=-mu[i];        }    }}LL sum(int n,int m){    LL s1=(LL)(n)*(LL)(n+1)%mod*ny%mod;    LL s2=(LL)(m)*(LL)(m+1)%mod*ny%mod;    return s1*s2%mod;}LL f(int x,int y){    int r,lim=mymin(x,y);LL ans=0;    for (int i=1;i<=lim;i=r+1)    {        r=mymin(x/(x/i),y/(y/i));        ans=(ans+sum(x/i,y/i)*(pre[r]-pre[i-1])%mod)%mod;    }while (ans<0) ans+=mod;    return ans;}int main(){    //freopen("a.in","r",stdin);    //freopen("a.out","w",stdout);    int n,m,r,i,lim;LL ans=0;    scanf("%d%d",&n,&m);    lim=mymin(n,m);mobius(lim);    pre[1]=mu[1];for (i=2;i<=lim;i++) pre[i]=(pre[i-1]+(LL)i*i%mod*mu[i]%mod)%mod;    for (i=1;i<=lim;i=r+1)    {        r=mymin(n/(n/i),m/(m/i));        ans=(ans+(LL)(r-i+1)*(r+i)%mod*ny%mod*f(n/i,m/i)%mod)%mod;    }while (ans<0) ans+=mod;    printf("%lld\n",ans);    return 0;}