Fibonacci primes in parallel
The problem is to find all primes in the sequence of rapidly growing Fibonacci numbers.
Miller-Rabin probabilistic primality test is used to check numbers in the sequence. In addition only the elements with prime indices in the sequence are considered due to known properties of Fibonacci numbers.
Parallel capabilities of ghc give almost twofold performance speedup if running the program on a dual core system.
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 | n `mod` 2 == 0 = 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)
In order to gain parallel benefits compile the program using
ghc --make –threaded parfibs.hs
then run by
parfibs +RTS -N2
For each found Fibonacci prime the program prints the time spent from the start. The following output was received on dual-core Athlon64 x2 4600+ microprocessor.
(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)
In about 10 hours the program found 32 Fibonacci primes from about 40 known so far. Of course, obtaining each next number takes more time than all preceding ones.
1 Optimising further
We can improve the performance somewhat by throwing more cores at the problem, and turning up the optimisation level. Compiling, with optimisations on, on a 4 core machine, also, replacing lazy pairs with strict ones, we get:
$ ghc -O2 -optc-O2 -fvia-C -threaded A.hs -no-recomp --make
We get the result slightly sooner, and with better cpu utilisation,
$ time ./A +RTS -N4 .... (69.869631s, 25, 9311) (88.357936s, 26, 9677) (297.578843s, 27, 14431)
Note that this program runs at around 300% cpu utilisation, despite having 4 cores, as we're contesting for cpu time on this machine.
2 Improving the algorithm
Are there other factorization threoems for fibonacci primes we can exploit?