【LOJ6229】这是一道简单的数学题

【LOJ6229】这是一道简单的数学题

by AmanoKumiko

Description

求

\[F(n)=Σ_{i=1}^nΣ_{i=1}^i\frac{lcm(i,j)}{gcd(i,j)} \]

Input

输入一行，一个正整数 n

Output

输出F(n)，对\(10^9+7\)取模

Sample Input

5

Sample Output

84

Data Constraint

\(1≤n≤10^9\)

Solution

《简单》

\[F(n)=\frac{2F(n)}{2}=\frac{1}{2}(Σ_{i=1}^nΣ_{i=1}^n\frac{lcm(i,j)}{gcd(i,j)}+Σ_{i=1}^n\frac{lcm(i,i)}{gcd(i,i)}) \]

\[F(n)=\frac{1}{2}(Σ_{i=1}^nΣ_{i=1}^n\frac{lcm(i,j)}{gcd(i,j)}+n) \]

\[S(n)=Σ_{i=1}^nΣ_{i=1}^n\frac{lcm(i,j)}{gcd(i,j)} \]

\[S(n)=Σ_{i=1}^nΣ_{i=1}^n\frac{ij}{gcd(i,j)^2} \]

\[S(n)=Σ_{d=1}^nΣ_{i=1}^nΣ_{j=1}^n[gcd(i,j)=d]\frac{ij}{d^2} \]

\[S(n)=Σ_{d=1}^nΣ_{i=1}^{\lfloor{\frac{n}{d}}\rfloor}Σ_{j=1}^{\lfloor{\frac{n}{d}}\rfloor}[i⊥j]ij \]

\[S(n)=Σ_{d=1}^nΣ_{i=1}^{\lfloor{\frac{n}{d}}\rfloor}φ(i)i^2 \]

\[S(n)=Σ_{i=1}^nφ(i)i^2{\lfloor{\frac{n}{i}}\rfloor} \]

对\(S\)数论分块，故现在要求解的是\(Σφ(i)i^2\)的部分，杜教筛即可

Code

#include<bits/stdc++.h>
#include<tr1/unordered_map>
using namespace std;
#define F(i,a,b) for(int i=a;i<=b;i++)
#define Fd(i,a,b) for(int i=a;i>=b;i--)
#define LL long long
#define mo 1000000007
#define N 1000000

tr1::unordered_map<int,int>h;
int n,tot,prime[N+10],vis[N+10];
LL phi[N+10],s[N+10],i2,ans;

LL calc(LL x){
	LL v[3];
	v[0]=x;v[1]=x+1;v[2]=2*x+1;
	F(i,0,1)if(!(v[i]&1))v[i]/=2;
	F(i,0,2)if(!(v[i]%3))v[i]/=3;
	return v[0]*v[1]%mo*v[2]%mo%mo;
}

LL sum(LL x){
	if(x<N)return s[x];
	if(h[x])return h[x];
	LL res=(x*(x+1)/2%mo)*(x*(x+1)/2%mo)%mo;
	for(LL l=2,r;l<=x;l=r+1){
		r=x/(x/l);
		LL v1=calc(l-1),v2=calc(r);
		res=(mo+res-(v2-v1+mo)%mo*sum(x/l)%mo)%mo;
	}
	return (h[x]=res);
}

LL mi(LL x,LL y){
	if(y==1)return x;
	return y%2?x*mi(x*x%mo,y/2)%mo:mi(x*x%mo,y/2);
}

int main(){
	scanf("%d",&n);
	vis[1]=phi[1]=1;s[1]=1;
	F(i,2,N){
		if(!vis[i])prime[++tot]=i,phi[i]=i-1;
		for(int j=1;j<=tot&&i*prime[j]<=N;j++){
			vis[i*prime[j]]=1;
			if(!(i%prime[j])){phi[i*prime[j]]=phi[i]*prime[j];break;}
			else phi[i*prime[j]]=phi[i]*(prime[j]-1);
		}
		s[i]=(s[i-1]+phi[i]*i%mo*i%mo)%mo;
	}
	i2=mi(2,mo-2);
	for(int l=1,r;l<=n;l=r+1){
		r=n/(n/l);
		LL v1=sum(l-1),v2=sum(r);
		(ans+=(v2-v1+mo)%mo*(n/l)%mo)%=mo;
	}
	printf("%lld",(ans+n)%mo*i2%mo);
	return 0;
}