51nod 1238 最小公倍数之和 V3

原题链接：
http://www.51nod.com/onlineJudge/questionCode.html#problemId=1238&noticeId=338278

题面错误。。。

题目的实际意思是：

G=∑i=1n∑j=1nlcm(i,j)

因为题面的错误 。反反复复推了好久。

按照一般套路：

G=∑i=1n∑j=1nlcm(i,j)=∑d=1n∑i=1n∑j=1nijd[gcd(i,j)=d]

这是一个套路吧。。（看到lcm就转gcd找容斥关系）

G=∑d=1n1d∑i=1n∑j=1nij[gcd(i,j)=d]

对于：

∑i=1n∑j=1nij[gcd(i,j)=d]

记

S1(n)=∑i=1niS2(n)=∑i=1ni2

反演有：

∑i=1n∑j=1nij[gcd(i,j)=d]=∑d|aμ(ad)a2S1(⌊na⌋)

所以：

G=∑d=1n1d∑d|aμ(ad)a2S1(⌊na⌋)=∑a=1nS1(⌊na⌋)a2∑d|aμ(ad)1d=∑a=1nS1(⌊na⌋)a∑d|aμ(d)d

令

h(n)=μ(n)n

令

f(n)=∑d|nh(d)d

显然：

h∗I=fh∗id1=e

卷积运算有：

f∗id1∗μ=h∗I∗id1∗μ=e

φ∗I=id1

所以：

φ=id1∗μ

所以：

f=h∗I=φ−

一个思路就是快速计算

af(a)

的前缀和函数。

令

F(n)=nf(n)

显然F逆函数F−

F−(n)=nφ(n)

对于：

∑i=1nF−(i)=∑i=1niφ(i)=∑i=1ni∑d|iμ(id)d=∑d=1nd∑d|iμ(id)i∑d=1nd2∑d|ih(id)

用Sh表示h的前缀和函数。其他函数一样。

则：

∑i=1nF−(i)=∑d=1nd2Sh(⌊nd⌋)

前面的分析有：

h∗id1=e

从h开始,反向计算需要的函数前缀和即可。

#include <algorithm>#include <string.h>#include <cmath>#include <stdio.h>#define MAXN 3111111using namespace std;typedef long long LL;const int P=1e9+7;int Iv[10]={0,1};int mu[MAXN]={0,-1};int phi[MAXN];int h[MAXN];int g[MAXN];int f[MAXN];int Sg[MAXN];int _Sg[MAXN];int Sh[MAXN];int _Sh[MAXN];int Sf[MAXN];int _Sf[MAXN];void init (){    for(int i=2;i<10;i++)Iv[i]=(P-(LL)(P/i)*Iv[P%i]%P)%P;    for(int i=1;i<MAXN;i++)    {        mu[i]=-mu[i];        phi[i]=i-phi[i];        for(int j=i<<1;j<MAXN;j+=i)        {            mu[j]+=mu[i];            phi[j]+=phi[i];        }        h[i]=mu[i]*i;        if(h[i]<0)h[i]+=P;        Sh[i]=Sh[i-1]+h[i];        if(Sh[i]>=P)Sh[i]-=P;        g[i]=(LL)i*phi[i]%P;        Sg[i]=Sg[i-1]+g[i];        if(Sg[i]>=P)Sg[i]-=P;    }    for(int i=1;i<MAXN;i++)    {        for(int j=i;j<MAXN;j+=i)        {            f[j]+=(P+mu[i]*i)%P;            if(f[j]>=P)f[j]-=P;        }        f[i]=(LL)f[i]*i%P;        Sf[i]=Sf[i-1]+f[i];        if(Sf[i]>=P)Sf[i]-=P;    }}int S1(LL n){    n%=P;    return n*(n+1)%P*Iv[2]%P;}int S2(LL n){    n%=P;    return n*(n+1)%P*((n<<1)+1)%P*Iv[6]%P;}void clat_h(LL n,int d){    int m=sqrt(n)+1,ans=0;    for(int L=1;L<m;L++)    {        int u=S1(n/L)-S1(n/(L+1));        u=(u+P)%P;        ans+=(LL)Sh[L]*u%P;        if(ans>=P)ans-=P;    }    for(int i=(int)(n/m);i>1;i--)    {        LL u=n/i;        if(u<MAXN)            u=Sh[u];        else            u=_Sh[i*d];        ans+=u*i%P;        if(ans>=P)ans-=P;    }    _Sh[d]=(1-ans+P)%P;}void slove_h(LL n){    if(n<MAXN)return;    for(int i=(int)(n/MAXN);i;i--)clat_h(n/i,i);}void clat_g(LL n,int d){    int m=sqrt(n)+1,ans=0;    for(int L=1;L<m;L++)    {        int u=S2(n/L)-S2(n/(L+1));        u=(u+P)%P;        ans+=(LL)Sh[L]*u%P;        if(ans>=P)ans-=P;    }    for(int i=(int)(n/m);i;i--)    {        int bu=(LL)i*i%P;        LL u=n/i;        if(u<MAXN)            u=Sh[u];        else            u=_Sh[i*d];        ans+=u*bu%P;        if(ans>=P)ans-=P;    }    _Sg[d]=ans;}void slove_g(LL n){    if(n<MAXN)return;    for(int i=(int)(n/MAXN);i;i--)clat_g(n/i,i);}void clat_f(LL n,int d){    int m=sqrt(n)+1,ans=0;    for(int L=1;L<m;L++)    {        LL f1=n/L,f2=n/(L+1);        if(f1<MAXN)            f1=Sg[f1];        else            f1=_Sg[L*d];        if(f2<MAXN)            f2=Sg[f2];        else            f2=_Sg[(L+1)*d];        ans+=(LL)Sf[L]*(f1-f2+P)%P;        if(ans>=P)ans-=P;    }    for(int i=(int)(n/m);i>1;i--)    {        LL u=n/i;        if(u<MAXN)            u=Sf[u];        else            u=_Sf[i*d];        ans+=u*g[i]%P;        if(ans>=P)ans-=P;    }    _Sf[d]=(1-ans+P)%P;}void slove_f(LL n){    if(n<MAXN)return;    for(int i=(int)(n/MAXN);i;i--)clat_f(n/i,i);}int main (){    init();    LL n;    scanf("%lld",&n);    slove_h(n);    slove_g(n);    slove_f(n);    int m=sqrt(n)+1,ans=0;    for(int L=1;L<m;L++)    {        LL f1=n/L,f2=n/(L+1);        if(f1<MAXN)            f1=Sf[f1];        else            f1=_Sf[L];        if(f2<MAXN)            f2=Sf[f2];        else            f2=_Sf[L+1];        LL u=S1(L);        u=u*u%P;        ans+=u*(f1-f2+P)%P;        if(ans>=P)ans-=P;    }    for(int i=(int)(n/m);i;i--)    {        LL u=S1(n/i);        u=u*u%P;        ans+=u*f[i]%P;        if(ans>=P)ans-=P;    }    printf("%d\n",ans);    return 0;}