Prime numbers

Simple Prime sieve

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

  primes :: [Integer]
  primes = sieve [2..]
    where sieve (p:xs) = p : sieve [x | x<-xs, x `mod` p /= 0]

Given an infinite list of prime numbers, we can implement primality tests and integer factorization:

  isPrime n = n > 1 && n == head (factorize n)

  factorize 1 = []
  factorize n = go primes n
     where
     go ps@(p:pt) n
        | p*p > n        = [n]
        | x `mod` p == 0 = p : go (n `div` p) ps
        | otherwise      = go n pt

Prime Wheels

TODO

Implicit Heap

The following is a more efficient prime generator, implementing the sieve of Eratosthenes. See also the message thread around Re: Code and Perf. Data for Prime Finders.

data People a = VIP a (People a) | Crowd [a]

mergeP :: Ord a => People a -> People a -> People a
mergeP (VIP x xt) ys                    = VIP x $ mergeP xt ys
mergeP (Crowd xs) (Crowd ys)            = Crowd $ merge  xs ys
mergeP xs@(Crowd ~(x:xt)) ys@(VIP y yt) = case compare x y of
    LT -> VIP x $ mergeP (Crowd xt) ys
    EQ -> VIP x $ mergeP (Crowd xt) yt
    GT -> VIP y $ mergeP xs yt


merge :: Ord a => [a] -> [a] -> [a]
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

foldTree :: (a -> a -> a) -> [a] -> a
foldTree f ~(x:xs) = f x . foldTree f . pairs $ xs
    where pairs ~(x: ~(y:ys)) = f x y : pairs ys

primes, nonprimes :: [Integer]
primes    = 2:3:diff [5,7..] nonprimes
nonprimes = serve . foldTree mergeP . map multiples $ tail primes
    where
    multiples p = vip [p*k | k <- [p,p+2..]]

    vip (x:xs)       = VIP x $ Crowd xs
    serve (VIP x xs) = x:serve xs
    serve (Crowd xs) = xs

nonprimes effectively implements a heap, exploiting lazy evaluation.

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

Yet another set of prime functions

I'm sure this isn't the most efficient, but it is good enough for many applications and it is simple enough for a beginning Haskeller like me to create. It only checks prime numbers as factors when trying to factor a new candidate and only checks odd candidates, otherwise it is pretty straightforward. Working more on efficiency wasn't worth it for the program that used it. Uses integral instead of int to allow for bigger numbers. Thanks to this, I now know that none of my (10-digit) phone numbers are prime numbers. Now to find the factors of my phone numbers...

--| Simple prime functions
--| Produces an infinite list of primes.
myprimelist :: Integral a => [a]
myprimelist = 2 : [x | x <- [3,5..], myprimecheck x (takeWhile (<x) myprimelist)]

--| simple nth prime using the list
mynthprime :: (Integral a, Integral b) => a -> b
mynthprime 1 = 2
mynthprime n = head $ drop  (fromIntegral (n-1)) myprimelist

--| Check a number for primeness.
myisprime :: Integral a => a -> Bool
myisprime n | n<2 = False
            | n<4 = True
            | even n = False
            | otherwise = myprimecheck n (takeWhile (<n) myprimelist)

--| Used by myisprime and myprimelist to check for primeness.
myprimecheck :: Integral a => a -> [a] -> Bool
myprimecheck n [] = True
myprimecheck n (f:fs)	= if f*f <= n 
                                then
                                        if (n `mod` f)==0 
                                                then False
                                                else myprimecheck n fs
                                else True

Miller-Rabin Primality Test

find2km :: Integral a => a -> (a,a)
find2km n = f 0 n
    where 
        f k m
            | r == 1 = (k,m)
            | otherwise = f (k+1) q
            where (q,r) = quotRem m 2        
 
millerRabinPrimality :: Integer -> Integer -> Bool
millerRabinPrimality n a
    | a <= 1 || a >= n-1 = 
        error $ "millerRabinPrimality: a out of range (" 
              ++ show a ++ " for "++ show n ++ ")" 
    | n < 2 = False
    | even n = False
    | b0 == 1 || b0 == n' = True
    | otherwise = iter (tail b)
    where
        n' = n-1
        (k,m) = find2km n'
        b0 = powMod n a m
        b = take (fromIntegral k) $ iterate (squareMod n) b0
        iter [] = False
        iter (x:xs)
            | x == 1 = False
            | x == n' = True
            | otherwise = iter xs
 
pow' :: (Num a, Integral b) => (a -> a -> a) -> (a -> a) -> a -> b -> a
pow' _ _ _ 0 = 1
pow' mul sq x' n' = f x' n' 1
    where 
        f x n y
            | n == 1 = x `mul` y
            | r == 0 = f x2 q y
            | otherwise = f x2 q (x `mul` y)
            where
                (q,r) = quotRem n 2
                x2 = sq x
 
mulMod :: Integral a => a -> a -> a -> a
mulMod a b c = (b * c) `mod` a
squareMod :: Integral a => a -> a -> a
squareMod a b = (b * b) `rem` a
powMod :: Integral a => a -> a -> a -> a
powMod m = pow' (mulMod m) (squareMod m)