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 . It seems relatively simple. Here is the code to do that:
#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 . Here is the code to do that:
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 
and 
for all r in the range 
] then return composite.
Well, this is quite easy to implement in Python:
UPDATED FUNCTION
#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:
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
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
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.
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.
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.
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.
[...] (which helped me a lot, thanks). And then I mixed it up with my miller rabin implementation and had a bit of fun [...]
This test is for testing whether a number is composite or not.If a number is composite it will always return “Composite”.Even if the number of “Composite” are more than “Prime” if there is at least one prime the number is prime,not composite.