神犇和蒟蒻
根据
μ
(
n
)
\mu(n)
μ(n) 的定义，当
n
n
n 为大于
1
1
1的完全平方数时
μ
(
n
)
=
0
\mu(n)=0
μ(n)=0，当
n
=
1
n=1
n=1 时
μ
(
n
)
=
1
\mu(n)=1
μ(n)=1。因此
∑
i
=
1
n
μ
(
i
2
)
=
1
\sum_{i=1}^n \mu(i^2)=1
∑i=1nμ(i2)=1。
对于
φ
(
n
)
\varphi(n)
φ(n)，根据经验，转化为约数的
φ
\varphi
φ 之和。即：
n
=
∑
i
∣
n
φ
(
i
)
n=\sum_{i|n}\varphi(i)
n=∑i∣nφ(i)
因此有：
∑
i
=
1
n
∑
j
∣
i
φ
(
j
)
i
=
∑
i
=
1
n
i
2
=
1
6
n
(
n
+
1
)
(
2
n
+
1
)
\sum_{i=1}^n \sum_{j|i}\varphi(j)i=\sum_{i=1}^n i^2=\frac{1}{6}n(n+1)(2n+1)
∑i=1n∑j∣iφ(j)i=∑i=1ni2=61n(n+1)(2n+1)
换一种表达方法，先枚举约数，再枚举倍数：
∑
i
=
1
n
∑
i
∣
j
φ
(
i
)
j
=
∑
i
=
1
n
φ
(
i
)
i
∑
j
=
1
[
n
/
i
]
j
=
1
6
n
(
n
+
1
)
(
2
n
+
1
)
\sum_{i=1}^n \sum_{i|j}\varphi(i)j=\sum_{i=1}^n \varphi(i)i \sum_{j=1}^{[n/i]}j=\frac{1}{6}n(n+1)(2n+1)
∑i=1n∑i∣jφ(i)j=∑i=1nφ(i)i∑j=1[n/i]j=61n(n+1)(2n+1)
不要忘了求解的目标：
∑
i
=
1
n
φ
(
i
)
i
\sum_{i=1}^n \varphi(i)i
∑i=1nφ(i)i
将其从式子中提出来，移项：
∑
i
=
1
n
φ
(
i
)
i
=
1
6
n
(
n
+
1
)
(
2
n
+
1
)
+
∑
i
=
1
[
n
/
2
]
φ
(
i
)
i
(
1
−
∑
j
=
1
[
n
/
i
]
j
)
\sum_{i=1}^n \varphi(i)i=\frac{1}{6}n(n+1)(2n+1)+\sum_{i=1}^{[n/2]}\varphi(i)i(1-\sum_{j=1}^{[n/i]}j)
∑i=1nφ(i)i=61n(n+1)(2n+1)+∑i=1[n/2]φ(i)i(1−∑j=1[n/i]j)
因此最终的答案为：
1
6
n
(
n
+
1
)
(
2
n
+
1
)
+
∑
i
=
1
[
n
/
2
]
φ
(
i
)
i
(
1
−
1
2
[
n
/
i
]
(
[
n
/
i
]
+
1
)
\frac{1}{6}n(n+1)(2n+1)+\sum_{i=1}^{[n/2]}\varphi(i)i(1-\frac{1}{2}[n/i]([n/i]+1)
61n(n+1)(2n+1)+∑i=1[n/2]φ(i)i(1−21[n/i]([n/i]+1)
根据式子跑杜教筛即可。
时空复杂度
O
(
n
2
3
)
O(n^{\frac{2}{3}})
O(n32)
#include<iostream>
#include<map>
using namespace std;
#define R register int
#define L long long
#define N 10000001
#define P 1000000007
int prime[664579],pre[N];
map<int,int>Q;
inline int Solve(int n){
if(n<N){
return pre[n];
}
if(Q.count(n)!=0){
return Q[n];
}
int res=(1ll+n)*n%P*(n<<1|1)%P*166666668%P,r,lt=0,cur;
for(R i=1;i<=n>>1;i=r+1){
int tem=n/i;
r=n/tem;
cur=Solve(r);
res+=(1-((1ll+tem)*tem>>1)%P)*(cur-lt)%P;
if(res<0){
res+=P;
}else if(res>=P){
res-=P;
}
lt=cur;
}
Q[n]=res;
return res;
}
int main(){
pre[1]=1;
int n,ct=0;
for(R i=2;i!=N;i++){
if(pre[i]==0){
pre[i]=i-1;
prime[ct]=i;
ct++;
}
for(R j=0;j!=ct&&i*prime[j]<N;j++){
if(i%prime[j]==0){
pre[i*prime[j]]=pre[i]*prime[j];
break;
}
pre[i*prime[j]]=pre[i]*(prime[j]-1);
}
pre[i]=((L)pre[i]*i+pre[i-1])%P;
}
cin>>n;
cout<<1<<endl<<Solve(n);
return 0;
}