Fibonacci primes in parallel

module Main where
 
import Data.Time.Clock.POSIX
import Data.Maybe
import Numeric
import Char
import Control.Parallel.Strategies
 
-- binary representation of an integer
binary :: Integer -> String 
binary = flip (showIntAtBase 2 intToDigit) []
 
-- returns True if (a) is a witness that odd n is compound
witness :: Integer -> Integer -> Bool
witness n a = pow (tail $ binary $ n-1) a
  where
    pow _ 1 = False
    pow [] _ = True
    pow xs d = pow' xs d $ (d*d) `mod` n
      where
        pow' _ d d2 | d2 == 1 && d /= (n-1) = True
        pow' ('0':xs) d d2 = pow xs d2
        pow' ('1':xs) d d2 = pow xs $ (d2*a) `mod` n
 
-- is (n) a prime with probability 4^(-k)
miller_rabin k n
  | even n    = n == 2
  | otherwise = not $ any (witness n) $ takeWhile (< n) $ take k primes
 
primes :: [Integer]
primes = 2:3:5: [ x | x <- candidates 7 11, isPrime x ]
    where
      candidates a1 a2 = a1 : a2 : candidates (a1+6) (a2+6)
 
-- simple primality test applied for indices only, not for Fibonacci numbers      
isPrime x = all ((0 /=) . (x `mod`)) $ takeWhile ((x>=).(^2)) primes
 
fib = 1 : 1 : [ a+b | (a,b) <- zip fib (tail fib) ]
 
-- indexed Fibonacci numbers
numfib = zip [1..] fib
 
isPrimeOr4 x = x /= 1 && x /=2 && (x == 4 || isPrime x)
 
-- Fibonacci numbers with primal indices
numfib' = filter (isPrimeOr4 . fst) numfib
 
isProbablyPrime = miller_rabin 10
 
maybeFibprimes = [if isProbablyPrime f then Just (n,f) else Nothing | (n,f) <- numfib' ]
 
-- probably you need to increase 10 if running the program
-- in more than two threads
fibprimes = catMaybes $ parBuffer 10 rnf maybeFibprimes
 
main = do
  start <- getPOSIXTime
  printEach fibprimes start 1
    where printEach (x:xs) start n = do
          t <- getPOSIXTime
          print (t - start, n, fst x)
          printEach xs start (n+1)

ghc --make –threaded parfibs.hs

parfibs +RTS -N2

(0s,1,3)
(0s,2,4)
(0s,3,5)
(0s,4,7)
(0s,5,11)
(0s,6,13)
(0s,7,17)
(0s,8,23)
(0s,9,29)
(0s,10,43)
(0s,11,47)
(0s,12,83)
(0.015625s,13,131)
(0.015625s,14,137)
(0.03125s,15,359)
(0.046875s,16,431)
(0.046875s,17,433)
(0.046875s,18,449)
(0.0625s,19,509)
(0.078125s,20,569)
(0.078125s,21,571)
(2.546875s,22,2971)
(10.890625s,23,4723)
(16.859375s,24,5387)
(104.390625s,25,9311)
(120.546875s,26,9677)
(464.359375s,27,14431)
(3368.578125s,28,25561)
(6501.296875s,29,30757)
(11361.453125s,30,35999)
(13182.71875s,31,37511)
(37853.765625s,32,50833)