-- from http://w...content-available-to-author-only...l.org/haskellwiki/Prime_numbers#Using_ST_Array
{-# OPTIONS -O2 -optc-O3 #-}
import Data.Word
import Data.Array.ST
import Data.Array.Unboxed
import Data.Array.Base
primesUA :: () -> [Word32]
primesUA () = do
let pgSZBTS = 262144 * 8
sieveUA low bps = do
buf
<- newArray
(0,pgSZBTS
- 1) True
:: ST s
(STUArray s
Int Bool) let cullUAbase i = do
let p = i + i + 3
strt = p * (i + 1) + i
let cull j = do
if j >= pgSZBTS then cullUAbase (i + 1)
else do
unsafeWrite buf j False
cull (j + p)
when (strt < pgSZBTS) $ do
e <- unsafeRead buf i
if e then cull strt else cullUAbase (i + 1)
let cullUA ~(p:t) = do
s = bp * (i + 1) + i
let cull j = do
if j >= pgSZBTS then cullUA t
else do
unsafeWrite buf j False
when (s < nxt) $ do
let strt = do
else do
let b
= (low
- s
) `
rem` bp
cull strt
if low <= 0 then cullUAbase 0 else cullUA bps
let sieveList low bps = do
[2 * ((fromIntegral i
) + low
) + 3 | (i
,True
) <- assocs
$ runSTUArray
$ sieveUA low bps
] let sieve low bps = do
(sieveList low bps
) ++ sieve
(low
+ (fromIntegral pgSZBTS
)) bps
let primes' = ((sieveList 0 []) ++ sieve (fromIntegral pgSZBTS) primes') :: [Word32]
2 : sieve 0 primes'
main = do
x <- read `fmap` getLine -- 1mln 2mln 10mln 100mln
-- 0.02 0.03 0.13 1.13 seconds
print (length (takeWhile ((>=) (fromIntegral x)) (primesUA ())))
LS0gZnJvbSBodHRwOi8vdy4uLmNvbnRlbnQtYXZhaWxhYmxlLXRvLWF1dGhvci1vbmx5Li4ubC5vcmcvaGFza2VsbHdpa2kvUHJpbWVfbnVtYmVycyNVc2luZ19TVF9BcnJheQoKey0jIE9QVElPTlMgLU8yIC1vcHRjLU8zICMtfQppbXBvcnQgRGF0YS5Xb3JkCmltcG9ydCBDb250cm9sLk1vbmFkCmltcG9ydCBDb250cm9sLk1vbmFkLlNUCmltcG9ydCBEYXRhLkFycmF5LlNUCmltcG9ydCBEYXRhLkFycmF5LlVuYm94ZWQKaW1wb3J0IERhdGEuQXJyYXkuQmFzZQoKcHJpbWVzVUEgOjogKCkgLT4gW1dvcmQzMl0KcHJpbWVzVUEgKCkgPSBkbwogICAgbGV0IHBnU1pCVFMgPSAyNjIxNDQgKiA4CiAgICBsZXQgc2lldmVVQSA6OiAoSW50ZWdyYWwgdCwgSW50ZWdyYWwgdDEpID0+IHQgLT4gW3QxXSAtPiBTVCBzIChTVFVBcnJheSBzIEludCBCb29sKQogICAgICAgIHNpZXZlVUEgbG93IGJwcyA9IGRvCiAgICAgICAgICAgIGxldCBueHQgPSAoZnJvbUludGVncmFsIGxvdykgKyAoZnJvbUludGVncmFsIHBnU1pCVFMpCiAgICAgICAgICAgIGJ1ZiA8LSBuZXdBcnJheSAoMCxwZ1NaQlRTIC0gMSkgVHJ1ZSA6OiBTVCBzIChTVFVBcnJheSBzIEludCBCb29sKQogICAgICAgICAgICBsZXQgY3VsbFVBYmFzZSBpID0gZG8KICAgICAgICAgICAgICAgIGxldCBwID0gaSArIGkgKyAzCiAgICAgICAgICAgICAgICAgICAgc3RydCA9IHAgKiAoaSArIDEpICsgaQogICAgICAgICAgICAgICAgbGV0IGN1bGwgaiA9IGRvCiAgICAgICAgICAgICAgICAgICAgaWYgaiA+PSBwZ1NaQlRTIHRoZW4gY3VsbFVBYmFzZSAoaSArIDEpCiAgICAgICAgICAgICAgICAgICAgZWxzZSBkbwogICAgICAgICAgICAgICAgICAgICAgICB1bnNhZmVXcml0ZSBidWYgaiBGYWxzZQogICAgICAgICAgICAgICAgICAgICAgICBjdWxsIChqICsgcCkKICAgICAgICAgICAgICAgIHdoZW4gKHN0cnQgPCBwZ1NaQlRTKSAkIGRvCiAgICAgICAgICAgICAgICAgICAgZSA8LSB1bnNhZmVSZWFkIGJ1ZiBpCiAgICAgICAgICAgICAgICAgICAgaWYgZSB0aGVuIGN1bGwgc3RydCBlbHNlIGN1bGxVQWJhc2UgKGkgKyAxKQogICAgICAgICAgICBsZXQgY3VsbFVBIH4ocDp0KSA9IGRvCiAgICAgICAgICAgICAgICBsZXQgYnAgPSAoZnJvbUludGVncmFsIHApCiAgICAgICAgICAgICAgICAgICAgaSA9IChicCAtIDMpIGBkaXZgIDIKICAgICAgICAgICAgICAgICAgICBzID0gYnAgKiAoaSArIDEpICsgaQogICAgICAgICAgICAgICAgbGV0IGN1bGwgaiA9IGRvCiAgICAgICAgICAgICAgICAgICAgaWYgaiA+PSBwZ1NaQlRTIHRoZW4gY3VsbFVBIHQKICAgICAgICAgICAgICAgICAgICBlbHNlIGRvCiAgICAgICAgICAgICAgICAgICAgICAgIHVuc2FmZVdyaXRlIGJ1ZiBqIEZhbHNlCiAgICAgICAgICAgICAgICAgICAgICAgIGN1bGwgKGogKyAoZnJvbUludGVncmFsIHApKQogICAgICAgICAgICAgICAgd2hlbiAocyA8IG54dCkgJCBkbwogICAgICAgICAgICAgICAgICAgIGxldCBzdHJ0ID0gZG8KICAgICAgICAgICAgICAgICAgICAgICAgaWYgcyA+PSBsb3cgdGhlbiBmcm9tSW50ZWdyYWwgKHMgLSBsb3cpCiAgICAgICAgICAgICAgICAgICAgICAgIGVsc2UgIGRvCiAgICAgICAgICAgICAgICAgICAgICAgICAgICBsZXQgYiA9IChsb3cgLSBzKSBgcmVtYCBicAogICAgICAgICAgICAgICAgICAgICAgICAgICAgaWYgYiA9PSAwIHRoZW4gMCBlbHNlIGZyb21JbnRlZ3JhbCAoYnAgLSBiKQogICAgICAgICAgICAgICAgICAgIGN1bGwgc3RydAogICAgICAgICAgICBpZiBsb3cgPD0gMCB0aGVuIGN1bGxVQWJhc2UgMCBlbHNlIGN1bGxVQSBicHMKICAgICAgICAgICAgcmV0dXJuIGJ1ZgogICAgbGV0IHNpZXZlTGlzdCBsb3cgYnBzID0gZG8KICAgICAgICBbMiAqICgoZnJvbUludGVncmFsIGkpICsgbG93KSArIDMgfCAoaSxUcnVlKSA8LSBhc3NvY3MgJCBydW5TVFVBcnJheSAkIHNpZXZlVUEgbG93IGJwc10KICAgIGxldCBzaWV2ZSBsb3cgYnBzID0gZG8KICAgICAgICAoc2lldmVMaXN0IGxvdyBicHMpICsrIHNpZXZlIChsb3cgKyAoZnJvbUludGVncmFsIHBnU1pCVFMpKSBicHMKICAgIGxldCBwcmltZXMnID0gKChzaWV2ZUxpc3QgMCBbXSkgKysgc2lldmUgKGZyb21JbnRlZ3JhbCBwZ1NaQlRTKSBwcmltZXMnKSA6OiBbV29yZDMyXQogICAgMiA6IHNpZXZlIDAgcHJpbWVzJwoKbWFpbiA9IGRvIAogICB4IDwtIHJlYWQgYGZtYXBgIGdldExpbmUgICAgICAtLSAgIDFtbG4gICAgMm1sbiAgICAxMG1sbiAgIDEwMG1sbgogICAgICAgICAgICAgICAgICAgICAgICAgICAgICAgICAtLSAgIDAuMDIgICAgMC4wMyAgICAgMC4xMyAgICAgMS4xMyAgc2Vjb25kcwogICBwcmludCAobGVuZ3RoICh0YWtlV2hpbGUgKCg+PSkgKGZyb21JbnRlZ3JhbCB4KSkgKHByaW1lc1VBICgpKSkp