# Prime numbers

### From HaskellWiki

(Difference between revisions)

CaleGibbard (Talk | contribs) |
DonStewart (Talk | contribs) (bitwise prime sieve with template haskell) |
||

Line 44: | Line 44: | ||

<hask>nonprimes</hask> effectively implements a heap, exploiting Haskell's lazy evaluation model. For another example of this idiom see the Prelude's <hask>ShowS</hask> type, which again exploits Haskell's lazy evaluation model | <hask>nonprimes</hask> effectively implements a heap, exploiting Haskell's lazy evaluation model. For another example of this idiom see the Prelude's <hask>ShowS</hask> type, which again exploits Haskell's lazy evaluation model | ||

to avoid explicitly coding efficient concatenable strings. This is generalized by the [http://hackage.haskell.org/cgi-bin/hackage-scripts/package/dlist-0.3 DList package]. | to avoid explicitly coding efficient concatenable strings. This is generalized by the [http://hackage.haskell.org/cgi-bin/hackage-scripts/package/dlist-0.3 DList package]. | ||

+ | |||

+ | |||

+ | == Bitwise prime sieve == | ||

+ | |||

+ | Count the number of prime below a given 'n'. Shows fast bitwise arrays, | ||

+ | and an example of Template Haskell to defeat your enemies. | ||

+ | |||

+ | <haskell> | ||

+ | {-# OPTIONS -O2 -optc-O -fbang-patterns #-} | ||

+ | |||

+ | module Primes (pureSieve) where | ||

+ | |||

+ | import Control.Monad.ST | ||

+ | import Data.Array.ST | ||

+ | import Data.Array.Base | ||

+ | import System | ||

+ | import Control.Monad | ||

+ | import Data.Bits | ||

+ | |||

+ | pureSieve :: Int -> Int | ||

+ | pureSieve n = runST ( sieve n ) | ||

+ | |||

+ | sieve n = do | ||

+ | a <- newArray (3,n) True :: ST s (STUArray s Int Bool) | ||

+ | let cutoff = truncate (sqrt (fromIntegral n)) + 1 | ||

+ | go a n cutoff 3 1 | ||

+ | go !a !m cutoff !n !c | ||

+ | | n >= m = return c | ||

+ | | otherwise = do | ||

+ | e <- unsafeRead a n | ||

+ | if e then | ||

+ | if n < cutoff | ||

+ | then let loop !j | ||

+ | | j < m = do | ||

+ | x <- unsafeRead a j | ||

+ | when x $ unsafeWrite a j False | ||

+ | loop (j+n) | ||

+ | |||

+ | | otherwise = go a m cutoff (n+2) (c+1) | ||

+ | |||

+ | in loop ( if n < 46340 then n * n else n `shiftL` 1) | ||

+ | else go a m cutoff (n+2) (c+1) | ||

+ | |||

+ | else go a m cutoff (n+2) c | ||

+ | |||

+ | |||

+ | </haskell> | ||

+ | |||

+ | And places in a module: | ||

+ | |||

+ | <haskell> | ||

+ | {-# OPTIONS -fth #-} | ||

+ | |||

+ | import Primes | ||

+ | |||

+ | main = print $( let x = pureSieve 10000000 in [| x |] ) | ||

+ | </haskell> | ||

+ | |||

+ | Run as: | ||

+ | |||

+ | <haskell> | ||

+ | $ ghc --make -o primes Main.hs | ||

+ | $ time ./primes | ||

+ | 664579 | ||

+ | ./primes 0.00s user 0.01s system 228% cpu 0.003 total | ||

+ | </haskell> | ||

[[Category:Code]] | [[Category:Code]] |

## Revision as of 13:24, 15 July 2007

The following is an elegant (and highly inefficient) way to generate a list of all the prime numbers in the universe:

primes = sieve [2..] where sieve (p:xs) = p : sieve (filter (\x -> x `mod` p > 0) xs)

With this definition made, a few other useful (??) functions can be added:

is_prime n = n `elem` (takeWhile (n >=) primes) factors n = filter (\p -> n `mod` p == 0) primes factorise 1 = [] factorise n = let f = head $ factors n in f : factorise (n `div` f)

takeWhile

The following is a more efficient prime generator, implementing the sieve of Eratosthenes:

merge xs@(x:xt) ys@(y:yt) = case compare x y of LT -> x : (merge xt ys) EQ -> x : (merge xt yt) GT -> y : (merge xs yt) diff xs@(x:xt) ys@(y:yt) = case compare x y of LT -> x : (diff xt ys) EQ -> diff xt yt GT -> diff xs yt primes, nonprimes :: [Int] primes = [2,3,5] ++ (diff [7,9..] nonprimes) nonprimes = foldr1 f . map g $ tail primes where f (x:xt) ys = x : (merge xt ys) g p = [ n*p | n <- [p,p+2..]]

nonprimes

ShowS

to avoid explicitly coding efficient concatenable strings. This is generalized by the DList package.

## Bitwise prime sieve

Count the number of prime below a given 'n'. Shows fast bitwise arrays, and an example of Template Haskell to defeat your enemies.

{-# OPTIONS -O2 -optc-O -fbang-patterns #-} module Primes (pureSieve) where import Control.Monad.ST import Data.Array.ST import Data.Array.Base import System import Control.Monad import Data.Bits pureSieve :: Int -> Int pureSieve n = runST ( sieve n ) sieve n = do a <- newArray (3,n) True :: ST s (STUArray s Int Bool) let cutoff = truncate (sqrt (fromIntegral n)) + 1 go a n cutoff 3 1 go !a !m cutoff !n !c | n >= m = return c | otherwise = do e <- unsafeRead a n if e then if n < cutoff then let loop !j | j < m = do x <- unsafeRead a j when x $ unsafeWrite a j False loop (j+n) | otherwise = go a m cutoff (n+2) (c+1) in loop ( if n < 46340 then n * n else n `shiftL` 1) else go a m cutoff (n+2) (c+1) else go a m cutoff (n+2) c

And places in a module:

{-# OPTIONS -fth #-} import Primes main = print $( let x = pureSieve 10000000 in [| x |] )

Run as:

$ ghc --make -o primes Main.hs $ time ./primes 664579 ./primes 0.00s user 0.01s system 228% cpu 0.003 total