使用斐波那契生成素数可能吗？

Question

我正在用以下方法从斐波那契数列中生成质数（使用 Python，并借助 mpmath 和 sympy 来处理任意精度）：

from mpmath import *

def GCD(a,b):
    while a:
        a, b = fmod(b, a), a
    return b

def generate(x):
    mp.dps = round(x, int(log10(x))*-1)
    if x == GCD(x, fibonacci(x-1)):
        return True
    if x == GCD(x, fibonacci(x+1)):
        return True

    return False

for x in range(1000, 2000)
    if generate(x)
        print(x)

这个算法相对简单，但似乎能生成所有的质数（除了 5，这又是另一个问题）。我说“似乎”是因为在生成的数中，有一小部分（在1000以下是0.5%，在10000以下是0.16%，而且这个比例越来越小）并不是质数。例如，在1000以下，323、377和442这些数也是被生成的，但它们并不是质数。

我的代码有什么问题吗？我试着通过将 .dps 设置与正在计算的数字关联起来，以确保精度。难道斐波那契数和质数之间真的有这么大的关系，但当细看时却又不是这样吗？ :)

Answer 1

这其实不是一个Python的问题，而是一个数学或算法的问题。你可能想在数学相关的论坛上问这个问题。

另外，计算这个问题时完全不需要用到非整数的运算：你要计算的是floor(log10(x))，这可以通过简单的整数运算来完成。使用任意精度的数学运算会大大减慢这个算法的速度，还可能引入一些奇怪的数字错误。

下面是一个简单的floor_log10(x)实现：

from __future__ import division # if using Python 2.x

def floor_log10(x):
  res = 0
  if x < 1:
    raise ValueError
  while x >= 1:
    x //= 10
    res += 1
  return res

Answer 2

对于这种问题，你可以看看 gmpy2 这个库。gmpy2 是一个可以使用GMP多精度库的工具，它里面有一些函数，比如gcd()和fib()，可以快速计算最大公约数和第n个斐波那契数，而且只用整数运算。

下面是用 gmpy2 重写的你的程序。

import gmpy2

def generate(x):
    if x == gmpy2.gcd(x, gmpy2.fib(x-1)):
        return True
    if x == gmpy2.gcd(x, gmpy2.fib(x+1)):
        return True
    return False

for x in range(7, 2000):
    if generate(x):
        print(x)

你不应该使用任何浮点运算。计算最大公约数可以直接用内置的 %（取余）运算符。

更新

正如其他人所说的，你在检查斐波那契伪素数。实际的测试和你的代码稍有不同。我们把要测试的数字叫做 n。如果 n 能被5整除，那么测试通过的条件是 n 能整除 fib(n)。如果 n 除以5的余数是1或4，那么测试通过的条件是 n 能整除 fib(n-1)。如果 n 除以5的余数是2或3，那么测试通过的条件是 n 能整除 fib(n+1)。你的代码没有正确区分这三种情况。

如果 n 能整除另一个数字，比如 x，那么余数就是0。这就等于说 x % n 是0。计算第 n 个斐波那契数的所有数字并不是必须的。测试只关心余数。你可以在每一步计算余数，而不是计算完整的斐波那契数。下面的代码只计算斐波那契数的余数。这个代码是基于@pts在 Python mpmath not arbitrary precision? 中提供的代码。

def gcd(a,b):
    while b:
        a, b = b, a % b
    return a

def fib_mod(n, m):
    if n < 0:
        raise ValueError

    def fib_rec(n):
        if n == 0:
            return 0, 1
        else:
            a, b = fib_rec(n >> 1)
            c = a * ((b << 1) - a)
            d = b * b + a * a
            if n & 1:
                return d % m, (c + d) % m
            else:
                return c % m, d % m

    return fib_rec(n)[0]

def is_fib_prp(n):
    if n % 5 == 0:
        return not fib_mod(n, n)
    elif n % 5 == 1 or n % 5 == 4:
        return not fib_mod(n-1, n)
    else:
        return not fib_mod(n+1, n)

这个代码是用纯Python写的，非常快。

大家熟知的斐波那契数列其实是一个更一般的卢卡斯数列 L(n) = p*L(n-1) - q*L(n-2) 的特例。通常的斐波那契数是通过 (p,q) = (1,-1) 生成的。gmpy2.is_fibonacci_prp() 可以接受任意的p和q值。gmpy2.is_fibonacci(1,-1,n) 应该和上面提到的 is_fib_pr(n) 的结果一致。

免责声明：我维护gmpy2。