51nod 1237 最大公约数之和 V3

51nod 1237 最大公约数之和 V3

原题链接：
https://www.51nod.com/onlineJudge/submitDetail.html#!judgeId=353365

题面错误。原意是计算： $$G=\sum_{i=1}^n\sum_{j=1}^ngcd(i,j)$$

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

应用经典反演有： $$f(d)=\sum_{i=1}^n\sum_{j=1}^n[gcd(i,j)=d]\\F(d)=\sum_{d|a}f(a)=\sum_{d|i,i\leq n}\sum_{d|j,j\leq n}1=\big\lfloor\frac{n}{d}\big\rfloor^2\\f(d)=\sum_{d|a}\mu\big(\frac{a}{d}\big)F(a)$$

方法1（复杂度较高）

此时不在化简。之间对 $$\sum_{d=1}^ndf(d)$$ 计算。

记$f(d)=A(\lfloor\frac{n}{d}\rfloor)$

这是因为： $$f(d)=\sum_{d|a}\mu\big(\frac{a}{d}\big)F(a)\\=\sum_{a=1}^{\lfloor\frac{n}{d}\rfloor}\mu(a)\Big\lfloor\frac{\lfloor\frac{n}{d}\rfloor}{a}\Big\rfloor^2$$

记，$A$的定义为： $$A(n)=\sum_{a=1}^{n}\mu(a)\Big\lfloor\frac{n}{a}\Big\rfloor^2$$

对$A$分段计算即可。需要计算$M$(梅藤斯函数)

$$M(n)=\sum_{i=1}^n\mu(i)$$

因为：$\mu*I=e$

所以： $$M(n)=1-\sum_{i=2}^nM\big(\Big\lfloor\frac n i\Big\rfloor\big)$$

方法2

$$G=\sum_{d=1}^nd\sum_{d|a}\mu\big(\frac{a}{d}\big)F(a)\\=\sum_{d=1}^nd\sum_{d|a}f(d)=\sum_{d|a}\mu\big(\frac{a}{d}\big)\Big\lfloor\frac n a\Big\rfloor^2\\=\sum_{a=1}^n\Big\lfloor\frac n a\Big\rfloor^2\sum_{d|a}\mu\big(\frac{a}{d}\big)d\\=\sum_{a=1}^n\Big\lfloor\frac n a\Big\rfloor^2\varphi(a)$$

其中： $$\varphi*I=N$$

所以： $$S_\varphi(n)=\frac{n(n+1)}{2}-\sum_{i=2}^nS_\varphi\big(\Big\lfloor\frac n i\Big\rfloor\big)$$

给出第二种方法的代码：

#include <algorithm>
#include <string.h>
#include <stdio.h>
#include <cmath>
#define MAXN 2111111
#define SQR 111111

using namespace std;
typedef long long LL;
const int P= 1e9+7;
const int IV= P-P/2;

int phi[MAXN];
int Sp[MAXN];
int _Sp[SQR];


void init()
{
    for(int i=1,j;i<MAXN;Sp[i]=(Sp[i-1]+phi[i])%P,i++)
        for(j=i<<1,phi[i]=i-phi[i];j<MAXN;j+=i)phi[j]+=phi[i];
}

int S1(LL n)
{
    n%=P;
    return n*(n+1)%P*IV%P;
}

void clat_P(LL n,int d)
{
    int m=sqrt(n)+1,ans=0;
    for(int L=1;L<m;L++)
    {
        LL u=(n/L-n/(L+1))%P;
        ans+=u*Sp[L]%P;
        if(ans>=P)ans-=P;
    }
    for(int i=(int)(n/m);i>1;i--)
    {
        LL u=n/i;
        if(u<MAXN)
            u=Sp[u];
        else
            u=_Sp[i*d];
        ans+=u;
        if(ans>=P)ans-=P;
    }
    _Sp[d]=(S1(n)-ans+P)%P;
}

void slove_P(LL n)
{
    if(n<MAXN)return;
    for(int i=(int)(n/MAXN);i;i--)clat_P(n/i,i);
}
int slove(LL n)
{
    slove_P(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=Sp[f1];
        else
            f1=_Sp[L];
        if(f2<MAXN)
            f2=Sp[f2];
        else
            f2=_Sp[L+1];
        LL u=L;
        u=u*u%P;
        ans+=(f1-f2+P)*u%P;
        if(ans>=P)ans-=P;
    }
    for(int i=(int)(n/m);i;i--)
    {
        LL u=n/i;
        u%=P;
        u=u*u%P;
        ans+=u*phi[i]%P;
        if(ans>=P)ans-=P;
    }
    return ans;
}
int main ()
{
    init();
    LL n;
    scanf("%lld",&n);
    printf("%d\n",slove(n));
}