Difference between revisions of "99 questions/Solutions/31"

From HaskellWiki
Jump to: navigation, search
(adding an interesting variant)
(Converted a discussion on-page to a paragraph. Please add discussions to the discussions page only.)
 
(19 intermediate revisions by 7 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:
   
However, we don't actually need to check all natural numbers <= sqrt P. We need only check the all natural primes <= sqrt P.
 
   
 
<haskell>
 
<haskell>
-- Infinite list of all prime numbers
 
  +
-- duplicates-removing union of two ordered increasing lists
allPrimes :: [Int]
 
  +
union (x:xs) (y:ys) = case (compare x y) of
allPrimes = filter (isPrime) [2..]
 
  +
LT -> x : union xs (y:ys)
  +
EQ -> x : union xs ys
  +
GT -> y : union (x:xs) ys
  +
</haskell>
   
isPrime :: Int -> Bool
 
  +
Here is another solution, intended to be extremely short while still being reasonably fast.
isPrime p
 
  +
| p < 2 = error "Number too small"
 
  +
<haskell>
| p ≡ 2 = True
 
  +
isPrime :: (Integral a) => a -> Bool
| p > 2 = all (\n -> p `mod` n /= 0) (getPrimes sqrtp)
 
  +
isPrime n | n < 4 = n > 1
where getPrimes z = takeWhile (≤ z) allPrimes
 
  +
isPrime n = all ((/=0).mod n) $ 2:3:[x + i | x <- [6,12..s], i <- [-1,1]]
sqrtp = floor∘sqrt $ fromIntegral p
+
where s = floor $ sqrt $ fromIntegral n
  +
</haskell>
  +
  +
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. 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:
  +
  +
<haskell>
  +
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
 
</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 07:06, 11 May 2016

(**) 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. 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