Lazy List of Prime Numbers

Question

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.

Yoshida

Something like #1764663? – Duncandunce

Consider hackage.haskell.org/package/primes – Comfy

Quite the contrary: it's a tricky task to create non-lazy prime numbers list in Haskell – Billetdoux

by walpen at codegolf: 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. – Titration

@WillNess; that yields all numbers here. Maybe you meant nubBy (((>1).).gcd) [2..]? – Ingathering

@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). – Titration

I’m on GHC-7.10 right now. See https://mcmap.net/q/534438/-nubby-is-not-working-as-expected for a rationale of the change. – Ingathering

@JoachimBreitner so it was working on the older versions (though shouldn't have). OK. – Titration

FWIW, it is possible to find one nth prime directly, independently of all the previous ones. (see "The fast way" section in that answer). – Titration

Answer 1

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.

Evaevacuant

Please read this and rethink your answer: cs.hmc.edu/~oneill/papers/Sieve-JFP.pdf – Hapsburg

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, it instead runs at about O(n^1.4..1.45) (in n primes produced), just like any other normal trial division. – Titration

Answer 2

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

Stumper

Answer 3

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.

Elysia

Answer 4

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.

Interrogation

You might want to use [3,5..] instead of [3..] – Grog

you might also want to use (\n -> n*n <= x) or ((<= n) . (^2))instead of (<= (floor . sqrt $ fromIntegral x) – Sparke

Answer 5

I like to keep things separate so that one might come in and ask: "is 1 a (unary) prime number?" ('unary' in the sense that it only needs itself to pass the primality test).

primes :: [Integer]
primes = 2 : 3 : potentialPrimes

Should the community's answer be 'yes', then the above code can start with 1, because the other prime numbers are going to be generated by the potentialPrimes call.

potentialPrimes :: [Integer]
potentialPrimes = 5 : 7 : [p | p <- [6*k+l | k <- [1..], l <- [5,7]], isPrime p]

What it does is to consider 5 and 7 as given, and then generate multiples of the perfect number 6 - its divisors, excluding 6 itself, sum up to it - with 5 or 7 added for potential twin primes. They're going to be checked by isPrime at every p.

isPrime :: Integer -> Bool
isPrime n = all (\p -> n `mod` p /= 0) (takeWhile (\p -> p*p <= n) potentialPrimes)

And last, but not the least, what this does is to consider if every number that has been fetched from the potential primes list, whose square is less than or equal to the one we're testing, satisfies the condition that they do not divide the tested number. It's true when none of them divide it (hence all), so p can be added to the primes' list, when a single exception makes it false and excludes the number.

Flybynight