Finding primes using a primes map with Haskell and Hugs98

[Joe English]

Shlomi Fish wrote:

> As some of you may know, a Haskell program that prints all the primes can be
> as short as the following:
>
> primes = sieve [2.. ] where
>          sieve (p:x) = p : sieve [ n | n <- x, n `mod` p > 0 ]
>
> Now, this program roughly corresponds to the following perl program:

[ ~20 line Perl program snipped ]

> The program can be more optimized for both speed and code size, but I wanted
> to make it as verbose as possible.
>
> There is a different algorithm which keeps a boolean map [...]
> The algorithm iterates over all the numbers from 2 to the square root
> of the desired bound, and if it encounters a prime number it marks all the
> numbers p*p, p*p+p, p*p+2*p, p*p+3*p, etc. as not prime.

[~40 line Perl implementation snipped]

> Now, I tried writing an equivalent Haskell program and the best I
> could do was the following:

[ ~45 line Haskell implementation snipped ]

Another way to do this is to compute the final array directly,
instead of computing successive versions of the array:

    import Array
    primes n = [ i | i <- [2 ..n], not (primesMap ! i)] where
	primesMap   = accumArray (||) False (2,n) multList
	multList    = [(m,True) | j <- [2 .. n `div` 2], m <- multiples j]
	multiples j = takeWhile (n>=) [k*j | k <- [2..]]

Now this version does a lot more work than the algorithm
described above -- it computes multiples of *all* the integers
less than n/2, not just the primes less than sqrt(n) -- but
it has the virtue of being short enough to reason about effectively
and is probably a better starting point for further optimization.

> The problem is that when running it on hugs98 on a Windows98 computer with
> 64MB of RAM, I cannot seem to scale beyond 30,000 or so, as my boundary. When
> entering how_much as 50,000 I get the following message:
>
> ERROR: Garbage collection fails to reclaim sufficient space

My implementation fares even worse under Hugs -- it runs out
of space around n = 4500 (Linux box, 64M RAM).  With GHC
it has no problem for n = 100,000, although the space usage
is still extremely poor.  It grows to consume all
available RAM at around n = 200,000.  (On the other hand,
it's considerably faster than the traditional 2-liner
listed above, up to the point where it starts paging).

I suspect the poor memory usage is due to the way accumArray
works -- it's building up a huge array of suspensions of the form

	(False && (False && ( ... && True)))

that aren't reduced until an array element is requested.

(A strict version of accumArray, analogous to "foldl_strict"
defined below, would solve this problem, but I don't
see any way to implement it in Standard Haskell).

> In perl I can scale beyond 100,000, and if I modify the code to use a bit
> vector (using vec) to much more. So my question is what am I or hugs are
> doing wrong and how I can write better code that implements this specific
> algorithm.
>
> From what I saw I used tail recursion, (and hugs98 has proper tail recursion
> right?), and there's only one primes_map present at each iteration (and thus,
> at all), so it shouldn't be too problematic.

Actually no; this is a common misconception.  In a strict
language like Scheme, tail call optimization works because
a tail call is the last thing a function does.  In Haskell
though the tail call is the *first* thing that gets evaluated
(more or less), leaving all the "earlier" work as an unevaluated
suspension.  Code that is space-efficient in a strict language
frequently suffers from awful space leaks in a lazy language.
For example:

sum_first_n_integers n = f n 0 where
	f 0 a = a
	f n a = f (n-1) (n+a)

quickly leads to a "Control Stack Overflow" error in Hugs.
BTW, the trick to fix it is to change the last line to:

	f n acc = f (n-1) $! (n+acc)

or to replace the whole thing with:

    foldl_strict (+) 0 [1..n]

where

    foldl_strict f a []     = a
    foldl_strict f a (x:xs) = (foldl_strict f $! f a x) xs


> Does it have to do with the way hugs98 implements and Int to Bool array?

Most likely yes.  Hugs is optimized for interactive use and quick
compilation, not for space usage.  Try it with GHC or HBC and
see how it does.


--Joe English

  jenglish@flightlab.com