Difference between revisions of "99 questions/Solutions/31"
< 99 questions | Solutions
Jump to navigation
Jump to 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 ( |
+ | (all ((/= 0).(p `rem`)) $ candidateFactors p) |
| + | |||
| ⚫ | |||
| + | candidateFactors p = takeWhile ((<= p).(^2)) [2..] |
||
</haskell> |
</haskell> |
||
| − | Well, a natural number ''p'' is a prime number |
+ | 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 |
+ | 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 |
||
| + | {-# OPTIONS_GHC -O2 -fno-cse #-} |
||
| − | allPrimes :: [Int] |
||
| + | candidateFactors p = let z = floor $ sqrt $ fromIntegral p + 1 in |
||
| − | allPrimes = filter (isPrime) [2..] |
||
| ⚫ | |||
| + | -- tree-merging Eratosthenes sieve |
||
| − | isPrime :: Int -> Bool |
||
| + | primesTME = 2 : gaps 3 (join [[p*p,p*p+2*p..] | p <- primes']) |
||
| − | isPrime p |
||
| + | where |
||
| − | | p < 2 = error "Number too small" |
||
| + | primes' = 3 : gaps 5 (join [[p*p,p*p+2*p..] | p <- primes']) |
||
| − | | p == 2 = True |
||
| − | + | join ((x:xs):t) = x : union xs (join (pairs t)) |
|
| + | pairs ((x:xs):ys:t) = (x : union xs ys) : pairs t |
||
| − | where getPrimes z = takeWhile (<= z) allPrimes |
||
| − | + | 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 |
||
| − | 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. |
||
| + | 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. |
||
| − | 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. |
||
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.