Prime numbers, Miller-Rabin

April 9th, 2008 · 5 Comments · Mathematics, python

EDIT 2 June 2009: Upon a mistake being pointed out to me by Justin, I changed a few lines of the code. This should work but please do tell me if it still breaks for any case.

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

Tags:

5 Comments so far ↓

Justin // Jun 2, 2009 at 1:39 am

Just a note for anyone else who stumbles on to this code snippet… this isn’t a good implementation. The rand_int generated is less than n (or s, according to 2^s*d), whereas it needs to be less than the original number.

This implementation is not a good implementation, as 221 for example, is reported as prime more often than composite. Accuracy would be around 50% or less, whereas a good rabin-miller implementation would have 99.999% accuracy.
Shriphani // Jun 2, 2009 at 4:20 pm

Thanks for pointing out that this is wrong. rand_int’s bounds are correct but I also felt that the core of the test could also do with some refactoring since I believe that was also not implemented correctly. Anyway, please let me know if the new snippet is still failing. I will be happy to fix it.
Matt // Jun 11, 2009 at 3:52 pm

Maybe you can help me, I’m trying to implement Miller-Rabin as described in SICP, and I’m getting a little lost. The algorithm uses ‘1 % n’, but I thought this would always be equal to 1? At least that’s what my calculator says. My math is a little rusty and I know I’m being dense here, but Google isn’t helping me.

BTW, if you’re looking for the fastest speed possible you may want to look at SICP’s expmod function (in section 1.2.6), they show a method for calculating a ^ b % n without actually calculating a ^ b. Don’t know if it’s any faster in practice, but worth a shot.
Shriphani // Jun 11, 2009 at 11:43 pm

Matt,

The algorithm uses a check if a^p is congruent to 1 % n. a^p congruent to 1 % n means a^p -1 is divisible by n.

Nice to know SICP has expmod. It uses the idea that a^b % n = ((a % n)^b) % n. I did write it some time back and it is on this blog. Thanks for pointing that out, I will try to integrate modular exponentiation into this script.
Bio-Informatics — Shriphani ‘PSP’ Palakodety // Sep 4, 2009 at 5:20 pm

[...] (which helped me a lot, thanks). And then I mixed it up with my miller rabin implementation and had a bit of fun [...]

Shriphani ‘PSP’ Palakodety

Weblog of an Aspiring Computer Scientist