Personal tools

99 questions/Solutions/31

From HaskellWiki

< 99 questions | Solutions(Difference between revisions)
Jump to: navigation, search
(putting in a better alternative for allPrimes)
Line 4: Line 4:
 
isPrime :: Integral a => a -> Bool
 
isPrime :: Integral a => a -> Bool
 
isPrime p = p > 1 &&
 
isPrime p = p > 1 &&
(all (\n -> p `mod` n /= 0 )
+
(all ((/= 0).(p `rem`)) $ candidateFactors p)
$ takeWhile (\n -> n*n <= p) [2..])
+
  +
candidateFactors p = takeWhile ((<= p).(^2)) [2..]
 
</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''.
+
Well, a natural number ''p'' is a prime number if it is larger than '''1''' and no natural number ''n >= 2'' with ''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 non-zero remainder concerning the division of ''p'' by ''n''.
   
However, we don't actually need to check all natural numbers <= sqrt P. We need only check the natural ''primes'' <= sqrt P.
+
However, we don't actually need to check all natural numbers ''<= sqrt P''. We need only check the ''primes <= sqrt P'':
   
 
<haskell>
 
<haskell>
 
-- Infinite list of all prime numbers
 
-- Infinite list of all prime numbers
allPrimes :: [Int]
+
{-# OPTIONS_GHC -O2 -fno-cse #-}
allPrimes = filter (isPrime) [2..]
+
candidateFactors p = let z = floor $ sqrt $ fromIntegral p + 1 in
  +
takeWhile (<= z) primesTME
   
isPrime :: Int -> Bool
+
-- tree-merging Eratosthenes sieve
isPrime p
+
primesTME = 2 : gaps 3 (join [[p*p,p*p+2*p..] | p <- primes'])
| p < 2 = error "Number too small"
+
where
| p == 2 = True
+
primes' = 3 : gaps 5 (join [[p*p,p*p+2*p..] | p <- primes'])
| p > 2 = all (\n -> p `mod` n /= 0) (getPrimes sqrtp)
+
join ((x:xs):t) = x : union xs (join (pairs t))
where getPrimes z = takeWhile (<= z) allPrimes
+
pairs ((x:xs):ys:t) = (x : union xs ys) : pairs t
sqrtp = floor . sqrt $ fromIntegral p
+
gaps k xs@(x:t) | k==x = gaps (k+2) t
</haskell>
+
| True = k : gaps (k+2) xs
   
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.
+
-- 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
  +
</haskell>
   
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.
+
The tree-merging Eratosthenes sieve here seems to strike a good balance between efficiency and brevity. More at [[Prime numbers]] haskellwiki page.

Revision as of 07:58, 31 May 2011

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

isPrime :: Integral a => a -> Bool
isPrime p = p > 1 && 
            (all ((/= 0).(p `rem`)) $ candidateFactors p)
 
candidateFactors p = takeWhile ((<= p).(^2)) [2..]

Well, a natural number p is a prime number if it is larger than 1 and no natural number n >= 2 with 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 non-zero remainder concerning the division of p by n.

However, we don't actually need to check all natural numbers <= sqrt P. We need only check the primes <= sqrt P:

-- Infinite list of all prime numbers
{-# OPTIONS_GHC -O2 -fno-cse #-}
candidateFactors p = let z = floor $ sqrt $ fromIntegral p + 1 in
                     takeWhile (<= z) primesTME
 
-- tree-merging Eratosthenes sieve
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
 
-- 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

The tree-merging Eratosthenes sieve here seems to strike a good balance between efficiency and brevity. More at Prime numbers haskellwiki page.