99 questions/Solutions/39

From HaskellWiki
Jump to navigation Jump to search

(*) A list of prime numbers.

Given a range of integers by its lower and upper limit, construct a list of all prime numbers in that range.

Solution 1:

primesR :: Integral a => a -> a -> [a]
primesR a b = filter isPrime [a..b]

If we are challenged to give all primes in the range between a and b we simply take all number from a up to b and filter the primes out.

Solution 2:

primes :: Integral a => [a]
primes = let sieve (n:ns) = n:sieve [ m | m <- ns, m `mod` n /= 0 ] 
         in sieve [2..]

primesR :: Integral a => a -> a -> [a]
primesR a b = takeWhile (<= b) $ dropWhile (< a) primes

Another way to compute the claimed list is done by using the Sieve of Eratosthenes. The primes function generates a list of all (!) prime numbers using this algorithm and primesR filter the relevant range out. [But this way is very slow and I only presented it because I wanted to show how nicely the Sieve of Eratosthenes can be implemented in Haskell :)]

Solution 3:

Use the proper Sieve of Eratosthenes from e.g. 31st question's solution (instead of the above sieve of Turner), adjusted to start its multiples production from the given start point:

{-# OPTIONS_GHC -O2 -fno-cse #-}
-- tree-merging Eratosthenes sieve, primesTME of haskellwiki/prime_numbers, 
--  adjusted to produce primes in a given range
primesR a b 
  | b<a || b<2 = []
  | otherwise  = 
     (if a <= 2 then [2] else []) ++
     gaps a' (join [[x,x+step..b] | p <- takeWhile (<= z) (tail primesTME)
                    , let q = p*p ; step = 2*p
                          x = if a' <= q then q else snapUp a' q step ])
  where
    a'      = if a<=3 then 3 else (if even a then a+1 else a)
    z       = floor $ sqrt $ fromIntegral b + 1
    join  (xs:t)    = union xs (join (pairs t))
    join  []        = []
    pairs (xs:ys:t) = (union xs ys) : pairs t
    pairs  t        = t
    gaps k xs@(x:t) | k==x  = gaps (k+2) t 
                    | True  = k : gaps (k+2) xs
    gaps k []       = [k,k+2..b]

    snapUp v origin step = let r = rem (v-origin) step
                           in if r==0 then v else v-r+step
    -- 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
    union  a      b     = a ++ b

(This turned out to be quite a project, with a few quite subtle points.) It should be much better then taking a slice of a full sequential list of primes, as it won't try to generate any primes between the square root of b and a. To wit,

> primesR 10100 10200                                            -- Sol.3
[10103,10111,10133,10139,10141,10151,10159,10163,10169,10177,10181,10193]
(6,038 reductions, 11,923 cells)

> takeWhile (<= 10200) $ dropWhile (< 10100) $ primesTME  -- TME: of Q.31
[10103,10111,10133,10139,10141,10151,10159,10163,10169,10177,10181,10193]
(140,313 reductions, 381,058 cells)

> takeWhile (<= 10200) $ dropWhile (< 10100) $ sieve [2..]       -- Sol.2
     where sieve (n:ns) = n:sieve [ m | m <- ns, m `mod` n /= 0 ]
[10103,10111,10133,10139,10141,10151,10159,10163,10169,10177,10181,10193]
(54,893,566 reductions, 79,935,263 cells, 6 garbage collections)

> filter isPrime [10100..10200]                                  -- Sol.1
[10103,10111,10133,10139,10141,10151,10159,10163,10169,10177,10181,10193]
(34,694 reductions, 55,146 cells)  -- isPrime: Q.31-Sol.1 using primesTME
(15,750 reductions, 29,292 cells)  -- isPrime: Q.31-Sol.2 using primesTME

(testing with Hugs of Nov 2002).