Prime numbers, Miller-Rabin

EDIT 2011 March: This implementation is buggy. Do not use it.

I was looking at a few interesting algorithms to figure out whether a given number was a prime, all because I wanted the fastest solutions possible while solving project euler sums involving prime numbers. And after that session, the sum is still not over and might never be as I have gotten around to post on this blog, the stuff I’ve gathered.The AKS algorithm is of extreme academic importance but has horrible speeds and hence is of no use to me whose main aim to get prime number with utter ease.

I then looked at the Miller-Rabin algorithm and then moved to try out an implementation on my own. I cheated by trying to look at implementations on the cookbook and understood nothing. So it was back to me to figure out something. The algorithm goes like this:

To check if “n” is prime, “n” should obviously be 2 or odd. Forget 2, because if you don’t, I will call you names ( and who needs miller-rabin to prove if 2 is prime ? ).

Since n is odd, n-1 is even. Try expressing n – 1 as {2}^{s}.d. It seems relatively simple. Here is the code to do that:

#!/usr/bin/python

#Miller-Rabin Primality Test.

#Author: Shriphani Palakodety a.k.a PSP#Now comes one of the world’s fastest primality tests:

def twoPower(n):

 factors = []

 s = 0

 while n > 0:

        if n % 2 == 0:

                factors.append(0)

                n = n / 2

                s += 1

        else:

                factors.append(n)

                d = n

                break

 return s, d

Once that is done, we need to generate a random number which lies in the range \left[1, n - 1\right]. Here is the code to do that:

#!/usr/bin/pythonfrom random import randint

def randGen(n):
 return randint(1, n)

Ok, I agree there was no need for that snippet. Let us move on anyway. The crux of this algorithm lies in the following statement:

if a^d \not\equiv 1\ \bmod\ n and a^{{2^r}d}\ \not\equiv -1 \bmod\ n for all r in the range [0, s - 1\ ] then return composite.

Well, this is quite easy to implement in Python:
UPDATED FUNCTION

#!/usr/bin/env python

#Modulo Check

def moduloCheck(rand_int, odd_no, s, n):
        a = True
        if (rand_int ** odd_no – 1) % n == 0:
                a = False
        for i in xrange(0, s):
                if (rand_int ** ((2**i) * odd_no) + 1) % n==0:
                        a = False
        if a:
                return "Composite"
        else:
                return "Prime"

Now comes the trouble:

number = 73n = twoPower(number- 1)[0]

odd_no = twoPower(number – 1)[1]

rand_int = randGen(n – 1)

print moduloCheck(rand_int, odd_no, n)

Have a look at what I get:

shriphani@psp-laptop:~/project_euler$ python miller_rabin.py

Prime

shriphani@psp-laptop:~/project_euler$ python miller_rabin.py

Prime

shriphani@psp-laptop:~/project_euler$ python miller_rabin.py

Prime

shriphani@psp-laptop:~/project_euler$ python miller_rabin.py

Composite

BOOM !

The trick to avoiding such problems is to just run the moduloCheck function as many times as possible and if we get more Primes than Composites, then the number is Prime, else Composite :)

I then decided to create a parameter “k” which would determine the number of times the moduloCheck function should run. Then I implemented a counter to check the frequency of the strings – “Prime” and “Composite” thrown out by the script. Here is the code:
UPDATED FUNCTION

k = 31

def freqCheck(k, odd_no, n):
        freq_list = []
        for i in range(0, k):
                freq_list.append(moduloCheck(randGen(number-1), odd_no, n))
        if freq_list.count("Prime") > freq_list.count("Composite"):
                return "Prime"
        else:
                return "Composite"

So, the entire script looks like this: http://shriphani.com/scripts/miller-rabin.py.html

About this entry

You’re currently reading “Prime numbers, Miller-Rabin,” an entry on Shriphani Palakodety

Published:
04.09.08 / 12pm
Category:
Mathematics, python