Weblog of an Aspiring Computer Scientist
Random header image... Refresh for more!

Prime numbers, Miller-Rabin

by Shriphani

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:

#!/usr/bin/env python

#Modulo Check
def moduloCheck(rand_int, odd_no, n):

 if rand_int ** odd_no == 1 % n:

        for i in range(0, n):

                if rand_int ** ( ( 2**i) * odd_no ) == -1 % n:

                        continue

                else:

                        break

                        return "Prime"

 else:

        return "Prime"

 return "Composite"

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:

k = randint(1, n)

def freqCheck(k, odd_no, n):
        freq_list = []
        for i in range(0, k):
                freq_list.append(moduloCheck(randGen(n - 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.html

0 comments

There are no comments yet...

Kick things off by filling out the form below.

Leave a Comment