Geek asks and answers

Prime Sieve at Haskell

I am very new to Haskell, and I'm just trying to find the sum of the first 2 million primes. I am trying to generate prime numbers using a sieve (I think Eratosthenes sieve?), But it is really very slow, and I don’t know why. Here is my code.

sieve (x:xs) = x:(sieve $ filter (\a -> a `mod` x /= 0) xs)
ans = sum $ takeWhile (<2000000) (sieve [2..])

Thanks in advance.

+5
source share
5 answers

This is very slow because the algorithm is a test division that does not stop at the square root.

If you carefully look at what the algorithm does, you will see that for each prime pits multiples that do not have smaller prime divisors are removed from the list of candidates (multiples with smaller simply divisors were removed earlier).

, , , .

, , n √n.

, , k th k-1 . m n, , ,

(1-1) + (2-1) + (3-1) + ... + (m-1) = m*(m-1)/2

. n n / log n ( log ). n * √n, n, , .

, , 10 10 . , .

, ,

isPrime n = go 2
  where
    go d
      | d*d > n        = True
      | n `rem` d == 0 = False
      | otherwise      = go (d+1)

primes = filter isPrime [2 .. ]

1,9 * 10 9 ( , isPrime n √n - 179492732, ) ( 1) . , ( 2) -, .

Eratosthenes - O(n * log (log n)), :

primeSum.hs:

module Main (main) where

import System.Environment (getArgs)
import Math.NumberTheory.Primes

main :: IO ()
main = do
    args <- getArgs
    let lim = case args of
                (a:_) -> read a
                _     -> 1000000
    print . sum $ takeWhile (<= lim) primes

10 :

$ ghc -O2 primeSum && time ./primeSum 10000000
[1 of 1] Compiling Main             ( primeSum.hs, primeSum.o )
Linking primeSum ...
3203324994356

real    0m0.085s
user    0m0.084s
sys     0m0.000s

1 ( Int):

$ ghc -O2 tdprimeSum && time ./tdprimeSum 1000000
[1 of 1] Compiling Main             ( tdprimeSum.hs, tdprimeSum.o )
Linking tdprimeSum ...
37550402023

real    0m0.768s
user    0m0.765s
sys     0m0.002s

Turner 100 000:

$ ghc -O2 tuprimeSum && time ./tuprimeSum 100000
[1 of 1] Compiling Main             ( tuprimeSum.hs, tuprimeSum.o )
Linking tuprimeSum ...
454396537

real    0m2.712s
user    0m2.703s
sys     0m0.005s

(1)

2000000
   ∑ √k ≈ 4/3*√2*10^9
 k = 1

. - - , .

:

   ∑ √p ≈ 2/3*N^1.5/log N
 p < N
p prime

N = 2000000 1,3 * 10 8. , ( 1 n 2 N > 10).

, , () √k , , , .

, , ,

N^1.5/(log N)^2

n . , , - .

+12

Eratoshenes:

P= {3,5,...}\ & bigcup; {{p 2 p 2 + 2p,...} | p P}

( 2).:) "", Haskell ,

primes = 2 : g (fix g)  where
   g xs = 3 : (gaps 5 $ unionAll [[p*p, p*p+2*p..] | p <- xs])

unionAll ((x:xs):t) = x : union xs (unionAll $ pairs t)  where
  pairs ((x:xs):ys:t) = (x : union xs ys) : pairs t 
fix g = xs where xs = g xs
union (x:xs) (y:ys) = case compare x y of LT -> x : union  xs (y:ys)
                                          EQ -> x : union  xs    ys 
                                          GT -> y : union (x:xs) ys
gaps k s@(x:xs) | k<x  = k:gaps (k+2) s    
                | True =   gaps (k+2) xs

augustss 1,9 200k 2,1 400k, O(n^1.12..1.15) vs O(n^1.4) . 2,6 1 . .

, . , .

sieve (x:xs) = x:sieve [y|y<-xs, rem y p/=0] , :

[2..] ==> sieve --> 2
[3..] ==> nomult 2 ==> sieve --> 3
[4..] ==> nomult 2 ==> nomult 3 ==> sieve 
[5..] ==> nomult 2 ==> nomult 3 ==> sieve --> 5
[6..] ==> nomult 2 ==> nomult 3 ==> nomult 5 ==> sieve 
[7..] ==> nomult 2 ==> nomult 3 ==> nomult 5 ==> sieve --> 7
[8..] ==> nomult 2 ==> nomult 3 ==> nomult 5 ==> nomult 7 ==> sieve

nomult p = filter (\y->rem y p/=0). 8 3, , 3^2 == 9, 5 7.

, , . .

+5

, , ; ( ). :

import Control.Monad (forM_, when)
import Control.Monad.ST
import Data.Array.ST
import Data.Array.Unboxed

sieve :: Int -> UArray Int Bool
sieve n = runSTUArray $ do
    let m = (n-1) `div` 2
        r = floor . sqrt $ fromIntegral n
    bits <- newArray (0, m-1) True
    forM_ [0 .. r `div` 2 - 1] $ \i -> do
        isPrime <- readArray bits i
        when isPrime $ do
            forM_ [2*i*i+6*i+3, 2*i*i+8*i+6 .. (m-1)] $ \j -> do
                writeArray bits j False
    return bits

primes :: Int -> [Int]
primes n = 2 : [2*i+3 | (i, True) <- assocs $ sieve n]

http://ideone.com/mu1RN.

+3

primes :: [Integer]
primes = 2 : filter (isPrime primes) [3,5..]
  where isPrime (p:ps) n = p*p > n || n `rem` p /= 0 && isPrime ps n

. , . ( .)

+3

, , , , .

... .. n/log (n) 1 n. , 2 32 . 2 , . , . O (n ^ 2). O (n * (log n) * log (log n))

, . http://en.literateprograms.org/Sieve_of_Eratosthenes_ (Haskell)

+1

More articles:


All Articles