51nod 1238 最小公倍数之和 V3

51nod 1238 最小公倍数之和 V3

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

题面错误。。。

题目的实际意思是：

$$G=\sum_{i=1}^n\sum_{j=1}^n lcm(i,j)$$

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

按照一般套路：

$$G=\sum_{i=1}^n\sum_{j=1}^n lcm(i,j)\\=\sum_{d=1}^n\sum_{i=1}^n\sum_{j=1}^n\frac{ij}{d}[gcd(i,j)=d]$$

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

$$G=\sum_{d=1}^n\frac{1}{d}\sum_{i=1}^n\sum_{j=1}^nij[gcd(i,j)=d]$$

对于： $$\sum_{i=1}^n\sum_{j=1}^nij[gcd(i,j)=d]$$

记 $$S_1(n)=\sum_{i=1}^ni\\S_2(n)=\sum_{i=1}^ni^2$$

反演有： $$\sum_{i=1}^n\sum_{j=1}^nij[gcd(i,j)=d]\\=\sum_{d|a}\mu\big(\frac{a}{d}\big)a^2S_1\big(\big\lfloor\frac{n}{a}\big\rfloor\big)$$

所以： $$G=\sum_{d=1}^n\frac{1}{d}\sum_{d|a}\mu\big(\frac{a}{d}\big)a^2S_1\big(\big\lfloor\frac{n}{a}\big\rfloor\big)\\=\sum_{a=1}^nS_1\big(\big\lfloor\frac{n}{a}\big\rfloor\big)a^2\sum_{d|a}\mu\big(\frac{a}{d}\big)\frac{1}{d}\\=\sum_{a=1}^nS_1\big(\big\lfloor\frac{n}{a}\big\rfloor\big)a\sum_{d|a}\mu(d)d$$

令 $$h(n)=\mu(n)n$$

令 $$f(n)=\sum_{d|n}h(d)d$$

显然： $$h*I=f\\h*id_1=e$$

卷积运算有： $$f*id_1*\mu=h*I*id_1*\mu=e$$

$$\varphi*I=id_1$$

所以： $$\varphi=id_1*\mu$$

所以： $$f=h*I=\varphi^-$$

一个思路就是快速计算 $$af(a)$$

的前缀和函数。

令 $$F(n)=nf(n)$$

显然$F$逆函数$F^-$ $$F^-(n)=n\varphi(n)$$

对于： $$\sum_{i=1}^nF^-(i)\\=\sum_{i=1}^ni\varphi(i)\\=\sum_{i=1}^ni\sum_{d|i}\mu\big(\frac{i}{d}\big)d\\=\sum_{d=1}^nd\sum_{d|i}\mu\big(\frac{i}{d}\big)i\\\sum_{d=1}^nd^2\sum_{d|i}h\big(\frac{i}{d}\big)$$

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

则： $$\sum_{i=1}^nF^-(i)=\sum_{d=1}^nd^2S_h\big(\big\lfloor\frac{n}{d}\big\rfloor\big)$$

前面的分析有： $$h*id_1=e$$

从 h <script type="math/tex" id="MathJax-Element-436">h</script>开始,反向计算需要的函数前缀和即可。

#include <algorithm>
#include <string.h>
#include <cmath>
#include <stdio.h>
#define MAXN 3111111
using 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;
}