021. 亲和数(Amicable numbers)

\(d(n)\)表示自然数\(n\)所有真因子(除开数\(n\)身的所有因子)的和,如果\(d(a)=b\)\(d(b)=a\),其中\(a\neq b\),那么\(a\)\(b\)便为亲和数对,其中的每个数称为亲和数。例如,220的真因子为1,2, 4, 5, 10, 11, 20, 22, 44, 55 和110,所以\(d(220)=284\),284的真因子为1,2, 4, 71与142,所以\(d(284)=220\),求10000以下所有亲和数的和。

分析:这道题的整体思路比接直接,编写一个求特定自然数所有真因子和的函数,然后在一至一万内遍历寻找满足题目要求的亲和数对,然后对这些亲和数对求和即可。这里的核心在在于如何更快的求出特定自然数N的真因子之和。一个比较直接的思路是在\([1,\sqrt{N}]\)遍历,当找到一个数\(d\)整除\(N\),则\(d\)\(N\)的因子,同时\(N/d\)是另外一个因子,将这些因子相加并减去数N本身即为数\(N\)的所有真因子之和。这种算法在小规模数据上表现良好,但仍然存在一个更好的算法。

\(\sigma(n)\)表示数\(n\)所有因子之和,即\(\sigma(n)=\sum_{d|n}d\),其中\(d|n\)表示\(d\)\(n\)的因子。则已知对于一个质数\(p\),其因子只有一和自身,则\(\sigma(p)=p+1\),而对于质数\(p\)\(a\)次方,其因子为:

$$ 1\ ,\ p\ ,\ p^2\ ,\ p^3\cdots p^a $$

\(\sigma(p^a)=1+p+p^2+p^3+\cdots p^a\),显然为一个等比数列,则根据等比数列公式:

$$ \sigma(p^a)=\frac{p^{a+1}-1}{p-1} $$

现在对于任何质数,我们可以求其真因子之和。对于合数可对其进行质因数分解:

$$ c=p_1^{e_1}p_2^{e_2}\cdots p_a^{e_a} $$

从上式易知,任意合数\(c\)可将其分解为两个互素的数\(a\)\(b\)的乘积,即\(c=ab,(a,b)=1\),易证:

$$ \sigma(c)=\sigma(ab)=\sigma(a)\cdot\sigma(b) $$

\(a\)\(b\)进行质因数分解,则有:

$$ \sigma \left( c \right) =\sigma \left( p_{1}^{e_1} \right) \sigma \left( p_{2}^{e_2} \right) \cdots \sigma \left( p_{a}^{e_a} \right) =\prod_{i=1}^a{\frac{p_i^{e_i+1}-1}{p_i-1}} $$

在上式的计算结果中减去数\(c\)即为其真因子之和。据此我们可以编写一个计算特定自然数所有真因子和的函数,然后在\([1,10^4]\)范围内寻找亲和数对然后加总即可。第二种思路的代码实现如下:

# time cost = 206 ms ± 2.44 ms

from sympy import factorint

def sum_of_divisors(n):
    divs = factorint(n)
    res = 1
    for p in divs:
        res *= (p**(divs[p]+1)-1)//(p-1)
    return res-n

def main(n=10000):
    arr = []
    for i in range(2,n+1):
        a = sum_of_divisors(i)
        if a != i and sum_of_divisors(a) == i:
            arr.append(i)
    return sum(arr)