Personal tools

99 questions/Solutions/39

From HaskellWiki

< 99 questions | Solutions(Difference between revisions)
Jump to: navigation, search
Line 6: Line 6:
 
<haskell>
 
<haskell>
 
primesR :: Integral a => a -> a -> [a]
 
primesR :: Integral a => a -> a -> [a]
primesR a b = filter isPrime [a..b]
+
primesR a b | even a = filter isPrime [a+1,a+3..b]
  +
| True = filter isPrime [a,a+2..b]
 
</haskell>
 
</haskell>
   
Line 26: Line 26:
   
 
'''Solution 3.'''
 
'''Solution 3.'''
  +
  +
Another way to compute the claimed list is done by using the ''Sieve of Eratosthenes''.
 
<haskell>
 
<haskell>
 
primesR :: Integral a => a -> a -> [a]
 
primesR :: Integral a => a -> a -> [a]
Line 32: Line 34:
 
</haskell>
 
</haskell>
   
Another way to compute the claimed list is done by using the ''Sieve of Eratosthenes''. The <code>sieve [2..]</code> function call generates a list of all (!) prime numbers using this algorithm and <code>primesR</code> filters 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 :)]
+
The <code>sieve [2..]</code> function call generates a list of all (!) prime numbers using this algorithm and <code>primesR</code> filters 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 :)]
   
 
''this is of course a famous case of (mislabeled) executable specification, with all the implied pitfalls of inefficiency when (ab)used as if it were an actual code''.
 
''this is of course a famous case of (mislabeled) executable specification, with all the implied pitfalls of inefficiency when (ab)used as if it were an actual code''.
Line 46: Line 48:
 
| b<a || b<2 = []
 
| b<a || b<2 = []
 
| otherwise =
 
| otherwise =
(if a <= 2 then [2] else []) ++
+
(if a <= 2 then [2] else []) ++
gaps a' (join [[x,x+step..b] | p <- takeWhile (<= z) primes'
+
(takeWhile (<= b) $ gaps a' $
, let q = p*p ; step = 2*p
+
join [[x,x+step..] | p <- takeWhile (<= z) primes' -- p^2 <= b
x = snapUp (max a' q) q step ])
+
, let q = p*p ; step = 2*p
  +
x = snap (max a' q) q step ])
 
where
 
where
 
primes' = tail primesTME -- external unbounded list of primes
 
primes' = tail primesTME -- external unbounded list of primes
a' = snapUp (max 3 a) 3 2
+
a' = snap (max 3 a) 3 2
 
z = floor $ sqrt $ fromIntegral b + 1
 
z = floor $ sqrt $ fromIntegral b + 1
 
join (xs:t) = union xs (join (pairs t))
 
join (xs:t) = union xs (join (pairs t))
Line 60: Line 62:
 
gaps k xs@(x:t) | k==x = gaps (k+2) t
 
gaps k xs@(x:t) | k==x = gaps (k+2) t
 
| True = k : gaps (k+2) xs
 
| True = k : gaps (k+2) xs
gaps k [] = [k,k+2..b]
+
snap v origin step = let r = rem (v-origin) step -- v >= origin
snapUp v origin step = let r = rem (v-origin) step -- rem OK if v>=origin
+
in if r==0 then v else v+(step-r)
in if r==0 then v else v+(step-r)
 
 
-- duplicates-removing union of two ordered increasing lists
 
-- duplicates-removing union of two ordered increasing lists
 
union (x:xs) (y:ys) = case (compare x y) of
 
union (x:xs) (y:ys) = case (compare x y) of
Line 67: Line 69:
 
EQ -> x : union xs ys
 
EQ -> x : union xs ys
 
GT -> y : union (x:xs) ys
 
GT -> y : union (x:xs) ys
union a b = a ++ b
 
 
</haskell>
 
</haskell>
 
''(This turned out to be quite a project, with some 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,
 
''(This turned out to be quite a project, with some 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,
Line 73: Line 74:
 
> primesR 10100 10200 -- Sol.4
 
> primesR 10100 10200 -- Sol.4
 
[10103,10111,10133,10139,10141,10151,10159,10163,10169,10177,10181,10193]
 
[10103,10111,10133,10139,10141,10151,10159,10163,10169,10177,10181,10193]
(5,497 reductions, 11,382 cells)
+
(4,656 reductions, 10,859 cells)
   
 
> takeWhile (<= 10200) $ dropWhile (< 10100) $ primesTME -- Sol.2
 
> takeWhile (<= 10200) $ dropWhile (< 10100) $ primesTME -- Sol.2
Line 84: Line 85:
 
(54,893,566 reductions, 79,935,263 cells, 6 garbage collections)
 
(54,893,566 reductions, 79,935,263 cells, 6 garbage collections)
   
> filter isPrime [10100..10200] -- Sol.1
+
> filter isPrime [10101,10103..10200] -- Sol.1
 
[10103,10111,10133,10139,10141,10151,10159,10163,10169,10177,10181,10193]
 
[10103,10111,10133,10139,10141,10151,10159,10163,10169,10177,10181,10193]
(15,750 reductions, 29,292 cells) -- isPrime: Q.31
+
(12,927 reductions, 24,703 cells) -- isPrime: Q.31
 
</haskell>
 
</haskell>
   

Revision as of 09:50, 5 June 2011

(*) 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 | even a = filter isPrime [a+1,a+3..b]
            | True   = filter isPrime [a,a+2..b]

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

This is good for very narrow ranges as Q.31's isPrime tests numbers by trial division using (up to\sqrt b) a memoized primes list produced by sieve of Eratosthenes to which it refers internally. So it'll be slower, but immediate, testing the numbers one by one.

Solution 2.

For very wide ranges, specifically when a < \sqrt{b}, we're better off just using the primes sequence itself, without any post-processing:

primes :: Integral a => [a]
primes = primesTME               -- of Q.31
 
primesR :: Integral a => a -> a -> [a]
primesR a b = takeWhile (<= b) $ dropWhile (< a) primes

Solution 3.

Another way to compute the claimed list is done by using the Sieve of Eratosthenes.

primesR :: Integral a => a -> a -> [a]
primesR a b = takeWhile (<= b) $ dropWhile (< a) $ sieve [2..]
  where sieve (n:ns) = n:sieve [ m | m <- ns, m `mod` n /= 0 ]

The sieve [2..] function call generates a list of all (!) prime numbers using this algorithm and primesR filters 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 :)]

this is of course a famous case of (mislabeled) executable specification, with all the implied pitfalls of inefficiency when (ab)used as if it were an actual code.

Solution 4.

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 Q.31, 
--  adjusted to produce primes in a given range (inclusive)
primesR a b 
  | b<a || b<2 = []
  | otherwise  = 
     (if a <= 2 then [2] else []) ++ 
     (takeWhile (<= b) $ gaps a' $
        join [[x,x+step..] | p <- takeWhile (<= z) primes'  -- p^2 <= b
                       , let q = p*p ; step = 2*p
                             x = snap (max a' q) q step ])
  where
    primes' = tail primesTME        -- external unbounded list of primes
    a'      = snap (max 3 a) 3 2
    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
    snap v origin step = let r = rem (v-origin) step        -- v >= origin
                         in if r==0 then v else v+(step-r)
    -- 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

(This turned out to be quite a project, with some 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.4
[10103,10111,10133,10139,10141,10151,10159,10163,10169,10177,10181,10193]
(4,656 reductions, 10,859 cells)
 
> takeWhile (<= 10200) $ dropWhile (< 10100) $ primesTME         -- Sol.2
[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.3
     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 [10101,10103..10200]                            -- Sol.1
[10103,10111,10133,10139,10141,10151,10159,10163,10169,10177,10181,10193]
(12,927 reductions, 24,703 cells)                        -- isPrime: Q.31

(testing with Hugs of Nov 2002).

This solution is faster but not immediate. It has a certain preprocessing stage but then goes on fast to produce the whole range. To illustrate, to produce the 49 primes in 1000-wide range above 120200300100 it takes about 18 seconds on my oldish notebook for the 1st version, with the first number produced almost immediately (~ 0.4 sec); but this version spews out all 49 primes in one go after just under 1 sec.