Difference between revisions of "Prime numbers"
Lannyripple (talk | contribs) (prime factorization and primality testing routines) |
(Fixed implicit heap implementation) |
||
Line 1: | Line 1: | ||
+ | == Simple Prime sieve == |
||
The following is an elegant (and highly inefficient) way to generate a list of all the prime numbers in the universe: |
The following is an elegant (and highly inefficient) way to generate a list of all the prime numbers in the universe: |
||
<haskell> |
<haskell> |
||
− | primes |
+ | primes :: [Integer] |
+ | primes = sieve [2..] |
||
− | sieve (p:xs) = p : sieve |
+ | where sieve (p:xs) = p : sieve [x | x<-xs, x `mod` p /= 0] |
</haskell> |
</haskell> |
||
Line 9: | Line 11: | ||
<haskell> |
<haskell> |
||
− | + | isPrime n = n `elem` (takeWhile (n >=) primes) |
|
factors n = filter (\p -> n `mod` p == 0) primes |
factors n = filter (\p -> n `mod` p == 0) primes |
||
Line 20: | Line 22: | ||
(Note the use of <hask>takeWhile</hask> to prevent the infinite list of primes requiring an infinite amount of CPU time and RAM to process!) |
(Note the use of <hask>takeWhile</hask> to prevent the infinite list of primes requiring an infinite amount of CPU time and RAM to process!) |
||
+ | |||
+ | == Implicit Heap == |
||
The following is a more efficient prime generator, implementing the sieve of |
The following is a more efficient prime generator, implementing the sieve of |
||
+ | Eratosthenes. See also the message thread around [http://thread.gmane.org/gmane.comp.lang.haskell.cafe/26426/focus=26493 Re: Code and Perf. Data for Prime Finders]. |
||
− | Eratosthenes: |
||
<haskell> |
<haskell> |
||
+ | 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 |
merge xs@(x:xt) ys@(y:yt) = case compare x y of |
||
− | LT -> x : |
+ | LT -> x : merge xt ys |
− | EQ -> x : |
+ | EQ -> x : merge xt yt |
− | GT -> y : |
+ | GT -> y : merge xs yt |
+ | |||
⚫ | |||
− | diff |
+ | diff xs@(x:xt) ys@(y:yt) = case compare x y of |
− | LT -> x : |
+ | LT -> x : diff xt ys |
− | EQ -> diff xt yt |
+ | EQ -> diff xt yt |
− | GT -> diff xs yt |
+ | GT -> diff xs yt |
+ | |||
+ | foldTree :: (a -> a -> a) -> [a] -> a |
||
+ | foldTree f ~(x:xs) = f x . foldTree f . pairs $ xs |
||
⚫ | |||
+ | |||
⚫ | |||
⚫ | |||
⚫ | |||
⚫ | |||
⚫ | |||
+ | vip (x:xs) = VIP x $ Crowd xs |
||
⚫ | |||
+ | serve (VIP x xs) = x:serve xs |
||
⚫ | |||
+ | serve (Crowd xs) = xs |
||
⚫ | |||
⚫ | |||
⚫ | |||
</haskell> |
</haskell> |
||
− | <hask>nonprimes</hask> effectively implements a heap, exploiting Haskell's lazy evaluation model. |
+ | <hask>nonprimes</hask> effectively implements a heap, exploiting 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]. |
||
Building on the prime generator just above, here's two functions for fast prime factorization and primality testing. |
Building on the prime generator just above, here's two functions for fast prime factorization and primality testing. |
||
Line 164: | Line 187: | ||
else True |
else True |
||
</haskell> |
</haskell> |
||
+ | |||
− | == Yet another test prime functions == |
||
+ | == Miller-Rabin Primality Test == |
||
<haskell> |
<haskell> |
||
find2km :: Integral a => a -> (a,a) |
find2km :: Integral a => a -> (a,a) |
Revision as of 12:54, 8 June 2008
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]
With this definition made, a few other useful (??) functions can be added:
isPrime 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!)
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 Haskell's lazy evaluation model.
Building on the prime generator just above, here's two functions for fast prime factorization and primality testing.
isPrime :: (Integral a) => a -> Bool
isPrime x
| x < 2 = False
| otherwise = x == head (ftoa x)
-- ftoa: fundamental theorem of arithmetic
ftoa :: (Integral a) => a -> [a]
ftoa x = ftoa' x primes
where
ftoa' x ps@(p:pt)
| x == 1 = []
| p * p > x = [x]
| mod x p == 0 = p : ftoa' (div x p) ps
| otherwise = ftoa' x pt
-- maybe primes isn't infinite?
ftoa' x [] = [x]
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)