基准时间限制:8 秒 空间限制:262144 KB 分值: 640 难度:8级算法题
出一个数N,输出小于等于N的所有数,两两之间的最小公倍数之和。
相当于计算这段程序(程序中的lcm(i,j)表示i与j的最小公倍数):
由于结果很大,输出Mod 1000000007的结果。
G=0;
for(i=1;i<=N;i++)
for(j=1;j<=N;j++)
{
G = (G + lcm(i,j)) % 1000000007;
}
Input
输入一个数N。(2 <= N <= 10^10)
Output
输出G Mod 1000000007的结果。
Input示例
4
Output示例
72
之前推过一个差不多的,但是那个的数量级是10的5次的,而且不是正方形,是给出的n和m,牵扯上了莫比乌斯函数的,但是这个就不能这样搞了,反正先推一推吧。
ans=∑i=1n∑j=1nijgcd(i,j)=∑d=1n∑i=1n∑j=1n[gcd(i,j)==d]ijd=∑d=1nd∑i=1⌊nd⌋∑j=1⌊nd⌋[gcd(i,j)==1]ij
a
n
s
=
∑
i
=
1
n
∑
j
=
1
n
i
j
g
c
d
(
i
,
j
)
=
∑
d
=
1
n
∑
i
=
1
n
∑
j
=
1
n
[
g
c
d
(
i
,
j
)
==
d
]
i
j
d
=
∑
d
=
1
n
d
∑
i
=
1
⌊
n
d
⌋
∑
j
=
1
⌊
n
d
⌋
[
g
c
d
(
i
,
j
)
==
1
]
i
j
设:
F(n)=∑i=1n∑j=1n[gcd(i,j)==1]ij
F
(
n
)
=
∑
i
=
1
n
∑
j
=
1
n
[
g
c
d
(
i
,
j
)
==
1
]
i
j
观察发现若i==j,且gcd为1,只有i和j均为1时,所以我们不妨考虑这个式子,设:
F′(x)=∑i=1n∑j=1i[gcd(i,j)==1]ij
F
′
(
x
)
=
∑
i
=
1
n
∑
j
=
1
i
[
g
c
d
(
i
,
j
)
==
1
]
i
j
易知 F(x)=2F′(x)−1 F ( x ) = 2 F ′ ( x ) − 1 ,所以我们先来考虑 F′(x) F ′ ( x )
那很显然啊…
F′(x)=12+∑i=1ni2ϕ(i)2
F
′
(
x
)
=
1
2
+
∑
i
=
1
n
i
2
ϕ
(
i
)
2
所以:
F(x)=∑i=1ni2ϕ(i)
F
(
x
)
=
∑
i
=
1
n
i
2
ϕ
(
i
)
.
感觉这个式子好简洁,那么就是求:
ans=∑i=1niF(⌊ni⌋)
a
n
s
=
∑
i
=
1
n
i
F
(
⌊
n
i
⌋
)
似乎没什么实质化简,用杜教搞一搞啦。
想到 ∑d|nϕ(d)=n ∑ d | n ϕ ( d ) = n ,看看能不能利用这个性质了。
设: f(x)=x2ϕ(x),g(x)=x2 f ( x ) = x 2 ϕ ( x ) , g ( x ) = x 2 ,我们拿这两个函数做狄利克雷卷积,有:
f∗g(n)=∑d|nd2ϕ(d)n2d2=n3=∑d|ng(d)f(nd)
f
∗
g
(
n
)
=
∑
d
|
n
d
2
ϕ
(
d
)
n
2
d
2
=
n
3
=
∑
d
|
n
g
(
d
)
f
(
n
d
)
对卷积和做一下前缀和:
∑i=1ni3=∑i=1n∑j|ig(j)f(ij)=∑j=1ng(j)∑k=1⌊nj⌋f(k)=∑j=1ng(j)F(⌊nj⌋)=n2(n+1)24
∑
i
=
1
n
i
3
=
∑
i
=
1
n
∑
j
|
i
g
(
j
)
f
(
i
j
)
=
∑
j
=
1
n
g
(
j
)
∑
k
=
1
⌊
n
j
⌋
f
(
k
)
=
∑
j
=
1
n
g
(
j
)
F
(
⌊
n
j
⌋
)
=
n
2
(
n
+
1
)
2
4
所以:
g(1)F(n)=n2(n+1)24−∑j=2ng(j)F(⌊nj⌋)
g
(
1
)
F
(
n
)
=
n
2
(
n
+
1
)
2
4
−
∑
j
=
2
n
g
(
j
)
F
(
⌊
n
j
⌋
)
推导部分就到这了,现在就是递归分块求了。
代码:
#include<iostream>
#include<cstdio>
#include<algorithm>
#include<cstring>
#include<map>
#include<cmath>
#define maxx 5000000
#define mod 1000000007
#define ll long long
using namespace std;
bool isP[maxx];
int prime[400000];
int cnt;
int phi[maxx];
ll F[maxx];
ll inv2=500000004,inv4=250000002,inv6=166666668;
ll get(ll x)
{
x%=mod;
return (x+1)*x%mod*(2*x+1)%mod*inv6%mod;
}
void init()
{
phi[1]=1;
for(int i=2;i<maxx;i++)
{
if(!isP[i]){prime[cnt++]=i;phi[i]=i-1;}
for(int j=0;j<cnt&&(ll)i*prime[j]<maxx;j++)
{
isP[i*prime[j]]=true;
if(i%prime[j])
phi[i*prime[j]]=phi[i]*(prime[j]-1);
else
{
phi[i*prime[j]]=phi[i]*prime[j];
break;
}
}
}
for(ll i=1;i<maxx;i++)
F[i]=i*i%mod*phi[i]%mod;
for(int i=1;i<maxx;i++)
F[i]=(F[i]+F[i-1])%mod;
}
map<ll,ll>M;
ll work(ll x)
{
if(x<maxx)return F[x];
if(M[x])return M[x];
ll t=x%mod;
ll ans=t*t%mod*(t+1)%mod*(t+1)%mod*inv4%mod;
for(ll i=2,last;i<=x;i=last+1)
{
last=x/(x/i);
ans=(ans-(get(last)-get(i-1))*work(x/i)%mod)%mod;
}
if(ans<0)ans=(ans+mod)%mod;
M[x]=ans;
return ans;
}
ll _get(ll a,ll b)
{
return (b-a+1)%mod*((a+b)%mod)%mod*inv2%mod;
}
int main()
{
init();
//cout<<cnt<<endl;
ll n;
cin>>n;
ll ans=0;
for(ll i=1,last;i<=n;i=last+1)
{
last=n/(n/i);
ans=(ans+_get(i,last)*work(n/i)%mod)%mod;
}
cout<<ans<<endl;
return 0;
}