Difference between revisions of "99 questions/Solutions/39"
Line 40: | Line 40: | ||
gaps a' (join [[x,x+step..b] | p <- takeWhile (<= z) primes' |
gaps a' (join [[x,x+step..b] | p <- takeWhile (<= z) primes' |
||
, let q = p*p ; step = 2*p |
, let q = p*p ; step = 2*p |
||
− | x = |
+ | x = snapUp (max a' q) q step ]) |
where |
where |
||
primes' = tail primesTME -- external unbounded list of primes |
primes' = tail primesTME -- external unbounded list of primes |
||
− | a' = |
+ | a' = snapUp (max 3 a) 1 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 52: | Line 52: | ||
| True = k : gaps (k+2) xs |
| True = k : gaps (k+2) xs |
||
gaps k [] = [k,k+2..b] |
gaps k [] = [k,k+2..b] |
||
− | snapUp v origin step = let r = rem (v-origin) step |
+ | snapUp v origin step = let r = rem (v-origin) step -- rem OK if v>=origin |
− | in if r==0 then v else v |
+ | 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 65: | Line 65: | ||
> primesR 10100 10200 -- Sol.3 |
> primesR 10100 10200 -- Sol.3 |
||
[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, |
+ | (5,497 reductions, 11,382 cells) |
> takeWhile (<= 10200) $ dropWhile (< 10100) $ primesTME -- TME of Q.31 |
> takeWhile (<= 10200) $ dropWhile (< 10100) $ primesTME -- TME of Q.31 |
Revision as of 10:51, 3 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 = filter isPrime [a..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 by trial division using (up to) a memoized primes list produced by sieve of Eratosthenes to which it refers internally. So it'll be slower, but immediate.
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 :)]
this is of course a famous case of executable specification, with all the implied pitfalls of inefficiency when (ab)used as if it were an actual code.
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 (inclusive)
primesR a b
| b<a || b<2 = []
| otherwise =
(if a <= 2 then [2] else []) ++
gaps a' (join [[x,x+step..b] | p <- takeWhile (<= z) primes'
, let q = p*p ; step = 2*p
x = snapUp (max a' q) q step ])
where
primes' = tail primesTME -- external unbounded list of primes
a' = snapUp (max 3 a) 1 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
gaps k [] = [k,k+2..b]
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)
-- 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 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.3
[10103,10111,10133,10139,10141,10151,10159,10163,10169,10177,10181,10193]
(5,497 reductions, 11,382 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]
(15,750 reductions, 29,292 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 up all 49 primes in one go after just under 1 sec.
Solution 4.
For very wide ranges, specifically when , 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