RE: Finding primes using a primes map with Haskell and Hugs98
| 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..]] This style is definitely the way to go. Haskell does badly if you update an array one index at a time. Remember that arrays can be recursive. Here's a definition of Fibonacci for example; you can probably adapt it for primes fibs :: Int -> Array Int Int -- If a = fibs n, then a!i is fib(i), for i<=n. fibs n = a where a = array (1,n) ([(1,1),(2,1)] ++ [(i,a!(i-1) + a!(i-2) | i <- [3..n]]) -- Notice that a is recursive Simon
On Tue, 19 Dec 2000, Simon Peyton-Jones wrote:
| 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..]]
This style is definitely the way to go. Haskell does badly if you update an array one index at a time.
Unfortunately, it seems that this style is not the way to go. This program cannot scale beyond 5000 while my second program scales beyond 30000. I'm not saying 30000 is a good limit, but 5000 is much worse. Anyway, somebody who contacted me in private suggested the following method. It is a similiar algorithm which uses a list instead of an array. primes :: Int -> [Int] primes how_much = sieve [2..how_much] where sieve (p:x) = p : (if p <= mybound then sieve (remove (p*p) x) else x) where remove what (a:as) | what > how_much = (a:as) | a < what = a:(remove what as) | a == what = (remove (what+step) as) | a > what = a:(remove (what+step) as) remove what [] = [] step = (if (p == 2) then p else (2*p)) sieve [] = [] mybound = ceiling(sqrt(fromIntegral how_much)) I optimized it quite a bit, but the concept remained the same. Anyway, this code can scale very well to 100000 and beyond. But it's not exactly the same algorithm. I also implemented this algorithm in perl, and I can send it in person if anybody requests it. I'll try to see how the two programs run in GHC and HBC. Regards, Shlomi Fish
Remember that arrays can be recursive. Here's a definition of Fibonacci for example; you can probably adapt it for primes
fibs :: Int -> Array Int Int -- If a = fibs n, then a!i is fib(i), for i<=n. fibs n = a where a = array (1,n) ([(1,1),(2,1)] ++ [(i,a!(i-1) + a!(i-2) | i <- [3..n]]) -- Notice that a is recursive
Simon
_______________________________________________ Haskell mailing list Haskell@haskell.org http://www.haskell.org/mailman/listinfo/haskell
---------------------------------------------------------------------- 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.
There are numerous ways of optimising sieving for primes, none of which have much to do with this list. For example, two suggestions: (1) for each k modulo 2*3*5*7, if k is divisible by 2/3/5 or 7, ignore, otherwise sieve separately for this k on higher primes. (Or you might use products of more or less primes, depending on memory and how high you were going.) (2) use bitwise arithmetic. If you look in the literature I think you'll find plenty more possibilities. I don't really see why any of this has anything to do with Haskell though. When it comes to seriously icky bit-twiddling algorithms I don't think Haskell has much to offer over C, especially as you'd have to make everything unboxed if you want comparable speed.
Hello! On Wed, Dec 20, 2000 at 04:02:23PM +0200, Shlomi Fish wrote:
[...]
primes :: Int -> [Int] primes how_much = sieve [2..how_much] where sieve (p:x) = p : (if p <= mybound then sieve (remove (p*p) x) else x) where remove what (a:as) | what > how_much = (a:as) | a < what = a:(remove what as) | a == what = (remove (what+step) as) | a > what = a:(remove (what+step) as) remove what [] = [] step = (if (p == 2) then p else (2*p)) sieve [] = [] mybound = ceiling(sqrt(fromIntegral how_much))
[...]
How about this: Yet another infinite list version, however only sieving numbers to the square of the current maximum: primes :: [Int] primes = 2:fil' 4 primes [3..] where fil' cutoff fl@(f:fs) rl@(r:rs) = if r < cutoff then r:fil' cutoff fl rs else fil' (square (head fs)) fs [y|y<-rl, (y `mod` f) /= 0] square x = x*x Comparison: my version (with a small main around it): hannah@mamba:~/src/haskell $ time ./primes 1000000 >/tmp/p1 real 0m31.928s user 0m31.690s sys 0m0.317s your version (with adapted small main around it): hannah@mamba:~/src/haskell $ time ./primes2 1000000 >/tmp/p2 real 0m42.993s user 0m42.023s sys 0m0.236s (measurements with ghc 4.08, -O2 -O2-for-C, Pentium I 200 MHz 64 MB RAM) The outputs (primes list + number of found primes) are the same. In contrast, take the simple two-liner primes :: [Int] primes = p' [2..] where p' (p:ps) = p:p' [y|y<-ps, (y `mod` p) /= 0] with the same small main: hannah@mamba:~/src/haskell $ time ./primes0 100000 >/tmp/p0 real 2m17.460s user 2m15.773s sys 0m1.170s (mine with 100000 limit) hannah@mamba:~/src/haskell $ time ./primes 100000 >/tmp/p1 real 0m1.666s user 0m1.607s sys 0m0.048s (yours with 100000 limit) hannah@mamba:~/src/haskell $ time ./primes2 100000 >/tmp/p2 real 0m1.903s user 0m1.827s sys 0m0.064s Kind regards, Hannah.
participants (4)
-
George Russell -
Hannah Schroeter -
Shlomi Fish -
Simon Peyton-Jones