044. 五边形数(pentagon numbers)

五边形数可以由以下公式生成:\(P_n=n(3n-1)/2\),前十个五边形数如下:

$$ 1, 5, 12, 22, 35, 51, 70, 92, 117, 145, ... $$
可以看出\(P_4+P_7=22+70=92=P_8\)。然而,这两个数的差\(70-22=48\)并不是五边形数。假设有一对五边形数\(P_j\)\(P_k\),这两个数的和和差的绝对值都是五边形数,假设\(D=|P_k-P_j|\),求\(D\)的最小值。

分析:这道题的难度并不在于求出一个\(D\),而在于证明求出来的这个\(D\)是最小的,很多网络文章上的算法都忽视了这一点,默认算出来的第一个\(D\)就是最小的,而没有给出令人信服的证明。我只用了十几分钟就求出了答案,但是花了几天时间不断修改算法,才最终让我相信求出来的第一个\(D\)必然是的最小的\(D\)。看到这个题的第一个想法就是生成一个足够长的五边形数的列表,将列表中的数两两组合求其和与差,当和与差都为五边形数即满足了题目要求,差的绝对值即为所求。经过尝试,\(P_{2167}\)\(P_{1020}\)是满足条件的第一对五边形数,两者的差值就是题目的答案,但问题是如何证明这个差值就是最小的?

很多网络文章的思路是这样的,观察到相邻两个五边形数之间的差值在逐渐增大,事实上相邻两个五边形数的差值\(m=P_k-P_{k-1}=3k-2\)确实会随着\(k\)的增加而增加。因此,他们认为如果我们求出了一个\(D_1\),那么之后的\(D_2,D_3\cdots D_n\)都应该比\(D_1\)大,但是这种论证思路是不成立的。假设我们求出了一个\(D=P_k-P_j\)(\(k>j\)),虽然我们知道\(P_k\)会随着\(k\)的增大而增大,但是这并不能排除一对满足要求的值\(P_m\)\(P_n\),其中\(P_m>P_k\),但\((m-n)<(k-j)\),从而使得求出的\(D_{mn}<D_{kj}\)。为了完全排除这种可能性只有一种方法,那就是让\(k\)变得足够大,以使得相邻两个五边形数的差值都大于\(D\),如果经过验算表明在小于\(k\)的范围内\(D\)确实是最小值,而大于\(k\)的范围内不存在差值会小于\(D\),则\(D\)必然应是全局最小值。在已知答案的情况下,我们可以试着反推出这个足够大的\(k\)值,我们设相邻两个五边形数的差值大于题目的答案也就是最小的\(D\)值,即\(3k-2>5482660 \Rightarrow k>1827554\),这是一个非常大的值,验证在小于这个的数的范围内是满足要求的最小值要花很长的时间,至少在我的电脑上,我没法在可以接受的时间内验证这个想法。

因此,我们需要换一个思路,不从数对\(P_k\)\(P_j\)入手,而是从两者的和与差值入手。因为按照题目要求满足要求的数其和与差也是五边形数,所以我们只需要在这些五边形数中从小到大依次确定一个\(D\)值,再反推出相应的和,如果和与差都是五边形数,则是满足条件的数。因为我们是从小到大顺序来验证\(D\)值的,所以我们遇到的第一个满足要求的\(D\)值就必然是全局最小值,上面介绍的难题可以迎刃而解。具体算法是这样的,假设\(k>j\)\(P_k+P_j=S\),且有\(P_k-P_j=D\),则有\(S=D+P_j+P_j=D+2P_j\),这样当我们确定了\(D\)\(P_j\)后,只需要进一步确定\(D+P_j\)\(D+2P_j\)也是五边形数即找到了题目要求的数值。需要注意的是,我们并没有充分的理由相信\(P_j<D\),所以我们做验证的时候会先假设\(P_j<D\)算出\(D+2P_j\)进行验证,然后将两个数反转过来计算\(P_j+2D\)并验证,这样可以确保我们考察了所有可能的\((D,P_j)\)数对。

上面的问题弄清楚以后,剩下的只是一些平凡优化的技巧了。例如,为了确定一个数是否是五边形数,我们求出五边形数通项的反函数,即\(n=(1+\sqrt{24x+1})/6\),所以只有当且仅当我们得到\(1+\sqrt{24x+1}=0(mod\ 6)\)时,数\(x\)才是五边形数,使用这种方法,可以更快的判断一个数是否是五边形数。

from math import sqrt

def is_pantagon(x):
    if (sqrt(24*x+1)+1) % 6 == 0:
        return True
    return False

def main():
    pantagon = lambda n : n*(3*n-1)/2
    pt = [pantagon(n) for n in range(1,2000)]
    for i in range(1,2000):
        for j in range(0,i):
            s1,s2 = pt[i] + 2*pt[j], 2*pt[i] + pt[j]
            pk = pt[i] + pt[j]
            if is_pantagon(pk):
                if is_pantagon(s1):
                    return pt[i]
                if is_pantagon(s2):
                    return pt[j]