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