My Weblog

Blog about programming and math

Faster Rabin Miller Test in Haskell

I wrote this program to learn IO in Haskell and its bit faster then mine previous implementation [Technically that implementation is not fully owned by me. A lot of code was borrowed from Haskell site]. The Rabin Miller test can be written as pure function i but i rather prefer impurity :P. You can use filterM isProbable [k1..k2] to get primes in the range [k1,k2]. Remember filterM is not lazy. You can also have a look here.

import Random

helpMod a d n= fun a d where 
	fun a 1 = mod a n
	fun a 2 = mod (a^2) n
	fun a d = let 
			p=fun (mod (a^2) n) (div d 2) 
			q=if odd d then mod (a*p) n else p
		  in mod q n

isProbable::Integer->IO Bool
isProbable n |n<=1 = return False
	     |n==2 = return True
	     |even n = return False
	     |otherwise	= 
		 let d=until (\x->mod x 2/=0) (`div` 2) (n-1) --(n-1=2^s*d)
		     s=until (\x->d*2^x==pred n) (+1) 0  --counts the power
	         rabinMiller 0 n s d

rabinMiller::Integer->Integer->Integer->Integer->IO Bool
rabinMiller cnt n s d | cnt>=5= return True
		      | otherwise = 
	 a<-getStdRandom $ randomR (2,n-2)
	 let x=helpMod a d n
	 --let l="s= "++ show s ++" d= "++ show d ++" a= "++ show a ++" x= "++ show x
	 --putStrLn l
	 case x==1 of
	   True -> rabinMiller (cnt+1) n s d
	   _ -> if any ( ==pred n) $ take (fromIntegral s) $ iterate (\e->mod (e^2) n)  x then rabinMiller (cnt+1) n s d 
		else return False

main = do 
	 let n=read l
	 k<-isProbable n
	 print k

February 5, 2011 Posted by | Programming | 1 Comment