Difference between revisions of "99 questions/Solutions/31"
(adding an interesting variant) 
(Converted a discussion onpage 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 fnocse #} 

+   treemerging 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 treemerging Eratosthenes sieve here seems to strike a good balance between efficiency and brevity. More at [[Prime numbers]] haskellwiki page. The semistandard <code>union</code> function is readily available from <hask>Data.List.Ordered</hask> 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 

+   duplicatesremoving 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]] 

−  +  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 fnocse #}
 treemerging 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 treemerging Eratosthenes sieve here seems to strike a good balance between efficiency and brevity. More at Prime numbers haskellwiki page. The semistandard union
function is readily available from Data.List.Ordered
package, put here just for reference:
 duplicatesremoving 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