Computing The First 10 Fibonacci Primes

Computational Number Theory Challenge

Fibonacci Numbers are 1, 1, 2, 3, 5, 8, 13, 21, 34, 55, 89,….

The next number in the Fibonacci sequence is obtained by adding together the two previous numbers. By convention the first two numbers in the sequence are each 1, then to get the third we do 1+1=2, to get the fourth we do 2+1=3, to get the fifth we do 3+2=5, and so on.

Prime Numbers are 2, 3, 5, 7, 11, 13, 17, …

This is any number greater than 1 which is divisible by ONLY itself and 1.

Fibonacci Primes are numbers that are BOTH Fibonacci and Prime. It turns out, such numbers are rare!

Indeed, there are only 10 Fibonacci Primes less than 1 billion.

This is quite interesting given there are 50,847,534 prime numbers less than 1 billion.

How could we compute the10 Fibonacci Primes less than 1 billion ourselves?

This is not computationally trivial since it might get resource-intensive to compute the prime and Fibonacci numbers up to 1 billion and then work out their intersection.

I will try this in Python.

Approach

My approach will be as follows:

1. Computing Fibonacci Numbers: Rather than try to store Fibonacci Numbers in a data structure, I’m going to use Python generators to allow us to iterate over these numbers up to 1 billion.
2. Computing Prime Numbers: I will construct a boolean array with 1 billion entries. The index represents the number and the value (True or False) represents whether the number is prime e.g. prime[0]=False, prime[1]=False, prime[2]=True, prime[3]=True, prime[4]=False, prime[5]=True and so on. The Sieve of Eratosthenes algorithm (discussed below) can be used to work out which numbers are prime.
3. Computing Fibonacci Primes: These can be obtained by iterating over the Fibonacci Numbers in 1) and then checking their index in 2) to see whether they are also prime.

Details and results are below.

Computing Fibonacci Numbers

Here is a script that allows us to get a generator object for Fibonacci Numbers. This should be self-explanatory:

def fibonacci(limit):
    a, b = 0, 1
    while a < limit:
        yield a
        a, b = b, a+b

fibonacci_generator = fibonacci(30)
for number in fibonacci_generator:
    print(number, end=' ')

This will print 1 1 2 3 5 8 13 21.

Computing Prime Numbers

The Sieve Of Eratosthenes is an ancient algorithm that can be used to work out the number of primes less than an integer N. It is referenced as far back as Ancient Greece to Eratosthenes of Cyrene.

It works like this: We take 2 as prime and cross off every multiple of 2 below N. We then check the next highest number after 2 which is not crossed off, this is 3, which is then prime. We then cross every multiple of 3. The next number not crossed off will be 5 which is then prime and so on…

Here is an animated image that illustrates this concept and works out the primes numbers less than 120:

Here is a script that implements this algorithm:

from math import isqrt
import numpy as np

upper_bound = 120
primes = np.full(upper_bound, True, dtype=bool)
primes[0] = False
primes[1] = False
for i in range(2, isqrt(upper_bound)+1):
    if primes[i]:
        for x in range(i*i, upper_bound, i):
            primes[x]=False

result = np.where(primes == True)
print(result)

This will print the prime numbers less than 120.

Here is a previous article that discusses details of this and how this is more efficient than more intuitive approaches. One clear drawback is that for very large N, large amounts of memory are required.

Fibonacci Primes < 1 Billion

Putting this together is straightforward and we get a script like this (Note: the timer function is making use of Python Decorators to calculate how long this takes to compute):

import time
from math import isqrt
import numpy

def timer(func):
    def wrapper(*args, **kwargs):
        before = time.time()
        result = func(*args, **kwargs)
        print('Function took:', time.time() - before, " seconds")
        return result
    return wrapper

def fibonacci(limit):
    a, b = 1, 1
    while a < limit:
        yield a
        a, b = b, a+b

@timer
def fibonacci_primes(upper_bound):

    primes = numpy.full(upper_bound, True, dtype=bool)
    primes[0] = False
    primes[1] = False
    for i in range(2, isqrt(upper_bound)+1):
        if primes[i]:
            for x in range(i*i, upper_bound, i):
                primes[x]=False

    fibonacci_generator = fibonacci(upper_bound)
    fibonacci_primes = []
    for i in fibonacci_generator:
        if primes[i]:
            fibonacci_primes.append(i)
    return fibonacci_primes

print(fibonacci_primes(1000000000))

This returns 2, 3, 5, 13, 89, 233, 1597, 28657, 514229 and 433494437 as the only Fibonacci Primes less than 1 billion!

On my machine, a Mac Mini (M1 2020 16GB Memory), this took 234 seconds or almost 4 minutes to compute.

Closing Remarks

It remains an open problem in Mathematics as to whether there are infinitely many Fibonacci Primes.

I’m interested in more efficient approaches to computing these numbers, please leave your comments below.

Revision: 5th November 2021

Here is a follow-up article which gets these numbers much faster and gets the time down to less than 0.01 seconds.

**Consider a Medium membership with my referral link here. You will be supporting my writing and there is no extra cost to you. Alternatively, you can buy me a coffee below. Many thanks!***