CODEWITHSUNDEEP
ocamlprimes
Xophmeister
Xophmeister

Reputation: 9211

Generating primes without blowing the stack

I'm learning OCaml (so forgive my style) and am trying to write a function that generates a list of prime numbers up to some upper bound. I've managed to do this in several different ways, all of which work until you scale them to a relatively high upper bound.

How can I change these (any of them) so that the recursion doesn't fill up the stack? I thought my while loop version would achieve this, but apparently not!

Generator

let primes max =
  let isPrime p x =
    let hasDivisor a = (x mod a = 0) in
    not (List.exists hasDivisor p) in

  let rec generate p test =
    if test < max then
      let nextTest = test + 2 in
      if isPrime p test then generate (test :: p) nextTest
                        else generate p nextTest
    else p in

  generate [5; 3; 2] 7;;

This has been my most successful solution insofar as it doesn't immediately overflow the stack when running primes 2000000;;. However it just sits there consuming CPU; I can only assume it will complete eventually! The following alternatives all have the stack overflow problem:

Recursive Sieve of Eratosthenes

let primes max =
  let rec sieve found toTest =
    let h = List.hd toTest
    and t = List.tl toTest in

    let newPrimes = h :: found
    and doesntDivide x = (x mod h <> 0) in

    let nonDivisors = List.filter doesntDivide t in
      if nonDivisors = [] then newPrimes
                          else sieve newPrimes nonDivisors in

  let rec range a b =
    if a > b then []
             else a :: range (a + 1) b in

  let p = range 2 max in

  sieve [] p;;

Recursive Sieve of Eratosthenes v2

let primes max =
  let rec sieve toTest =
    let h = List.hd toTest
    and t = List.tl toTest in
    let doesntDivide x = (x mod h <> 0) in
    let nonDivisors = List.filter doesntDivide t in
      if nonDivisors = [] then [h]
                          else (h :: sieve nonDivisors) in

  let rec range a b =
    if a > b then []
             else a :: range (a + 1) b in

  let p = range 2 max in

  sieve p;;

While Loop Sieve of Eratosthenes

let primes max =
  let rec range a b =
    if a > b then []
             else a :: range (a + 1) b in

  let tail = ref (range 2 max)
  and p    = ref [] in

  while !tail <> [] do
    let h = List.hd !tail
    and t = List.tl !tail in
    let doesntDivide x = (x mod h <> 0) in
    let newTail = ref (List.filter doesntDivide t) in

    tail := !newTail;
    p := h :: !p
  done;

  !p;;

Upvotes: 2

Views: 2030

Answers (2)

seanmcl
seanmcl

Reputation: 9946

The stack overflows occur because your range function is not tail recursive. One that works is, e.g. 

  let rec range store a b =
    if a > b then store
    else range (a :: store) (a + 1) b
  in

  let p = List.rev (range [] 2 max) in

With that definition, and a format line, gives

$ ocamlopt -o primes2 primes2.ml
$ ./primes2
2
3
5
7
11
13
17
19
23
29
31
37
41
43
47
53
59
...

Since you're learning, I'll give you some unsolicited style comments as well :)

  • Don't use hd and tl. Prefer pattern matching. Then the compiler can tell you the cases you missed. E.g.

    let rec sieve found toTest = let h = List.hd toTest and t = List.tl toTest in

would be

let rec sieve found = function
  | h :: t -> ...
  | [] -> Error handling...

  • Don't use x = []. Use pattern patching.

    match x with | [] -> ... | h::t -> ...

  • Use anonymous functions rather than short (i.e. <= 1 line) named single use functions:

    let doesntDivide x = (x mod h <> 0) in let nonDivisors = List.filter doesntDivide t in

    let nonDivisors = List.filter (fun x -> (x mod h <> 0)) t in

  • Use imperative features sparingly.

Upvotes: 9

user448810
user448810

Reputation: 17866

Your algorithms that you claim are the Sieve of Eratosthenes actually are not; they use trial division instead of sieving, which is easy to spot by looking for a comparison of a remainder (the mod operator) to zero. Here's a simple implementation of the Sieve of Eratosthenes, in pseudocode instead of Ocaml because it's been a long time since I wrote Ocaml code:

function primes(n)
    sieve := makeArray(2..n, True)
    for p from 2 to n
        if sieve[p]
            output p
            for i from p*p to n step p
                sieve[i] := False

That can be optimized further, though for small limits like n = 2000000 there is little point in doing so; in any case, a sieve will be very much faster than the trial division that you are using. If you're interested in programming with prime numbers, I modestly recommend this essay at my blog.

Upvotes: 4

Related Questions