连分数分解法寻找整数的因子(Python)

　　首先，先讲个小故事。
　　1903年10月，科尔（Cole）在纽约参加美国数学（AMS）的一个会议，议程里有他的一篇题目平淡的论文：“关于一个大数的分解”。大会主席请他宣读论文时，一向寡言的科尔走上黑板，一言不发，拿起粉笔就开始算2的67次方。然后，他小心地减去1.接着，还是不说一个字，他移到黑板的空白处，老老实实地计算

$$193707721\times761838257287$$ 两个计算结果一致……美国数学会的听众们破天荒地第一次也是唯一一次为他们听到的报告热烈鼓掌欢呼。科尔还是一言不发，回到自己的位置。也没人向他提任何问题。
　　几年过后，1911年，贝尔（E.T.Bell）问科尔用了多长时间来分解M67（Mersenne素数）。科尔回答：“三年里的星期天。”

---

　　在初等数论中，我们可以利用连分数分解法来寻找大整数的因子，具体的算法可以参考《 初等数论及其应用 美第五版》。以下我们给出Python代码：

# -*- coding: utf-8 -*-
from math import sqrt,floor,gcd
#判断是否为平方数
def is_square(n):
    t = int(sqrt(n))
    return t*t == n
def main():
    n = 2**79-1 #要分解的数
    print('n=%s'%('2**67-1'))
    seq_P, seq_Q=0, 1
    a0 = floor(sqrt(n))
    seq_a= [a0]
    seq_p = [0,a0]
    i = 1
    while True:
        #P_{k},Q_{k}序列
        seq_P = seq_a[0]*seq_Q-seq_P
        seq_Q = divmod(n-seq_P**2,seq_Q)[0]
        t = (seq_P+sqrt(n))/seq_Q
        seq_a[0] = floor(t)
        #p_{k}序列
        if i == 1:
            seq_p.append(a0*seq_a[0]+1)
        else:
            seq_p[0] = seq_p[1]
            seq_p[1] = seq_p[2]
            seq_p[2] = seq_a[0]*seq_p[1]+seq_p[0]
        #分解因式
        if i%2 == 0 and is_square(seq_Q):
            s = int(sqrt(seq_Q))
            factor = [gcd(seq_p[1]-s,n),gcd(seq_p[1]+s,n)]
            if factor[0] != 1 and factor[1] != 1:
                print('第%d项：%d,%d'%(i,factor[0],factor[1]))
                print('程序运行完毕！')
                break
        if i%1000 == 0:
            print('已运行到第%d项'%i)
        i += 1

main()

运行结果如下：

因此，M67有因子193707721和761838257287，因此M67不是梅森素数。

---

　　按照这个思路，我们再给出其他几个非梅森素数的 $M_{k}$，当然读者也可以自己测试。注意：连分数分解法不一定能找到整数的因子。
　　
　　
　　
　　
注意：本人现已开通两个微信公众号： 用Python做数学（微信号为：python_math）以及轻松学会Python爬虫（微信号为：easy_web_scrape）， 欢迎大家关注哦~~