Fibonacci Sequence: Finding Composites Problem

Question

I am looking to find a pair of numbers with a GCD (Greatest Common Denominator) of 1, that the first N terms of the sequence X0, X1, ... XN are all composite.

For my code, for some reason, it gets stuck when i == 15, j == 878, and k == 78. It gets stuck when running is_prime() on the two last items in the list.

import math


def is_prime(num):
    if num < 2:
        return False
    for x in range(2, math.floor(math.sqrt(num)) + 1):
        if num % x == 0:
            return False
    return True


# create list containing a range of composite numbers
numbers = []
for i in range(4, 200):
    if not is_prime(i):
        numbers.append(i)

for i in numbers:
    found = False
    for j in numbers:
        if math.gcd(i, j) == 1:
            # print(i, "-", j, end=" | ")
            fibonacci = [i, j]
            contains_prime = False
            for k in range(2, 500):
                if is_prime(fibonacci[-2] + fibonacci[-1]):
                    contains_prime = True
                    break
                else:
                    fibonacci = [fibonacci[-1], fibonacci[-2] + fibonacci[-1]]
            if not contains_prime:
                print(i, j)
                found = True
        if found:
            break
    if found:
        break
    if i == numbers[-1]:
        print("No Possibilities Exist.")

It's getting stuck because the numbers you're checking are really large, and large numbers take a _really_ long time to factorize. — hyper-neutrino, Jun 15 '21 at 20:14
trial division is the simplest and weakest primality test algorithm there is, because is the slowest of them all for big numbers, you need a more advance algorithm such as Miller–Rabin or Baillie-PSW primality test — Copperfield, Jun 15 '21 at 20:57
your code get stuck at testing if 51452069511074530979 is prime or not, which is an absurdly big number to test with trial division (its square root is 7173009794) — Copperfield, Jun 15 '21 at 21:08

score 2 · Answer 1 · 2021-06-15T22:27:47.763

The problem is that your is_prime function is to slow, instead of checking if every number is a prime inside of your for loop. Why not generate a list of primes, lets say the first 1 million, store them in a list. Then too check if your number is prime, just check if it is inside of the list.

import math

def gen_prime(n):
    D = {}
    q = 2

    for i in range(n):
        if q not in D:
            yield q
            D[q * q] = [q]

        else:
            for p in D[q]:
                D.setdefault(p + q, []).append(p)
            del D[q]

        q += 1

primes = [i for i in gen_prime(1_000_000)]

def is_prime(num):
    if num in primes:
        return True
    else:
        return False

# create list containing a range of composite numbers
numbers = []
for i in range(4, 200):
    if not is_prime(i):
        numbers.append(i)

for i in numbers:
    found = False
    for j in numbers:
        if math.gcd(i, j) == 1:
            # print(i, "-", j, end=" | ")
            fibonacci = [i, j]
            contains_prime = False
            for k in range(2, 500):
                if is_prime(fibonacci[-2] + fibonacci[-1]):
                    contains_prime = True
                    break
                else:
                    fibonacci = [fibonacci[-1], fibonacci[-2] + fibonacci[-1]]
            if not contains_prime:
                print(i, j)
                found = True
        if found:
            break
    if found:
        break
    if i == numbers[-1]:
        print("No Possibilities Exist.")

when running this code I get the result

18 187
>>>

EDIT: I decided to do some research into prime algorithms, I came across Miller–Rabin.

The Miller–Rabin primality test or Rabin–Miller primality test is a probabilistic primality test: an algorithm which determines whether a given number is likely to be prime, similar to the Fermat primality test and the Solovay–Strassen primality test.

More information about it can be found on Wikipedia.

import random, math

# miller_rabin algorithm
def is_prime(n, k = 40):
    if n == 2:
        return True

    if n % 2 == 0:
        return False

    r, s = 0, n - 1
    while s % 2 == 0:
        r += 1
        s //= 2

    for i in range(k):
        a = random.randrange(2, n - 1)
        x = pow(a, s, n)
        if x == 1 or x == n - 1:
            continue

        for i in range(r - 1):
            x = pow(x, 2, n)
            if x == n - 1:
                break

        else:
            return False

    return True

# create list containing a range of composite numbers
numbers = []
for i in range(4, 200):
    if not is_prime(i):
        numbers.append(i)

for i in numbers:
    found = False
    for j in numbers:
        if math.gcd(i, j) == 1:
            # print(i, "-", j, end=" | ")
            fibonacci = [i, j]
            contains_prime = False
            for k in range(2, 500):
                if is_prime(fibonacci[-2] + fibonacci[-1]):
                    contains_prime = True
                    break
                else:
                    fibonacci = [fibonacci[-1], fibonacci[-2] + fibonacci[-1]]
            if not contains_prime:
                print(i, j)
                found = True
        if found:
            break
    if found:
        break
    if i == numbers[-1]:
        print("No Possibilities Exist.")

Here is the full code, as you can see when run it returns the result

143 142
>>>

the code eventually ask for the primality of number significantly bigger than `10**6` your prime testing function will give false positive/negative for those — Copperfield, Jun 15 '21 at 21:18
Just increase the the counter `primes = [i for i in gen_prime(1_000_000)]` instead of a million put a billion etc. `gen_prime(1_000_000_000)` or how ever big of a number you want to go up too. Storing the primes in RAM is a **much faster solution**, then looping through and individually checking each number using the OP's original code. — , Jun 15 '21 at 21:30
one of the number it test is 51452069511074530979, generate and store all the primes up to that number (`10**20`) to properly test it does not look like a good idea to me... — Copperfield, Jun 15 '21 at 21:41
@Jared Thank you Jared! Copperfield's solution runs much faster as I scale up the numbers, so I decided to accept his answer. — Daniel, Jun 16 '21 at 03:18

score 2 · Accepted Answer · answered Jun 15 '21 at 22:10

as I mention in the comments you need a better primality test algorithm, those algorithm can be hard to understand or implement on your own, so you can use a library that offer it instead like sympy for example, just install it with pip and use it as:

import math
from sympy.ntheory import isprime 

# create list containing a range of composite numbers
numbers = []
for i in range(4, 200):
    if not isprime(i):
        numbers.append(i)

for i in numbers:
    found = False
    for j in numbers:
        if math.gcd(i, j) == 1:
            fibonacci = [i, j]
            contains_prime = False
            for k in range(2, 500):
                if isprime(fibonacci[-2] + fibonacci[-1]):
                    contains_prime = True
                    break
                else:
                    fibonacci = [fibonacci[-1], fibonacci[-2] + fibonacci[-1]]
            if not contains_prime:
                print(i, j)
                found = True
        if found:
            break
    if found:
        break
    if i == numbers[-1]:
        print("No Possibilities Exist.")

now is quickly and run without getting stuck and print as result: 143 142

It worked flawlessly in my jupyter-notebook, but when I run it in my PyCharm IDE, it does not recognize sympy even though I have it installed. Do you know why this happens? — Daniel, Jun 16 '21 at 03:11
If those use their own separate python installation, that might be the cause, otherwise no idea — Copperfield, Jun 16 '21 at 03:25
Found a temporary fix [here](https://stackoverflow.com/questions/31235376/pycharm-doesnt-recognise-installed-module#:~:text=Solution%3A%20just%20install%20the%20package,you%20will%20see%20it%20pycharm%20.) — Daniel, Jun 16 '21 at 03:34

Fibonacci Sequence: Finding Composites Problem

2 Answers2