## Programming With Prime Numbers: Source Code In Haskell

```-- Programming With Prime Numbers -- http://programmingpraxis.files.wordpress.com/ -- 2012/09/programmingwithprimenumbers.pdf```

``` import Control.Monad (forM_, when) import Control.Monad.ST import Data.Array.ST import Data.Array.Unboxed import Data.List (sort) ancientSieve :: Int -> UArray Int Bool ancientSieve n = runSTUArray \$ do     bits <- newArray (2, n) True     forM_ [2 .. n] \$ \p -> do         isPrime <- readArray bits p         when isPrime \$ do             forM_ [2*p, 3*p .. n] \$ \i -> do                 writeArray bits i False     return bits ancientPrimes :: Int -> [Int] ancientPrimes n = [p | (p, True) <-                   assocs \$ ancientSieve n] sieve :: Int -> UArray Int Bool sieve n = runSTUArray \$ do     let m = (n-1) `div` 2         r = floor . sqrt \$ fromIntegral n     bits <- newArray (0, m-1) True     forM_ [0 .. r `div` 2 - 1] \$ \i -> do         isPrime <- readArray bits i         when isPrime \$ do             let a = 2*i*i + 6*i + 3                 b = 2*i*i + 8*i + 6             forM_ [a, b .. (m-1)] \$ \j -> do                 writeArray bits j False     return bits primes :: Int -> [Int] primes n = 2 : [2*i+3 | (i, True) <- assocs \$ sieve n] tdPrime :: Int -> Bool tdPrime n = prime (2:[3,5..])     where prime (d:ds)             | n < d * d = True             | n `mod` d == 0 = False             | otherwise = prime ds tdFactors :: Int -> [Int] tdFactors n = facts n (2:[3,5..])     where facts n (f:fs)             | n < f * f = [n]             | n `mod` f == 0 =                   f : facts (n `div` f) (f:fs)             | otherwise = facts n fs powmod :: Integer -> Integer -> Integer -> Integer powmod b e m =     let times p q = (p*q) `mod` m         pow b e x             | e == 0 = x             | even e = pow (times b b)                               (e `div` 2) x             | otherwise = pow (times b b)                               (e `div` 2) (times b x)     in pow b e 1 isSpsp :: Integer -> Integer -> Bool isSpsp n a =     let getDandS d s =             if even d then getDandS (d `div` 2) (s+1)                       else (d, s)         spsp (d, s) =             let t = powmod a d n             in if t == 1 then True else doSpsp t s         doSpsp t s             | s == 0 = False             | t == (n-1) = True             | otherwise = doSpsp ((t*t) `mod` n) (s-1)     in spsp \$ getDandS (n-1) 0 isPrime :: Integer -> Bool isPrime n =     let ps = [2,3,5,7,11,13,17,19,23,29,31,37,41,              43,47,53,59,61,67,71,73,79,83,89,97]     in n `elem` ps || all (isSpsp n) ps rhoFactor :: Integer -> Integer -> Integer rhoFactor n c =     let f x = (x*x+c) `mod` n         fact t h             | d == 1 = fact t' h'             | d == n = rhoFactor n (c+1)             | isPrime d = d             | otherwise = rhoFactor d (c+1)             where                 t' = f t                 h' = f (f h)                 d = gcd (t' - h') n     in fact 2 2 rhoFactors :: Integer -> [Integer] rhoFactors n =     let facts n             | n == 2 = [2]             | even n = 2 : facts (n `div` 2)             | isPrime n = [n]             | otherwise = let f = rhoFactor n 1                           in f : facts (n `div` f)     in sort \$ facts n ```

```main = do     print \$ ancientPrimes 100     print \$ primes 100     print \$ length \$ primes 1000000     print \$ tdPrime 716151937     print \$ tdFactors 8051     print \$ powmod 437 13 1741     print \$ isSpsp 2047 2     print \$ isPrime 600851475143     print \$ isPrime 2305843009213693951     print \$ rhoFactors 600851475143```