Prime numbers

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