感性理解后面部分是关于 \(i\) 的积性函数（设为 \(f(i)\)），然后考虑质数幂处的值。
\[ \begin{aligned} f(p^c) &= \sum\limits_{j=1}^{p^c} \sum\limits_{k=1}^{p^c} {\rm lcm}(\gcd(p^c,j),\gcd(p^c,k)) \\ &= \sum\limits_{i=0}^c p^i \sum\limits_{j=1}^{p^c} \sum\limits_{k=1}^{p^c} [{\rm lcm}(\gcd(p^c,j),\gcd(p^c,k)) = p^i] \\ &= \sum\limits_{i=0}^c p^i \sum\limits_{j=1}^{p^c} \sum\limits_{k=1}^{p^c} ([{\rm lcm}(\gcd(p^c,j),\gcd(p^c,k)) \mid p^i] - [{\rm lcm}(\gcd(p^c,j),\gcd(p^c,k)) \mid p^{i-1}]) \\ &= \sum\limits_{i=0}^c p^i \left[\left(\sum\limits_{d=0}^i \varphi(p^{c - d})\right)^2 - \left(\sum\limits_{d=0}^{i-1} \varphi(p^{c-d})\right)^2\right] \\ &= \sum\limits_{i=0}^c p^i \left[(p^c-[i < c]p^{c-i-1})^2 - (p^c-p^{c-i})^2\right] \\ &= (2c+1)(p^{2c}-p^{2c-1})+p^{c-1} \end{aligned} \]
然后 Min_25 筛即可。
代码： 1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
using namespace std;
const int MX = 1e5;
long long n;
int lim;
int vis[MX + 5],cnt,prime[MX + 5];
unsigned Gprime[MX + 5][2];
int tot,le[MX + 5],ge[MX + 5];
long long lis[2 * MX + 5];
unsigned G[2 * MX + 5][3],Fprime[2 * MX + 5];
inline int &id(long long x)
{
return x <= lim ? le[x] : ge[n / x];
}
unsigned F(int k,long long n)
{
if(n < prime[k] || n <= 1)
return 0;
unsigned ret = Fprime[id(n)] - (3 * Gprime[k - 1][1] - 3 * Gprime[k - 1][0] + (k - 1));
for(register int i = k;i <= cnt && (long long)prime[i] * prime[i] <= n;++i)
{
long long pw0 = 1,pw = prime[i],pw2 = (long long)prime[i] * prime[i];
for(register int c = 1;pw2 <= n;pw0 = pw,pw = pw2,pw2 *= prime[i],++c)
ret += ((2 * c + 1) * (pw * pw - pw * pw0) + pw0) * F(i + 1,n / pw) + ((2 * c + 3) * (pw2 * pw2 - pw2 * pw) + pw);
}
return ret;
}
int main()
{
freopen("sum.in","r",stdin),freopen("sum.out","w",stdout);
scanf("%lld",&n),lim = sqrt(n);
for(register int i = 2;i <= lim;++i)
{
if(!vis[i])
prime[++cnt] = i,Gprime[cnt][0] = Gprime[cnt - 1][0] + i,Gprime[cnt][1] = Gprime[cnt - 1][1] + i * i;
for(register int j = 1;j <= cnt && i * prime[j] <= lim;++j)
{
vis[i * prime[j]] = 1;
if(!(i % prime[j]))
break;
}
}
for(register long long l = 1,r;l <= n;l = r + 1)
{
r = n / (n / l);
lis[id(n / l) = ++tot] = n / l;
G[tot][0] = n / l - 1;
if(n / l & 1)
G[tot][1] = (n / l + 1) / 2 * (n / l) - 1;
else
G[tot][1] = (n / l) / 2 * (n / l + 1) - 1;
long long tmp[3];
tmp[0] = n / l,tmp[1] = n / l + 1,tmp[2] = 2 * (n / l) + 1;
for(register int i = 0;i < 3;++i)
if(!(tmp[i] % 2))
{
tmp[i] /= 2;
break;
}
for(register int i = 0;i < 3;++i)
if(!(tmp[i] % 3))
{
tmp[i] /= 3;
break;
}
G[tot][2] = tmp[0] * tmp[1] * tmp[2] - 1;
}
for(register int k = 1;k <= cnt;++k)
{
int p = prime[k];
long long s = (long long)p * p;
for(register int i = 1;lis[i] >= s;++i)
G[i][0] -= G[id(lis[i] / p)][0] - (k - 1),
G[i][1] -= (long long)p * (G[id(lis[i] / p)][1] - Gprime[k - 1][0]),
G[i][2] -= s * (G[id(lis[i] / p)][2] - Gprime[k - 1][1]);
}
for(register int i = 1;i <= tot;++i)
Fprime[i] = 3 * G[i][2] - 3 * G[i][1] + G[i][0];
printf("%u\n",F(1,n) + 1U);
}