Taqi H
Taqi H

Reputation: 27

Why is my implementation of the Miller-Rabin algorithm not able to detect some primes?

I am trying to implement the Miller-Rabin primality checker for some project am working. However, the algorithm doesn't work for primes such as 101, 103, 107, 109... I can't figure out where the problem is. Thanks in advance for all the help.

def miller_rabin_is_prime(number, k=10):

    if number < 2:
        return False
    elif number <= 3:
        return True

    else:
        odd_num, power_of_two, factor_out = 0, 0, number - 1

        while number != (2 ** power_of_two)*odd_num + 1:
            if factor_out / 2 == int(factor_out / 2):
                power_of_two += 1
                factor_out /= 2
            else:
                odd_num = (number - 1) / (2 ** power_of_two)

        for _ in range(k):
            random = randint(2, number - 2)
            checker = (random**odd_num) % number

            if (checker == 1) or (checker == number - 1):
                continue
            try:
                for loop in range(power_of_two - 1):
                    checker = (checker**2) % number
                    if checker == number - 1:
                        raise TypeError
            except TypeError:
                continue
            return False

        return True

I expect the output for 101 to be True, but the actual output is False.

Upvotes: 0

Views: 174

Answers (1)

John Coleman
John Coleman

Reputation: 52008

If you replace

 odd_num = (number - 1) / (2 ** power_of_two)

by

 odd_num = (number - 1) // (2 ** power_of_two)

your code will work -- but fairly slowly for larger numbers. To improve the code:

  1. Use a simpler method of computing odd_num and power_of_two
  2. Use pow() for the modular exponentiation.

Something like:

from random import randint
def miller_rabin_is_prime(number, k=10):

    if number < 2:
        return False
    elif number <= 3:
        return True

    else:
        odd_num = number - 1
        power_of_two = 0

        while odd_num % 2 == 0:
            power_of_two += 1
            odd_num //= 2

        for _ in range(k):
            random = randint(2, number - 2)
            checker = pow(random,odd_num, number)

            if (checker == 1) or (checker == number - 1):
                continue
            try:
                for loop in range(power_of_two - 1):
                    checker = pow(checker,2,number)
                    if checker == number - 1:
                        raise TypeError
            except TypeError:
                continue
            return False

        return True

Then, for example, miller_rabin_is_prime(1000003) will evaluate to True almost instantly, whereas your original code (even after / is replaced by //) would take about 15 seconds because of the non-modular exponentiation.

As a final remark, you are using error handling for non-error conditions (clearly there is no type error when checker == number - 1). It would be much cleaner to refactor your main loop so that it doesn't use try--except. Error handling is not meant for ordinary control flow.

Upvotes: 1

Related Questions