description
给出一个数N(<= 10^10),输出小于等于N的所有数,两两之间的最大公约数之和。
相当于计算这段程序(程序中的gcd(i,j)表示i与j的最大公约数):
由于结果很大,输出Mod 1000000007的结果。
G=0;
for(i=1;i< N;i++)
for(j=1;j<=N;j++)
{
G = (G + gcd(i,j)) % 1000000007;
}
Solution
这道题想出来后被人告知是杜教筛……
前面的都很顺利,非常快的推出式子。
ans=∑i=1n∑j=1ngcd(i,j)=∑d=1nd∑i=1n/d⌊nid⌋∗⌊nid⌋∗μ(i)=∑T=1n⌊nT⌋∗⌊nT⌋∗∑d|Td∗μ(T/d)
推到这里,我们发现一个很喜欢的式子 ∑d|Td∗μ(T/d) ,因为我们知道 ϕ(n)=∑d|nd∗μ(n/d) 。
ans=∑T=1n⌊nT⌋∗⌊nT⌋∗ϕ(T)。
然后考虑分块,所以我们要求一段数的欧拉函数,这和之前做的一道 51nod1239 欧拉函数之和简直一模一样,这里就不再推导了,直接上结论。
设f(x)=∑xi=1ϕ(i)
那么
f(x)=(1+x)∗x/2−∑T=2xf(⌊x/T⌋)。
所以ans中欧拉函数再分块中用
f(j)−f(i−1)
算一下即可。复杂度O(
n2/3
)。
Code
#include<iostream>
#include<cmath>
#include<cstring>
#include<cstdio>
#include<algorithm>
#define ll long long
using namespace std;
const ll maxn=1e6+5,mo=1e9+7,mo2=5e8+4;
ll d[maxn],bz[maxn],p[maxn],h[maxn*10],f[maxn*10];
ll n,i,t,j,k,l,x,y,z,ans;
int hash(ll x){
ll t=x%maxn;
while (h[t] && h[t]!=x) t=(t+1)%maxn;
return t;
}
ll dg(ll n){
if (n<=maxn) return p[n];
ll t=n%mo,k=t*(t+1)%mo*mo2%mo,i=2,x=hash(n),y;
if (h[x]) return f[x];
while (i<=n){
t=n/(n/i);
k-=dg(n/i)*(t-i+1)%mo;i=t+1;
}
k=(k%mo+mo)%mo;h[x]=n;f[x]=k;return k;
}
int main(){
//freopen("data.in","r",stdin);
scanf("%lld",&n);p[1]=1;
for (i=2;i<=maxn;i++){
if (!bz[i]) d[++d[0]]=i,p[i]=i-1;
for (j=1;j<=d[0];j++){
if (i*d[j]>maxn) break;
bz[i*d[j]]=1;
if (i%d[j]==0){
p[i*d[j]]=p[i]*d[j];break;
}else p[i*d[j]]=p[i]*p[d[j]];
}
}
for (i=1;i<=maxn;i++)
p[i]+=p[i-1];
i=1;
while (i<=n){
t=n/(n/i);
k=(dg(t)-dg(i-1)+mo)%mo;l=n/i%mo;
ans=(ans+l*l%mo*k%mo)%mo;
i=t+1;
}
printf("%lld\n",ans);
}