012. 高度可除的三角数(highly divisible triangular number)

三角数即由依次排列的自然数的和构成,所以第7个三角数是\(1 + 2 + 3 + 4 + 5 + 6 + 7 = 28\),前十个三角数是:\(1,3, 6, 10, 15, 21, 28, 36, 45, 55,\cdots\),让我们列出前七个三角数的因子:

$$ \begin{equation} \begin{aligned} 1&:1\\ 3&:1,3\\ 6&: 1,2,3,6\\ 10&: 1,2,5,10\\ 15&: 1,3,5,15\\ 21&: 1,3,7,21\\ 28&: 1,2,4,7,14,28 \end{aligned} \end{equation} $$
可以看出28是第一个因子超过5的三角数,求第一个因子超过五百万的三角数。

分析:首先根据初等数论的知识,我们知道一个数的因子个数等于其各质因数的幂分别加一并相乘,如\(20=2^2\times5^1\),则其因子数量为\((2+1)(1+1)=6\)。其次,我们知道三角数的通项是\(T_n=n(n+1)/2\),因为相邻的两个数\(n\)\(n+1\)必然是互质的,即两者没有共同的质因数,则有:

$$ n=p_1^{e_1}p_2^{e_2}\cdots p_i^{e_i},\ n+1=q_1^{f_1}q_2^{f_2}\cdots q_j^{f_j} $$

因此\(n\times(n+1)\)的质因数分解形式即为:

$$ n(n+1)=p_1^{e_1}p_2^{e_2}\cdots p_i^{e_i}\times q_1^{f_1}q_2^{f_2}\cdots q_j^{f_j} $$

因此其因子数为:

$$ (e_1+1)(e_2+1)\cdots(e_i+1)\times(f_1+1)(f_2+1)\cdots(f_j+1) $$

考虑到三角数还需要在\(n(n+1)\)的基础上再除以二,而\(n\)\(n+1\)中必然有一个为偶数,则二的幂应该减去一,则三角数的因子数应该为:

$$ e_1(e_2+1)\cdots(e_i+1)\times(f_1+1)(f_2+1)\cdots(f_j+1) $$

因此为了求三角数的因子数,我们只需要求\(n\)\(n+1\)的因子数再相乘即可,因为这两个数明显要小于它们的乘积形成的三角数,所以可以更快的进行质因数分解。我们从\(n=7\)开始迭代,每次加一,计算相应的三角数的因子数。因为这次迭代中的\(n+1\)会成为下次迭代中的\(n\),因此可以把中间值保存下来避免重复计算。当迭代计算出来的因子数量超过了五百时,迭代终止,此时返回相对应的三角数即为所求。

# time cost = 127 ms ± 1.09 ms

from sympy import factorint

def number_of_divisor(n):
    d = factorint(n)
    res = 1
    for i in d:
        res *= (d[i]+1) if i%2==1 else d[i]
    return res

def main(): 
    n = 7
    nd,nnd = number_of_divisor(n),number_of_divisor(n+1)
    while nd*nnd <= 500:
        n += 1
        nd,nnd = nnd,number_of_divisor(n+1)
    return n*(n+1)//2