Prime number generation using Fibonacci possible?

Question

I'm generating prime numbers from Fibonacci as follows (using Python, with mpmath and sympy for arbitrary precision):

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)

It's a rather small algorithm but seemingly generates all primes (except for 5 somehow, but that's another question). I say seemingly because a very little percentage (0.5% under 1000 and 0.16% under 10K, getting less and less) isn't prime. For instance under 1000: 323, 377 and 442 are also generated. These numbers are not prime.

Is there something off in my script? I try to account for precision by relating the .dps setting to the number being calculated. Can it really be that Fibonacci and prime numbers are seemingly so related, but then when it's get detailed they aren't? :)

Answer 1

For this type of problem, you may want to look at the gmpy2 library. gmpy2 provides access to the GMP multiple-precision library which includes gcd() and fib() functions which calculate the greatest common divisor and the n-th fibonacci numbers quickly, and only using integer arithmetic.

Here is your program re-written to use 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)

You shouldn't be using any floating-point operations. You can calculate the GCD just using the builtin % (modulo) operator.

Update

As others have commented, you are checking for Fibonacci pseudoprimes. The actual test is slightly different than your code. Let's call the number being tested n. If n is divisible by 5, then the test passes if n evenly divides fib(n). If n divided by 5 leaves a remainder of either 1 or 4, then the test passes if n evenly divides fib(n-1). If n divided by 5 leaves a remainder of either 2 or 3, then the test passes if n evenly divides fib(n+1). Your code doesn't properly distinguish between the three cases.

If n evenly divides another number, say x, it leaves a remainder of 0. This is equivalent to x % n being 0. Calculating all the digits of the n-th Fibonacci number is not required. The test just cares about the remainder. Instead of calculating the Fibonacci number to full precision, you can calculate the remainder at each step. The following code calculates just the remainder of the Fibonacci numbers. It is based on the code given by @pts in 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)

It's written in pure Python and is very quick.

The sequence of numbers commonly known as the Fibonacci numbers is just a special case of a general Lucas sequence L(n) = p*L(n-1) - q*L(n-2). The usual Fibonacci numbers are generated by (p,q) = (1,-1). gmpy2.is_fibonacci_prp() accepts arbitrary values for p,q. gmpy2.is_fibonacci(1,-1,n) should match the results of the is_fib_pr(n) given above.

Disclaimer: I maintain gmpy2.

Answer 2

This isn't really a Python problem; it's a math/algorithm problem. You may want to ask it on the Math StackExchange instead.

Also, there is no need for any non-integer arithmetic whatsoever: you're computing floor(log10(x)) which can be done easily with purely integer math. Using arbitrary-precision math will greatly slow this algorithm down and may introduce some odd numerical errors too.

Here's a simple floor_log10(x) implementation:

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