Your algorithm seems to be based on the following idea: calculate the non-primes and derive the primes from them by calculating the set difference of the natural numbers and the non-primes. A naive implementation of this idea can be found as primes' in the attachached file. The function uses no multiplication or division and though performs 6 times worse than the sieve in calculating the first 30000 primes. The complexity for finding the next i'th prime with this naive implementation is about O(i). In comparison to this, the sieve provides a good optimization because only those natural numbers are tested against the i'th prime which have run through all other sieves. Nevertheless, your algorithm is promising when the non-primes are merged efficiently enough into a single sorted list which can be easily subtracted from the natural numbers. I think the deployment of an array is basically a way to efficiently merge the multiples of the primaries into a sorted list (where even duplicates are removed), thus hoping to reduce the number of the operations better than the optimization that is provided by the sieve. However, to use arrays this way, you probably need destructive array updates, because the array must be incrementally updated when new primes are found. I think that standard haskell arrays don't do the job very well. An implementation of the "merging" idea in Haskell is primes'' in the attached file. It is 15% faster then the sieve in calculating the 30000 first primes. The algorithm is realized as two mutually recursive functions noprimes and primes'', the latter calculating the set difference between the non-primes and the natural numbers, the former merging the all multiples of all primes into a sorted list. It should be possible to substantially optimize the merging operation. primes''' is an efficient variant of primes'. Instead of a list it uses a binary tree for the management of the lists of multiples of the already found primes, and thus requires some additional programming effort. The complexity is reduced from O(i) to something like O(Log(i)). Compared with the sieve, primes''' needs only half the time to calculate the first 30000 primes. (Tests with ghc 4.08, 64m heap) Best, Elke. On 15-Dec-00 Shlomi Fish wrote:
Hi!
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:
###### SNIP SNIP ##### #!/usr/bin/perl
use strict;
my (@primes, $a, $p); @primes = (2); MAIN_LOOP: for($a = 3; $a < 1000; $a++) { foreach $p (@primes) { if ($a % $p == 0) { next MAIN_LOOP; } } push @primes, $a; } print join(", ", @primes); ####### SNIP SNIP #####
The program can be more optimized for both speed and code size, but I wanted to make it as verbose as possible.
The algorithm keeps a list of the primes, and for each new number checks if it is divisable by any of them and if not it adds it to the list.
There is a different algorithm which keeps a boolean map which tells whether the number at that position is prime or not. At start it is initialized to all trues. 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. It is generally considered a better algorithm than the previous one, because it uses less costier operations (multiplications and additions instead of modulos.)
The perl program that implements that algorithm is this:
#### SNIP SNIP ##### #!/usr/bin/perl
use strict;
sub primes { my $how_much = shift;
my (@array, $bound, $a, $b, @primes);
@array = (1) x $how_much;
$bound = int(sqrt($how_much))+1;
for($a=2;$a<=$bound;$a++) { if ($array[$a]) { for($b=$a*$a;$b<$how_much;$b+=$a) { $array[$b] = 0; } push @primes, $a; } } for(;$a<$how_much;$a++) { if ($array[$a]) { push @primes, $a; } }
return @primes; }
print join(", ", primes(1000)); ##### SNIP SNIP ######
Now, I tried writing an equivalent Haskell program and the best I could do was the following:
---- SNIP SNIP ----- module Primes where
import Prelude import Array
how_much :: Int how_much = 1000
initial_primes_map :: Array Int Bool initial_primes_map = array (1, how_much) [ (i,True) | i <- [1 .. how_much] ]
mybound :: Int mybound = ceiling(sqrt(fromInteger(toInteger(how_much))))
next_primes_map :: Int -> Array Int Bool -> Array Int Bool next_primes_map a primes_map = if (a == mybound) then primes_map else next_primes_map (a+1) ( if primes_map!a then primes_map // [ (i*a, False) | i <- [a .. (prime_bound a)] ] else primes_map )
prime_bound :: Int -> Int prime_bound a = (floor(fromInteger(toInteger(how_much))/fromInteger(toInteger(a))))
get_primes_map :: Array Int Bool get_primes_map = (next_primes_map 2 initial_primes_map)
list_primes :: Array Int Bool -> Int -> [Int] list_primes primes_map n = if (n > how_much) then [] else ( if primes_map!n then n:(list_primes primes_map (n+1)) else list_primes primes_map (n+1) )
show_primes = show (list_primes get_primes_map 2) ---- SNIP SNIP -----
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
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. Does it have to do with the way hugs98 implements and Int to Bool array?
Regards,
Shlomi Fish
---------------------------------------------------------------------- Shlomi Fish shlomif@vipe.technion.ac.il Home Page: http://t2.technion.ac.il/~shlomif/ Home E-mail: shlomif@techie.com
The prefix "God Said" has the extraordinary logical property of converting any statement that follows it into a true one.
_______________________________________________ Haskell mailing list Haskell@haskell.org http://www.haskell.org/mailman/listinfo/haskell
--- Elke Kasimir Skalitzer Str. 79 10997 Berlin (Germany) fon: +49 (030) 612 852 16 mail: elke.kasimir@catmint.de> see: <http://www.catmint.de/elke> for pgp public key see: <http://www.catmint.de/elke/pgp_signature.html>