Personal tools

Fibonacci primes in parallel

From HaskellWiki

(Difference between revisions)
Jump to: navigation, search
m
 
(5 intermediate revisions by 2 users not shown)
Line 1: Line 1:
 
The problem is to find all primes in the sequence of rapidly growing Fibonacci numbers.
 
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.
+
[http://en.wikipedia.org/wiki/Miller-Rabin_primality_test 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.
+
The [http://www.haskell.org/ghc/dist/current/docs/users_guide/lang-parallel.html#id437410 parallel capabilities] of [http://haskell.org/ghc GHC] give almost twofold performance speedup if running the program on a dual core system.
   
 
<haskell>
 
<haskell>
Line 32: Line 32:
 
-- is (n) a prime with probability 4^(-k)
 
-- is (n) a prime with probability 4^(-k)
 
miller_rabin k n
 
miller_rabin k n
| n `mod` 2 == 0 = n == 2
+
| even n = n == 2
 
| otherwise = not $ any (witness n) $ takeWhile (< n) $ take k primes
 
| otherwise = not $ any (witness n) $ takeWhile (< n) $ take k primes
   
Line 119: Line 119:
 
== Optimising further ==
 
== 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, we can find the 25th prime in:
+
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
 
$ ghc -O2 -optc-O2 -fvia-C -threaded A.hs -no-recomp --make
Line 127: Line 127:
 
$ time ./A +RTS -N4
 
$ time ./A +RTS -N4
 
....
 
....
(79.675504s,25,9311)
+
(69.869631s, 25, 9311)
(89.850673s,26,9677)
+
(88.357936s, 26, 9677)
A: interrupted
+
(297.578843s, 27, 14431)
./A +RTS -N4 499.81s user 0.76s system 291% cpu 2:51.54 total
 
   
Note only 291% cpu utilisation, as we're contesting for cpu time on this machine.
+
Note that this program runs at around 300% cpu utilisation, despite having 4 cores, as we're contesting for cpu time on this machine.
  +
  +
== Improving the algorithm ==
  +
  +
Are there other [http://home.att.net/~blair.kelly/mathematics/fibonacci/factthms.html factorization theorems] for fibonacci primes we can exploit?
   
 
[[Category:Code]]
 
[[Category:Code]]

Latest revision as of 20:32, 14 December 2009

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.

The 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
  | 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)

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.

[edit] 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.

[edit] 2 Improving the algorithm

Are there other factorization theorems for fibonacci primes we can exploit?