http://www.51nod.com/onlineJudge/questionCode.html#!problemId=1238
设(A(n)=sumlimits_{i=1}^nfrac{in}{(i,n)}),则(ans=sumlimits_{i=1}^nleft(2A(i)-i
ight))
[egin{aligned}
A(n)=&nsum_{d|n}sum_{i=1}^{frac nd}ileft[left(i,frac nd
ight)=1
ight]\
=&frac n2sum_{d|n}left(varphi(d)d+[d=1]
ight)\
=&frac n2sum_{d|n}varphi(d)d+frac n2\
ans=&sum_{i=1}^nleft(2A(i)-i
ight)\
=&sum_{i=1}^nisum_{d|i}varphi(d)d\
=&sum_{i=1}^nisum_{d=1}^{leftlfloorfrac ni
ight
floor}varphi(d)d^2\
end{aligned}
]
设(S(n)=sumlimits_{i=1}^nvarphi(i)i^2),求出(Oleft(sqrt n
ight))不同下取整取值的(S)就可以算出答案了,所以现在重点是杜教筛(S(n))。
先让(f(n)=varphi(n)n^2)卷上(f(n)=n^2)
[sum_{i=1}^{n}sum_{d|i}varphi(d)d^2left(frac id
ight)^2=sum_{i=1}^ni^2sum_{d|i}varphi(d)=sum_{i=1}^ni^3
]
[egin{aligned}
sum_{i=1}^ni^3=&sum_{i=1}^{n}sum_{d|i}varphi(d)d^2left(frac id
ight)^2\
=&sum_{i=1}^ni^2sum_{d=1}^{leftlfloorfrac ni
ight
floor}varphi(d)d^2\
=&sum_{i=1}^ni^2Sleft(leftlfloorfrac ni
ight
floor
ight)\
end{aligned}
]
如果要求(S(n)),(S(n)=sumlimits_{i=1}^ni^3-sumlimits_{i=2}^ni^2Sleft(leftlfloorfrac ni
ight
floor
ight)),时间复杂度(Oleft(n^{frac 23}
ight))。
PS:(sumlimits_{i=1}^ni^3=left(frac{n(n+1)}{2}
ight)^2)
#include<cmath>
#include<cstdio>
#include<cstring>
#include<algorithm>
using namespace std;
typedef long long ll;
const int N = 100003;
const int B = 4641589;
const int p = 1000000007;
const int ni2 = 500000004;
const int ni6 = 166666668;
ll n;
bool notp[B + 1];
int phi[B + 1], sum[B + 1], num = 0, prime[B + 1];
void Euler_shai() {
sum[1] = phi[1] = 1;
for (int i = 2; i <= B; ++i) {
if (!notp[i]) prime[++num] = i, phi[i] = i - 1;
for (int j = 1; j <= num && prime[j] * i <= B; ++j) {
notp[prime[j] * i] = true;
if (i % prime[j] == 0) {
phi[prime[j] * i] = prime[j] * phi[i];
break;
} else
phi[prime[j] * i] = (prime[j] - 1) * phi[i];
}
sum[i] = (1ll * phi[i] * i % p * i % p + sum[i - 1]) % p;
}
}
int ps[B];
ll ndx;
int S(ll x) {return (ndx = n / x) <= B ? sum[ndx] : ps[x];}
void DJ_shai() {
for (ll i = n, y; i >= 1; i = n / (y + 1)) {
y = n / i;
if (y <= B) continue;
int &s = ps[i];
s = y % p * (y % p) % p * ((y + 1) % p) % p * ((y + 1) % p) % p * ni2 % p * ni2 % p;
for (ll j = 2, pre = 1, spre = 1, sj; j <= y; spre = sj, pre = j, ++j) {
j = y / (y / j);
sj = j % p * ((j + 1) % p) % p * ((j * 2 + 1) % p) % p * ni6 % p;
((s -= (sj - spre + p) % p * S(i * j) % p) += p) %= p;
}
}
}
main() {
scanf("%lld", &n);
Euler_shai();
int ans = 0;
DJ_shai();
for (ll i = 1, pre = 0; i <= n; pre = i, ++i) {
i = n / (n / i);
(ans += (i + pre + 1) % p * ((i - pre) % p) % p * ni2 % p * S(i) % p) %= p;
}
printf("%d
", ans);
return 0;
}