Numpy 寻找质因数

质因数是可以整除某个整数的质数。通过寻找质因数的方式破解密码看似几乎不可能，但如果使用正确的算法——Fermat因式分解法和NumPy，这将变得很容易。基本思路是使用如下公式把整数N分解为c和d两个数：

递归地应用这个因式分解法，直至得到需要的质因数。

具体步骤

该算法需要我们为上述公式中的a尝试选择若干取值。

1. 创建尝试值数组。

创建一个NumPy数组，避免使用循环语句，这是一个合理的做法。但要注意，不能创建太大的、消耗内存太多的数组。在我用的系统中，使用包含一百万个元素的数组看似是没问题的。

a = numpy.ceil(numpy.sqrt(n))
lim = min(n, LIM)
a = numpy.arange(a, a + lim)
b2 = a ** 2 - n 

这里使用ceil函数对其输入参数的数组元素向上取整。

2. 得到数组b的小数部分。

b = numpy.sqrt(b2)

我们将检查数组b2中的元素是否是某个整数的平方。使用NumPy中的modf函数，可以获得数组b的小数部分。

fractions = numpy.modf(numpy.sqrt(b2))[0]

3. 查找小数部分为0的数组元素。

使用NumPy中的where函数，获得数组fractions中取值为0的元素的索引值。

indices = numpy.where(fractions == 0)

4. 找到第一个小数部分为0的数组元素。

实际上，我们只需要第一个小数部分为0的数组元素。首先调用NumPy的take函数，把上一步骤得到的数组indices作为参数，获得数组a中小数部分为0的数组元素。然后需要使用NumPy的ravel函数，把得到的结果展开为一维数组。

a = numpy.ravel(numpy.take(a, indices))[0]

如果需要解决的问题是查找数字600851475143的最大质因数，完整代码如下。

import numpy

#13195的质因数是5、7、13和29。

#数字600851475143的最大质因数是多少？

N = 600851475143
LIM = 10 ** 6

def factor(n):
    #1. 创建尝试值数组
    a = numpy.ceil(numpy.sqrt(n))
    lim = min(n, LIM)
    a = numpy.arange(a, a + lim)
    b2 = a ** 2 - n

    #2. 检查数组b2中的元素是否是某个整数的平方
    fractions = numpy.modf(numpy.sqrt(b2))[0]

    #3. 查找小数部分为0的数组元素
    indices = numpy.where(fractions == 0)

    #4. 找到第一个小数部分为0的数组元素
    a = numpy.ravel(numpy.take(a, indices))[0]
    a = int(a)
    b = numpy.sqrt(a ** 2 - n) 
    b = int(b)
    c = a + b
    d = a - b

    if c == 1 or d == 1:
        return

    print c, d
    factor(c)
    factor(d)

factor(N)

这段代码的执行结果如下。

1234169 486847
1471    839
6857    71

小结

我们采用了递归的方式应用Fermat因式分解法，用到了NumPy中的ceil、modf、where、ravel和take函数。这些函数的功能描述如下。