CODEWITHSUNDEEP
pythoncryptographyprimes
thelastone
thelastone

Reputation: 575

Probability of finding a prime (using miller-rabin test)

I've implemented Miller-Rabin primality test and every function seems to be working properly in isolation. However, when I try to find a prime by generating random numbers of 70 bits my program generates in average more than 100000 numbers before finding a number that passes the Miller-Rabin test (10 steps). This is very strange, the probability of being prime for a random odd number of less than 70 bits should be very high (more than 1/50 according to Hadamard-de la Vallée Poussin Theorem). What could be wrong with my code? Would it be possible that the random number generator throws prime numbers with very low probability? I guess not... Any help is very welcome.

import random


def miller_rabin_rounds(n, t):
    '''Runs miller-rabin primallity test t times for n'''

    #  First find the values r and s such that 2^s * r = n - 1
    r = (n - 1) / 2
    s = 1
    while r % 2 == 0:
        s += 1
        r /= 2

    #  Run the test t times
    for i in range(t):
        a = random.randint(2, n - 1)
        y = power_remainder(a, r, n)

        if y != 1 and y != n - 1:
            #  check there is no j for which (a^r)^(2^j) = -1 (mod n)
            j = 0
            while j < s - 1 and y != n - 1:
                y = (y * y) % n
                if y == 1:
                    return False
                j += 1
            if y != n - 1:
                return False

    return True


def power_remainder(a, k, n):
    '''Computes (a^k) mod n efficiently by decomposing k into binary'''
    r = 1
    while k > 0:
        if k % 2 != 0:
            r = (r * a) % n
        a = (a * a) % n
        k //= 2
    return r


def random_odd(n):
    '''Generates a random odd number of max n bits'''
    a = random.getrandbits(n)
    if a % 2 == 0:
        a -= 1
    return a

if __name__ == '__main__':
    t = 10  # Number of Miller-Rabin tests per number
    bits = 70  # Number of bits of the random number
    a = random_odd(bits)
    count = 0
    while not miller_rabin_rounds(a, t):
        count += 1
        if count % 10000 == 0:
            print(count)
        a = random_odd(bits)
    print(a)

Upvotes: 2

Views: 1350

Answers (1)

nbryans
nbryans

Reputation: 1517

The reason this works in python 2 and not python 3 is that the two handle integer division differently. In python 2, 3/2 = 1, whereas in python 3, 3/2=1.5.

It looks like you should be forcing integer division in python 3 (rather than float division). If you change the code to force integer division (//) as such:

#  First find the values r and s such that 2^s * r = n - 1
r = (n - 1) // 2
s = 1
while r % 2 == 0:
    s += 1
    r //= 2

You should see the correct behaviour regardless of what python version you use.

Upvotes: 3

Related Questions