Personal tools

99 questions/Solutions/31

From HaskellWiki

< 99 questions | Solutions(Difference between revisions)
Jump to: navigation, search
(Faster alternative)
 
(16 intermediate revisions by 5 users not shown)
Line 1: Line 1:
 
(**) Determine whether a given integer number is prime.
 
(**) Determine whether a given integer number is prime.
  +
  +
Well, a natural number ''k'' is a prime number if it is larger than '''1''' and no natural number ''n >= 2'' with ''n^2 <= k'' is a divisor of ''k''. However, we don't actually need to check all natural numbers ''n <= sqrt k''. We need only check the '''''primes''' p <= sqrt k'':
   
 
<haskell>
 
<haskell>
 
isPrime :: Integral a => a -> Bool
 
isPrime :: Integral a => a -> Bool
isPrime p = p > 1 && (all (\n -> p `mod` n /= 0 ) $ takeWhile (\n -> n*n <= p) [2..])
+
isPrime k = k > 1 &&
  +
foldr (\p r -> p*p > k || k `rem` p /= 0 && r)
  +
True primesTME
 
</haskell>
 
</haskell>
   
Well, a natural number p is a prime number iff it is larger than 1 and no natural number n with n >= 2 and n^2 <= p is a divisor of p. That's exactly what is implemented: we take the list of all integral numbers starting with 2 as long as their square is at most p and check that for all these n there is a remainder concerning the division of p by n.
+
This uses
  +
  +
<haskell>
  +
{-# OPTIONS_GHC -O2 -fno-cse #-}
  +
-- tree-merging Eratosthenes sieve
  +
-- producing infinite list of all prime numbers
  +
primesTME = 2 : gaps 3 (join [[p*p,p*p+2*p..] | p <- primes'])
  +
where
  +
primes' = 3 : gaps 5 (join [[p*p,p*p+2*p..] | p <- primes'])
  +
join ((x:xs):t) = x : union xs (join (pairs t))
  +
pairs ((x:xs):ys:t) = (x : union xs ys) : pairs t
  +
gaps k xs@(x:t) | k==x = gaps (k+2) t
  +
| True = k : gaps (k+2) xs
  +
</haskell>
  +
  +
  +
The tree-merging Eratosthenes sieve here seems to strike a good balance between efficiency and brevity. More at [[Prime numbers]] haskellwiki page. The semi-standard <code>union</code> function is readily available from <hask>Data.List.Ordered</hask>&nbsp;package, put here just for reference:
   
A faster alternative that does not test if ''all'' numbers between 2 and n are ''not'' divisors of p. It rather tests if ''at least one'' number between 2 and n is a divisor of p and when this is the case p can immediately be dismissed and the next number can be tested:
 
   
 
<haskell>
 
<haskell>
isPrime :: Integral a => a -> Bool
+
-- duplicates-removing union of two ordered increasing lists
isPrime p = p > 1 && (not $ any (\n -> p `mod` n == 0) $ takeWhile (\n -> n*n <= p) [2..])
+
union (x:xs) (y:ys) = case (compare x y) of
  +
LT -> x : union xs (y:ys)
  +
EQ -> x : union xs ys
  +
GT -> y : union (x:xs) ys
 
</haskell>
 
</haskell>
   
However, we don't actually need to check all natural numbers <= sqrt P. We need only check the all natural primes <= sqrt P.
+
Here is another solution, intended to be extremely short while still being reasonably fast.
   
 
<haskell>
 
<haskell>
-- Infinite list of all prime numbers
+
isPrime :: (Integral a) => a -> Bool
allPrimes :: [Int]
+
isPrime n | n < 4 = n > 1
allPrimes = filter (isPrime) [2..]
+
isPrime n = all ((/=0).mod n) $ 2:3:[x + i | x <- [6,12..s], i <- [-1,1]]
  +
where s = floor $ sqrt $ fromIntegral n
  +
</haskell>
   
isPrime :: Int -> Bool
+
This one does not go as far as the previous, but it does observe the fact that you only need to check numbers of the form 6k +/- 1 up to the square root. And according to some quick tests (nothing extensive) this version can run a bit faster in some cases, but slower in others; depending on optimization settings and the size of the input.
isPrime p
+
| p < 2 = error "Number too small"
+
''There is a subtle bug in the version above. I'm new here (the wiki and the language) and don't know how corrections are best made (here, or on discussion?). Anyway, the above version will fail on 25, because the bound of s is incorrect. It is x+i that is bounded by the sqrt of the argument, not x. This version will work correctly:''
| p == 2 = True
+
| p > 2 = all (\n -> p `mod` n /= 0) (getPrimes sqrtp)
+
<haskell>
where getPrimes z = takeWhile (<= z) allPrimes
+
isPrime n | n < 4 = n /= 1
sqrtp = floor . sqrt $ fromIntegral p
+
isPrime n = all ((/=0) . mod n) $ takeWhile (<= m) candidates
  +
where candidates = (2:3:[x + i | x <- [6,12..], i <- [-1,1]])
  +
m = floor . sqrt $ fromIntegral n
 
</haskell>
 
</haskell>
   
Note that the mutual dependency of allPrimes and isPrime would result in an infinite loop if we weren't careful. But since we limit our observation of allPrimes to <= sqrt x, we avoid infinite recursion.
 
   
While the mutual dependency is interesting, this second version is not necessarily more efficient than the first. Though we avoid checking all natural numbers <= sqrt P in the isPrime method, we instead check the primality of all natural numbers <= sqrt P in the allPrimes definition.
+
[[Category:Programming exercise spoilers]]

Latest revision as of 19:42, 18 January 2014

(**) Determine whether a given integer number is prime.

Well, a natural number k is a prime number if it is larger than 1 and no natural number n >= 2 with n^2 <= k is a divisor of k. However, we don't actually need to check all natural numbers n <= sqrt k. We need only check the primes p <= sqrt k:

isPrime :: Integral a => a -> Bool
isPrime k = k > 1 &&
   foldr (\p r -> p*p > k || k `rem` p /= 0 && r)
      True primesTME

This uses

{-# OPTIONS_GHC -O2 -fno-cse #-}
-- tree-merging Eratosthenes sieve
--  producing infinite list of all prime numbers
primesTME = 2 : gaps 3 (join [[p*p,p*p+2*p..] | p <- primes'])
  where
    primes' = 3 : gaps 5 (join [[p*p,p*p+2*p..] | p <- primes'])
    join  ((x:xs):t)        = x : union xs (join (pairs t))
    pairs ((x:xs):ys:t)     = (x : union xs ys) : pairs t
    gaps k xs@(x:t) | k==x  = gaps (k+2) t 
                    | True  = k : gaps (k+2) xs


The tree-merging Eratosthenes sieve here seems to strike a good balance between efficiency and brevity. More at Prime numbers haskellwiki page. The semi-standard union function is readily available from
Data.List.Ordered
 package, put here just for reference:


-- duplicates-removing union of two ordered increasing lists
union (x:xs) (y:ys) = case (compare x y) of 
           LT -> x : union  xs  (y:ys)
           EQ -> x : union  xs     ys 
           GT -> y : union (x:xs)  ys

Here is another solution, intended to be extremely short while still being reasonably fast.

isPrime :: (Integral a) => a -> Bool
isPrime n | n < 4 = n > 1
isPrime n = all ((/=0).mod n) $ 2:3:[x + i | x <- [6,12..s], i <- [-1,1]]
            where s = floor $ sqrt $ fromIntegral n

This one does not go as far as the previous, but it does observe the fact that you only need to check numbers of the form 6k +/- 1 up to the square root. And according to some quick tests (nothing extensive) this version can run a bit faster in some cases, but slower in others; depending on optimization settings and the size of the input.

There is a subtle bug in the version above. I'm new here (the wiki and the language) and don't know how corrections are best made (here, or on discussion?). Anyway, the above version will fail on 25, because the bound of s is incorrect. It is x+i that is bounded by the sqrt of the argument, not x. This version will work correctly:

isPrime n | n < 4 = n /= 1 
isPrime n = all ((/=0) . mod n) $ takeWhile (<= m) candidates 
        where candidates = (2:3:[x + i | x <- [6,12..], i <- [-1,1]])
              m = floor . sqrt $ fromIntegral n