# Difference between revisions of "Prime numbers"

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

(Note the use of `takeWhile`

to prevent the infinite list of primes requiring an infinite amount of CPU time and RAM to process!)

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`

effectively implements a heap, exploiting Haskell's lazy evaluation model. For another example of this idiom see the Prelude's `ShowS`

type, which again exploits Haskell's lazy evaluation model
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
```