[51Nod1237]最大公约数之和-V3

题意

  给定n$n$，求∑ni=1∑nj=1gcd(i,j)$\sum _{i=1}^n\sum _{j=1}^{n}gcd(i,j)$

解法

  直接化式子：

∑i=1n∑j=1ngcd(i,j)=∑d=1nd∑i=1n∑j=1n[gcd(i,j)=d]$$\sum _{i=1}^n\sum _{j=1}^{n}gcd(i,j)=\sum _{d=1}^{n}d\sum _{i=1}^n\sum _{j=1}^n[gcd(i,j)=d]$$

  下面有两种（或许不止两种）推导的方法，一种是莫比乌斯反演一下，这样做的复杂度会是O(n$n$)的，可以拿到一定的部分分，第二种方法就是对∑ni=1∑nj=1[gcd(i,j)=d]$\sum _{i=1}^n\sum _{j=1}^n[gcd(i,j)=d]$进行变化：

∑i=1n∑j=1n[gcd(i,j)=d]=∑i=1⌊nd⌋∑j=1⌊nd⌋[gcd(i,j)=1]=2∗∑i=1⌊nd⌋φ(i)−1=2∗ϕ(⌊nd⌋)−1$$\sum _{i=1}^n\sum _{j=1}^n[gcd(i,j)=d]=\sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\sum _{j=1}^{\lfloor \frac{n}{d}\rfloor }[gcd(i,j)=1]=2\ast \sum _{i=1}^{\lfloor \frac{n}{d}\rfloor }\phi (i)-1=2\ast \varphi (\lfloor \frac{n}{d}\rfloor )-1$$

  然后式子就变成了：∑nd=1d∗(2∗ϕ(⌊nd⌋)−1)$\sum _{d=1}^{n}d\ast (2\ast \varphi (\lfloor \frac{n}{d}\rfloor )-1)$，显然后面的部分有O(n‾√$\sqrt{n}$)种取值，那么对d$d$进行数论分块，然后求ϕ$\varphi$的过程就可以利用杜教筛了

复杂度

O(n23$n^{\frac{2}{3}}$)

代码

#include<iostream>
#include<cstdlib>
#include<cstdio>
#include<map>
#define Rint register int
#define Lint long long int
using namespace std;
const int m2=500000004;
const int N=20001000;
const int M=1e9+7;
bool vis[N];
int pri[N/10],cnt;
Lint phi[N];
map<Lint,int> f;
void Prepare()
{
    phi[1]=1;
    for(int i=2;i<N;i++)
    {
        if( !vis[i] )   pri[++cnt]=i,phi[i]=i-1;
        for(int j=1;j<=cnt;j++)
        {
            int x=pri[j]*i;
            if( x>=N )   break ;
            vis[x]=1;
            if( i%pri[j] )   phi[x]=phi[i]*phi[pri[j]];
            else
            {
                phi[x]=phi[i]*pri[j];
                break ;
            }
        }
    }
    for(int i=1;i<N;i++)   phi[i]=(phi[i]+phi[i-1])%M;
}
int cal(Lint n)
{
    if( n<N )   return phi[n];
    if( f.count(n) )   return f[n];
    Lint x=2,ret=n%M*(n%M+1)%M*m2%M;
    while( x<=n )
    {
        Lint y=n/(n/x);
        ret=(ret-1ll*(y-x+1)%M*cal(n/x)%M+M)%M,x=y+1;
    }
    return f[n]=ret;
}
int sum(Lint l,Lint r)
{
    l%=M,r%=M;
    return (l+r)%M*(r-l+1+M)%M*m2%M;
}
int main()
{
    Prepare();
    Lint n,ans=0,x=1;
    scanf("%lld",&n);
    while( x<=n )
    {
        Lint y=n/(n/x);
        ans=(ans+1ll*(2*cal(n/x)%M-1+M)%M*sum(x,y)%M)%M,x=y+1;
    }
    printf("%lld\n",ans);
    return 0;
}