Description
求
∑i=1n∑j=1ngcd(i,j)
Solution
用反演,过程省略,
Ans=∑d=1nd∗∑i=1⌊nd⌋μ(i)⌊nid⌋2
设
T=id
Ans=∑T=1n⌊nT⌋2∑d|Tμ(d)∗Td
我们发现反演出了 φ ( 证明点这里)
Ans=∑T=1n⌊nT⌋2φ(T)
用一个分块加上杜教筛即可。
欧拉函数的杜教筛点这里
复杂度: O(n23)
Code
#include <cstdio>
#include <cstdlib>
#define fo(i,a,b) for(int i=a;i<=b;i++)
#define min(q,w) ((q)<(w)?(q):(w))
using namespace std;
typedef long long LL;
const int N=5000500,M=2544657,mo=1e9+7;
const LL eni=500000004;
LL n,ans;
bool prz[N+10];
int pr[N+10];
int phi[N+10];
int Hx[M+1][2];
int HX(LL q)
{
int i=q%M;
while(Hx[i][0]&&Hx[i][0]!=q)i=(i+1)%M;
return i;
}
LL Gphi(LL q)
{
if(q<=N)return phi[q];
int t=HX(q);
if(Hx[t][0])return Hx[t][1];
Hx[t][0]=q;
LL ans=0;
for(LL i=2,nx;i<=q;i=nx+1)
{
nx=q/(q/i);
ans=(ans+(nx-i+1)%mo*Gphi(q/i))%mo;
}
q%=mo;
return Hx[t][1]=(q*(q+1)%mo*eni%mo-ans+mo)%mo;
}
int main()
{
phi[1]=1;
fo(i,2,N)
{
if(!prz[i])pr[++pr[0]]=i,phi[i]=i-1;
fo(j,1,pr[0])
{
int t=pr[j]*i;
if(t>N)break;
prz[t]=1;
phi[t]=phi[i]*pr[j];
if(i%pr[j]==0)break;
phi[t]=phi[i]*(pr[j]-1);
}
}
fo(i,2,N)phi[i]=(phi[i]+phi[i-1])%mo;
scanf("%lld",&n);
for(LL i=1,nx;i<=n;i=nx+1)
{
nx=n/(n/i);
ans=(ans+(n/i)%mo*((n/i)%mo)%mo*(Gphi(nx)-Gphi(i-1))%mo)%mo;
}
printf("%lld\n",(ans+mo)%mo);
return 0;
}