27

How would one implement a list of prime numbers in Haskell so that they could be retrieved lazily?

I am new to Haskell, and would like to learn about practical uses of the lazy evaluation functionality.

Josh Lee
  • 149,877
  • 34
  • 253
  • 263
Mantas Vidutis
  • 15,230
  • 20
  • 73
  • 91
  • 2
    Something like http://stackoverflow.com/questions/1764163/help-explain-this-chunk-of-haskell-code-that-outputs-a-stream-of-primes? – kennytm Aug 29 '10 at 20:44
  • 1
    Consider http://hackage.haskell.org/package/primes – Don Stewart Aug 29 '10 at 23:17
  • Quite the contrary: it's a tricky task to create non-lazy prime numbers list in Haskell – vpolozov Aug 31 '10 at 18:31
  • [by walpen at codegolf](http://codegolf.stackexchange.com/a/6316/5021): `nubBy (((==0).).rem) [2..]`. To try it out in GHCi first bring up the `Data.List` module with `Prelude> :m +Data.List`. But lazyness plays no role here, except allowing for the *unbounded* definition. `[2..10000]` could be used as well and evaluated strictly. – Will Ness Sep 22 '12 at 08:44
  • @WillNess; that yields all numbers here. Maybe you meant `nubBy (((>1).).gcd) [2..]`? – Joachim Breitner Apr 19 '16 at 08:12
  • @JoachimBreitner yes; the former code was working on older versions of GHC (it works e.g. on 7.8.3). The newer ones flipped the order of arguments to the nubBy's function, IIRC. or is it the other way around and you're on an *older* version? (because the code with `rem` also works fine on tryhaskell.org). – Will Ness Apr 19 '16 at 13:47
  • 1
    I’m on GHC-7.10 right now. See http://stackoverflow.com/a/33533257/946226 for a rationale of the change. – Joachim Breitner Apr 20 '16 at 07:14
  • @JoachimBreitner so it was working on the *older* versions (though shouldn't have). OK. – Will Ness Apr 21 '16 at 19:25
  • FWIW, it *is* possible to [find one *n*th prime directly](https://stackoverflow.com/a/9704912/849891), independently of all the previous ones. (see **"The fast way"** section in that answer). – Will Ness Jan 06 '18 at 11:03

4 Answers4

24

Here's a short Haskell function that enumerates primes from Literate Programs: 

primes :: [Integer]
primes = sieve [2..]
  where
    sieve (p:xs) = p : sieve [x|x <- xs, x `mod` p > 0]

Apparently, this is not the Sieve of Eratosthenes (thanks, Landei). I think it's still an instructive example that shows you can write very elegant, short code in Haskell and that shows how the choice of the wrong data structure can badly hurt efficiency.

Will Ness
  • 62,652
  • 8
  • 86
  • 167
Niki
  • 15,188
  • 5
  • 41
  • 72
  • 19
    Please read this and rethink your answer: http://www.cs.hmc.edu/~oneill/papers/Sieve-JFP.pdf – Landei Aug 30 '10 at 08:19
  • 3
    The "wrong data structure" (i.e. list) has nothing to do with that code's *extreme* inefficiency ( O(n^2), *in `n` primes produced* ), which is instead the result of premature firing up of filters on each newly found prime instead of on its *square*. With filters creation [correctly postponed](http://www.haskell.org/haskellwiki/Prime_numbers#Postponed_Filters), it instead runs at about O(n^1.4..1.45) (in `n` primes produced), just like any other normal trial division. – Will Ness Feb 12 '12 at 00:44
10

There are a number of solutions for lazy generation of prime sequences right in the haskell wiki. The first and simplest is the Postponed Turner sieve: (old revision ... NB)

primes :: [Integer]
primes = 2: 3: sieve (tail primes) [5,7..]
 where 
  sieve (p:ps) xs = h ++ sieve ps [x | x <- t, x `rem` p /= 0]  
                                -- or:  filter ((/=0).(`rem`p)) t
                  where (h,~(_:t)) = span (< p*p) xs
Will Ness
  • 62,652
  • 8
  • 86
  • 167
sleepynate
  • 7,296
  • 2
  • 24
  • 38
4

The accepted answer from @nikie is not very efficient, is gets relatively slow after some thousands, but the answer of @sleepynate is much better. It took me some time to understand it, therefore here is the same code, but just with variables named more clearly:

lazyPrimes :: [Integer]
lazyPrimes = 2: 3: calcNextPrimes (tail lazyPrimes) [5, 7 .. ]
  where
    calcNextPrimes (p:ps) candidates =
      let (smallerSquareP, (_:biggerSquareP)) = span (< p * p) candidates in
      smallerSquareP ++ calcNextPrimes ps [c | c <- biggerSquareP, rem c p /= 0]

The main idea is that the candidates for the next primes already contain no numbers that are divisible by any prime less than the first prime given to the function. So that if you call

calcNextPrimes (5:ps) [11,13,17..]

the candidate list contains no number, that is divisible by 2 or 3, that means that the first non-prime candidate will be 5 * 5, cause 5* 2 and 5 * 3 and 5 * 4 are already eliminated. That allows you to take all candidates, that are smaller than the square of 5 and add them straight away to the primes and sieve the rest to eliminate all numbers divisible by 5.

Daria
  • 61
  • 5
3
primes = 2 : [x | x <- [3..], all (\y -> x `mod` y /= 0) 
                   (takeWhile (<= (floor . sqrt $ fromIntegral x)) primes)]

With 2 in the list initially, for each integer x greater than 2, check if for all y in primes such that y <= sqrt(x), x mod y != 0 holds, which means x has no other factors except 1 and itself.

Will Ness
  • 62,652
  • 8
  • 86
  • 167
Tianren Li
  • 31
  • 1