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\)次方，其因子为：

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

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

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

将\(a\)与\(b\)进行质因数分解，则有：

在上式的计算结果中减去数\(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)