Testing Fibonacci conjecture

Question

The sequence of Fibonacci numbers is defined as F1 = 1, Fn = Fn−2 + Fn−1. It has been conjectured that for any Fibonacci number F, F² + 41 is composite.

... [T]ask is to either prove the conjecture to be true or find a counter-example that demonstrates it is false.

(source)

Here is my working code for this problem - although I don't know the conjecture is either true or false in general.

I tested it for first 10⁷ natural numbers, also for first 1000 Fibonacci numbers - and it seems conjecture is true. but I takes a lot of time (80-90 seconds) to spit out a result.

As you see in last three functions, I have 3 method to test conjecture.

• First one (fibonacci_conjecture) generates a sequence of Fib²+41s and returns result of primality test.

• Second one is similar to first, but it first checks if numbers in a given range are Fibonacci numbers or not.

• Third one seems more optimal. As we know from definition, for a given Fibonacci number f, 5f²+4 or 5f²-4 is perfect square. If we add 41 to both of them we'll have 5f²+45 or 5f²+37. The first one is always divisible by 5, so there's no need to check these values.

How can this code be written in an optimal way, for faster results? Or, which techniques should I follow (either within the code or within the math)?

import math

def is_perfect_square(x):
    return int(math.sqrt(x)) == math.sqrt(x)

def is_fibo(x):
    c = 5*x*x
    return is_perfect_square(c+4) or is_perfect_square(c-4)

def is_prime(x):
    c = int(math.sqrt(x))+2
    if (x%2==0 and x!=2) or is_perfect_square(x): return False
    else:
        for i in range(2, c):
            if x%i==0: return False
    return True

def minus_four(fib):
    if is_fibo(fib) and fib>0: return is_perfect_square(5*fib*fib-4)

def fibo(x):
    if x<2: return 1
    else: return fibo(x-1)+fibo(x-2)

def fibonacci_conjecture(end):
    return (True not in map(is_prime, [fibo(i)*fibo(i)+41 for i in range(end+1)]))

def second_way(end):
    return (True not in map(is_prime, [i**2+41 for i in range(end+1) if is_fibo(i)]))

def third_way(end):
    return (True not in map(is_prime, [fibo(i)*fibo(i)+41 for i in range(end+1) if minus_four(i)]))

There is also a brute-force like implementation. that checks if 5x²+37 is prime, then checks if that number is a Fibonacci number.

from fibonacci_conjecture import * # functions above
import time
from sys import argv

script, start, end = argv

_list = []

start_time = time.time()
for i in range(int(start), int(end)):
    if is_prime(5*i*i+37) and is_fibo(5*i*i+37): _list.append(i)
print _list
print time.time() - start_time

Results calling the last one from Terminal:

$ python timings.py 0 1000
[]
0.0163550376892
$ python timings.py 1000 2000
[]
0.0392298698425
$ python timings.py 2000 12000
[]
1.66012907028
$ python timings.py 12000 42000
[]
35.4742548466
$ python timings.py 42000 43000
[]
2.01177597046
$ python timings.py 43000 48000
[]
11.1937119961

Answer 1

As always, the first option is to use PyPy:

>>> time python2 p.py 0 30000
[]
14.4747481346
python2 p.py 0 30000  14.50s user 0.01s system 100% cpu 14.500 total

>>> time pypy p.py 0 30000 
[]
0.990810871124
pypy p.py 0 30000  1.06s user 0.03s system 99% cpu 1.098 total

You shouldn't use star imports. Instead directly import the things you want:

from fibonacci_conjecture import is_fibo, is_prime

You should wrap the main body in a function called main.

The iteration would be better as a generator:

def get_falsifiers(start, end):
    for i in range(int(start), int(end)):
        if is_prime(5*i*i+37) and is_fibo(5*i*i+37):
            yield i

This would be run as:

_list = list(get_falsifiers(start, end))

Note that _list is a poor naming choice since it refers to the type and not the value. Something like falsifiers would be better.

These:

def fibonacci_conjecture(end):
    return (True not in map(is_prime, [fibo(i)*fibo(i)+41 for i in range(end+1)]))

def second_way(end):
    return (True not in map(is_prime, [i**2+41 for i in range(end+1) if is_fibo(i)]))

def third_way(end):
    return (True not in map(is_prime, [fibo(i)*fibo(i)+41 for i in range(end+1) if minus_four(i)]))

would be better as

def fibonacci_conjecture(end):
    return not any(map(is_prime, [fibo(i)*fibo(i)+41 for i in range(end+1)]))

def second_way(end):
    return not any(map(is_prime, [i**2+41 for i in range(end+1) if is_fibo(i)]))

def third_way(end):
    return not any(map(is_prime, [fibo(i)*fibo(i)+41 for i in range(end+1) if minus_four(i)]))

Further one should flatten the list comprehensions and map:

def fibonacci_conjecture(end):
    return not any([is_prime(fibo(i)*fibo(i)+41) for i in range(end+1)])

def second_way(end):
    return not any([is_prime(i**2+41) for i in range(end+1) if is_fibo(i)])

def third_way(end):
    return not any([is_prime(fibo(i)*fibo(i)+41) for i in range(end+1) if minus_four(i)])

You also then can just use a generator comprehension:

def fibonacci_conjecture(end):
    return not any(is_prime(fibo(i)*fibo(i)+41) for i in range(end+1))

def second_way(end):
    return not any(is_prime(i**2+41) for i in range(end+1) if is_fibo(i))

def third_way(end):
    return not any(is_prime(fibo(i)*fibo(i)+41) for i in range(end+1) if minus_four(i))

Your fibo(i)*fibo(i) requires two calls to fibo - this can be one with fibo(i) ** 2. Consider also using PEP 8 spacing. You can also remove the if in place and an and:

def fibonacci_conjecture(end):
    return not any(is_prime(fibo(i) ** 2 + 41) for i in range(end+1))

def second_way(end):
    return not any(is_fibo(i) and is_prime(i ** 2 + 41) for i in range(end+1))

def third_way(end):
    return not any(minus_four(i) and is_prime(fibo(i) ** 2 + 41) for i in range(end+1))

Some of these names (eg. minus_four) are particularly cryptic.

Your fibo actually takes \$\mathcal{O}(e^n)\$ time - what you really want as a basis is a generator to let you do

return all(not is_prime(f ** 2 + 41) for f in fibonacci_numbers(count))

This generator is typically written something like

def fibonacci_numbers(count):
    prev = 1
    curr = 0

    for _ in range(count):
        yield curr
        prev, curr = curr, prev + curr

Now each number takes \$\mathcal{O}(1)\$ (approximately), which is much faster.

Doing this actually quickly leads to an

OverflowError: long int too large to convert to float

since math.sqrt(x) needs to convert

400624556078135004801005131049868376277487458971454665625744044750241868509910118684223898375755770833486910716679881483356186192835924621611785526885928263384281394311273430963137957057246385921466246089044680129333479765771922731594997403327473526588890161488739900887481571279985972340702546216307221018066

to a float.

One other problem is also precision issues. For large numbers your is_perfect_square fails:

is_perfect_square(100000000**2 + 1)

since

float(100000000**2 + 1) == 1e+16

which is rounded. So the results are actually meaningless.