「Project Euler 530」GCD of Divisors

Problem

Description

定义 f(n)

f(n)=\sum_{d | n}\left(d,\frac{n}{d}\right)

定义 F(n)

F(n)=\sum_{i=1}^{n}f(i)

F\left(10\right)=32, F\left(1000\right)=12776
F\left(10^{15}\right)

Solution

Analysis

答案即为:

\sum_{i=1}^{n}\sum_{d | i}\left(d,\frac{i}{d}\right)

考虑枚举 \left(d,\frac{i}{d}\right)

\begin{aligned} &\sum_{d=1}^{n}d\sum_{i=1}^{n / d^{2}}\sum_{d' | i}\left[\left(d',\frac{i}{d'}\right)=1\right]\\ =&\sum_{d=1}^{n}d\sum_{i=1}^{n / d^{2}}\sum_{d' | i}\sum_{d''=1}^{n / d^{2}}\left[d'' | d'\right]\left[d'' |\frac{i}{d'}\right]\mu(d'')\\ =&\sum_{d=1}^{n}\sum_{d''=1}^{n}d\mu(d'')\sum_{i=1}^{n / d^{2}}\sum_{d' | i}\left[d'' | d'\right]\left[d'' |\frac{i}{d'}\right]\\ =&\sum_{d=1}^{n}\sum_{d''=1}^{n}d\mu(d'')\sum_{i=1}^{n / (dd'')^{2}}\sum_{d' | i}1\\ =&\sum_{d=1}^{n}\sum_{d''=1}^{n}d\mu(d'')\sum_{i=1}^{n / (dd'')^{2}}\tau(i) \end{aligned}

发现 dd'' 作为一个整体对后半部分式子产生影响,考虑枚举 dd''

\begin{aligned} &\sum_{d=1}^{n^{2}}\left(\sum_{d' | d}d'\mu\left(\frac{d}{d'}\right)\right)\sum_{i=1}^{n / d^{2}}\tau(i)\\ =&\sum_{d=1}^{n^{2}}\left(I\ast\mu\right)(d)\sum_{i=1}^{n / d^{2}}\tau(i)\\ =&\sum_{d=1}^{n^{2}}\varphi(d)\sum_{i=1}^{n / d^{2}}\tau(i) \end{aligned}

注意到后半部分式子只在 0<n / d^{2} 时有贡献,即有贡献的 d 最多为 \left\lfloor\sqrt{n}\right\rfloor

故答案为

\sum_{d=1}^{\left\lfloor\sqrt{n}\right\rfloor}\varphi(d)\sum_{i=1}^{n / d^{2}}\tau(i)

计算前 \sqrt{n} \varphi(n) 可以直接线性筛,时间复杂度 O\left(\sqrt{n}\right)
考虑计算 \sum_{i=1}^{n}\tau(i) ,枚举约数 d 可得:

\sum_{i=1}^{n}\tau(i)=\sum_{d=1}^{n}(n / d)

\sum_{i=1}^{n}\tau(i) 可以在 O\left(\sqrt{n}\right) 的时间复杂度内计算出来。
则暴力计算的复杂度为:

O\left(\sqrt{n}+\sum_{d=1}^{\left\lfloor\sqrt{n}\right\rfloor}\sqrt{\frac{n}{d^{2}}}\right)=O\left(\sqrt{n}+\sum_{d=1}^{\left\lfloor\sqrt{n}\right\rfloor}\frac{\sqrt{n}}{d}\right)=O\left(\sqrt{n}\ln{n}\right)

Code

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
#include<cstdio>

using i64=long long;
#define N 32000001

int phi[N];
bool notPrime[N];
int Prime[N>>3];
int cnt;

void Sieve(const int&Limit){
phi[1]=1;
for(int i=2;i<=Limit;++i){
if(!notPrime[i])
Prime[++cnt]=i,phi[i]=i-1;
for(int j=1,phii=phi[i];j<=cnt;++j){
int to=i*Prime[j];
if(Limit<to)
break;
notPrime[to]=true;
if(i%Prime[j])
phi[to]=phii*(Prime[j]-1);
else{
phi[to]=phii*Prime[j];
break;
}
}
}
}

inline i64 sum_sigma(const i64&n){
i64 Ans=0;
for(i64 d=1,nxd;d<=n;d=nxd+1){
nxd=n/(n/d);
Ans+=(n/d)*(nxd-d+1);
}return Ans;
}

int main(){
const i64 n=1000000000000000LL;
Sieve(N-1);
i64 Ans=0;
for(i64 d=1;d*d<=n;++d)
Ans+=(i64)phi[d]*sum_sigma(n/(d*d));
return(printf("%lld",Ans),0);
}