如前所述，使用可变数组具有最佳性能。以下代码源自此“TemplateHaskell”版本，该版本转换回更符合直接可变数组解决方案的内容，因为“TemplateHaskell”似乎没有任何区别，并进行了一些进一步的优化。我相信它比通常的可变未装箱数组版本更快，因为进一步优化，特别是由于使用了“unsafeRead”和“unsafeWrite”函数，避免了数组的范围检查，可能也在内部使用指针进行数组访问：
{-# OPTIONS -O2 -optc-O3 #-}
import Control.Monad
import Control.Monad.ST
import Data.Array.ST
import Data.Array.Unboxed
import Data.Array.Base
primesToUA :: Word32-> [Word32]
primesToUA top = do
let sieveUA top = runSTUArray $ do
let m = ((fromIntegral top) - 3) `div` 2 :: Int
buf <- newArray (0,m) True -- :: ST s (STUArray s Int Bool)
let cullUA i = do
let p = i + i + 3
strt = p * (i + 1) + i
let cull j = do
if j > m then cullUA (i + 1)
else do
unsafeWrite buf j False
cull (j + p)
if strt > m then return ()
else do
e <- unsafeRead buf i
if e then cull strt else cullUA (i + 1)
cullUA 0; return buf
if top > 1 then 2 : [2 * (fromIntegral i) + 3 | (i,True) <- assocs $ sieveUA top]
else []
main = do
x <- read `fmap` getLine -- 1mln 2mln 10mln 100mln
print (last (primesToUA x)) -- 0.01 0.02 0.09 1.26 seconds
编辑：上面的代码已被更正，下面的时间被编辑以反映更正，以及将分页与非分页版本进行比较的评论已被编辑。
将其运行到指示的最高范围的时间如上述源代码底部的评论表所示，由ideone.com测量，并且比@WillNess 发布的答案快大约五倍，也是在ideone测量的.com. 这段代码将素数剔除到 200 万只需要很短的时间，而将素数剔除到一亿只需要 1.26 秒。在我的 i7 (3.5 GHz) 上以 0.44 秒运行到一亿时，这些时间快了大约 2.86 倍，运行到十亿需要 6.81 秒。前者的内存使用量刚刚超过 6 兆字节，后者的内存使用量为 60 兆字节，这是巨大（位压缩）数组使用的内存。该数组还解释了非线性性能，因为当数组大小超过 CPU 缓存大小时，每个复合数表示剔除的平均内存访问时间会变差。
编辑_添加： 页面分段筛更有效，因为当缓冲区大小保持小于 L1 或 L2 CPU 缓存时，它具有更好的内存访问效率，并且具有无界的优点，因此上限不必是预先指定的内存占用量要小得多，即基素数小于所用范围的平方根加上页面缓冲内存。以下代码已编写为分页实现，比非分页版本要快一些；它还提供了一个优点，即可以将代码顶部的输出范围规范从“Word32”更改为“Word64”，因此它不限于 32 位数字范围，只需要花费一点处理时间（对于 32 位编译代码）适用于任何通用范围。代码如下：
-- from http://www.haskell.org/haskellwiki/Prime_numbers#Using_ST_Array
{-# OPTIONS -O2 -optc-O3 #-}
import Data.Word
import Control.Monad
import Control.Monad.ST
import Data.Array.ST
import Data.Array.Unboxed
import Data.Array.Base
primesUA :: () -> [Word32]
primesUA () = do
let pgSZBTS = 262144 * 8
let sieveUA low bps = runSTUArray $ do
let nxt = (fromIntegral low) + (fromIntegral pgSZBTS)
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
when (strt < pgSZBTS) $ do
e <- unsafeRead buf i
if e then do
let loop j = do
if j < pgSZBTS then do
unsafeWrite buf j False
loop (j + p)
else cullUAbase (i + 1)
loop strt
else cullUAbase (i + 1)
let cullUA ~(p:t) = do
let bp = (fromIntegral p)
i = (bp - 3) `div` 2
s = bp * (i + 1) + i
when (s < nxt) $ do
let strt = do
if s >= low then fromIntegral (s - low)
else do
let b = (low - s) `rem` bp
if b == 0 then 0 else fromIntegral (bp - b)
let loop j = do
if j < pgSZBTS then do
unsafeWrite buf j False
loop (j + (fromIntegral p))
else cullUA t
loop strt
if low <= 0 then cullUAbase 0 else cullUA bps
return buf
let sieveList low bps = do
[2 * ((fromIntegral i) + low) + 3 | (i,True) <- assocs $ 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 ())))
由于需要以不同于后续页面的方式从第一页数组中剔除复合数表示，因此上面的代码比非分页的情况多出几行代码。该代码还具有修复程序，因此不会由于基本素数列表和输出列表现在不是同一个列表而导致内存泄漏（从而避免将整个列表保留在内存中）。
请注意，由于剔除缓冲区的大小恒定小于 L2 CPU 缓存，此代码需要接近线性时间（在筛选范围内），因为范围变得更大。内存使用量是非分页版本使用量的一小部分，1 亿只使用不到 600 KB，10 亿只使用超过 600 KB，略有增加只是基素数所需的额外空间小于平方根的范围列表。
在ideone.com 上，此代码在大约 1.13 秒内生成一亿个素数，大约 12 秒到十亿个（32 位设置）。可能车轮分解和肯定的多核处理会使它在多核 CPU 上更快。在我的 i7 (3.5 GHz) 上，筛分到一亿只需要 0.44 秒，筛到十亿只需要 4.7 秒，随着范围的增加，性能大致呈线性。似乎在 ideone.com 上运行的 GHC 版本中存在某种非线性开销，对于 i7 不存在的较大范围有一些性能损失，这可能与不同的垃圾收集有关，如页面正在为每个新页面创建新的缓冲区。 END_EDIT_ADD
EDIT_ADD2：似乎上述页面分段代码的大部分处理时间都用于（惰性）列表处理，因此相应地重新编写了代码，并进行了如下几个改进：
实现了一个不使用列表处理并使用“popCount”表查找来一次计算 32 位字中“一个”位的数量的素数计数功能。这样，与实际筛分时间相比，找到结果的时间是微不足道的。
将基本素数存储为位打包页面段列表，这比存储素数列表更节省空间，并且根据需要将页面段转换为素数的时间并不是很大的计算开销。
调整素数段生成函数，使得对于初始零页段，它使用自己的位模式作为源页，从而使复合数剔除代码更短更简单。
然后代码变成如下：
{-# OPTIONS -O3 -rtsopts #-} -- -fllvm ide.com doesn't support LLVM
import Data.Word
import Data.Bits
import Control.Monad
import Control.Monad.ST
import Data.Array.ST (runSTUArray)
import Data.Array.Unboxed
import Data.Array.Base
pgSZBTS = (2^18) * 8 :: Int -- size of L2 data cache
type PrimeType = Word32
type Chunk = UArray PrimeType Bool
-- makes a new page chunk and culls it
-- if the base primes list provided is empty then
-- it uses the current array as source (for zero page base primes)
mkChnk :: Word32 -> [Chunk] -> Chunk
mkChnk low bschnks = runSTUArray $ do
let nxt = (fromIntegral low) + (fromIntegral pgSZBTS)
buf <- nxt `seq` newArray (fromIntegral low, fromIntegral nxt - 1) True
let cull ~(p:ps) =
let bp = (fromIntegral p)
i = (bp - 3) `shiftR` 1
s = bp * (i + 1) + i in
let cullp j = do
if j >= pgSZBTS then cull ps
else do
unsafeWrite buf j False
cullp (j + (fromIntegral p)) in
when (s < nxt) $ do
let strt = do
if s >= low then fromIntegral (s - low)
else do
let b = (low - s) `rem` bp
if b == 0 then 0 else fromIntegral (bp - b)
cullp strt
case bschnks of
[] -> do bsbf <- unsafeFreezeSTUArray buf
cull (listChnkPrms [bsbf])
_ -> cull $ listChnkPrms bschnks
return buf
-- creates a page chunk list starting at the lw value
chnksList :: Word32 -> [Chunk]
chnksList lw =
mkChnk lw basePrmChnks : chnksList (lw + fromIntegral pgSZBTS)
-- converts a page chunk list to a list of primes
listChnkPrms :: [Chunk] -> [PrimeType]
listChnkPrms [] = []
listChnkPrms ~(hdchnk@(UArray lw _ rng _):tlchnks) =
let nxtp i =
if i >= rng then [] else
if unsafeAt hdchnk i then
(case ((lw + fromIntegral i) `shiftL` 1) + 3 of
np -> np) : nxtp (i + 1)
else nxtp (i + 1) in
(hdchnk `seq` lw `seq` nxtp 0) ++ listChnkPrms tlchnks
-- the base page chunk list used to cull the higher order pages,
-- note that it has special treatment for the zero page.
-- It is more space efficient to store this as chunks rather than
-- as a list of primes or even a list of deltas (gaps), with the
-- extra processing to convert as needed not too much.
basePrmChnks :: [Chunk]
basePrmChnks = mkChnk 0 [] : chnksList (fromIntegral pgSZBTS)
-- the full list of primes could be accessed with the following function.
primes :: () -> [PrimeType]
primes () = 2 : (listChnkPrms $ chnksList 0)
-- a quite fast prime counting up to the given limit using
-- chunk processing to avoid innermost list processing loops.
countPrimesTo :: PrimeType -> Int
countPrimesTo limit =
let lmtb = (limit - 3) `div` 2 in
let sumChnks acc chnks@(chnk@(UArray lo hi rng _):chnks') =
let cnt :: UArray PrimeType Word32 -> Int
cnt bfw =
case if lmtb < hi then fromIntegral (lmtb - lo) else rng of
crng -> case crng `shiftR` 5 of
rngw ->
let cnt' i ac =
ac `seq` if i >= rngw then
if (i `shiftL` 5) >= rng then ac else
case (-2) `shiftL` fromIntegral (lmtb .&. 31) of
msk -> msk `seq`
case (unsafeAt bfw rngw) .&.
(complement msk) of
bts -> bts `seq` case popCount bts of
c -> c `seq` case ac + c of nac -> nac
else case ac + (popCount $ unsafeAt bfw i) of
nacc -> nacc `seq` cnt' (i + 1) (nacc)
in cnt' 0 0 in
acc `seq` case runST $ do -- make UArray _ Bool into a UArray _ Word32
stbuf <- unsafeThawSTUArray chnk
stbufw <- castSTUArray stbuf
bufw <- unsafeFreezeSTUArray stbufw
return $ cnt bufw of
c -> c `seq` case acc + c of
nacc -> nacc `seq` if hi >= lmtb then nacc
else sumChnks nacc chnks' in
if limit < 2 then 0 else if limit < 3 then 1 else
lmtb `seq` sumChnks 1 (chnksList 0)
main = do
x <- read `fmap` getLine -- 1mln 2mln 10mln 100mln 1000mln
-- 0.02 0.03 0.06 0.45 4.60 seconds
-- 7328 7328 8352 8352 9424 Kilobytes
-- this takes 14.34 seconds and 9424 Kilobytes to 3 billion on ideone.com,
-- and 9.12 seconds for 3 billion on an i7-2700K (3.5 GHz).
-- The above ratio of about 1.6 is much better than usual due to
-- the extremely low memory use of the page segmented algorithm.
-- It seems thaat the Windows Native Code Generator (NCG) from GHC
-- is particularly poor, as the Linux 32-bit version takes
-- less than two thirds of the time for exactly the same program...
print $ countPrimesTo x
-- print $ length $ takeWhile ((>=) x) $ primes () -- the slow way to do this
代码中给出的时间和内存要求与在 ideone.com 上运行时观察到的一样，在我的 i7-2700K (3.5 GHz) 上运行一、二、十、一百、一千（十亿）和三千（三十亿）需要 0.02、0.03、0.05、0.30、3.0 和 9.1 秒) 百万范围，具有几乎恒定的内存占用量，随着基本素数的数量小于所需范围的平方根而缓慢增加。当使用 LLVM 编译器后端编译时，这些时间分别变为 0.01、0.02、0.02、0.12、1.35 和 4.15 秒，因为它更有效地使用了寄存器和机器指令；最后一个非常接近与使用 64 位编译器而不是使用 32 位编译器编译时相同的速度，因为寄存器的有效使用意味着额外寄存器的可用性没有太大区别。
正如代码中所评论的那样，由于没有受到内存访问瓶颈的限制，我的真实机器和 ideone.com 服务器上的性能之间的比率变得比内存浪费更多的算法要小得多，因此速度限制主要是比率CPU 时钟速度以及每个时钟周期的 CPU 处理效率。然而，正如那里所评论的那样，当在 Windows（32 位编译器）下运行时，GHC 本地代码生成器（NCG）存在一些奇怪的低效率，因为运行时间比在 Linux 下运行（作为 ideone. com 服务器使用）。AFAIK 他们都有相同的 Haskell GHC 版本的通用代码库，唯一的分歧在于使用的链接器（它也与 LLVM 后端一起使用，不受影响），因为 GHC NCG 不使用 GCC，而只使用 mingw32 汇编器,
请注意，当使用 LLVM 编译器后端编译此代码时，其速度与为高度优化的“C/C++”实现编写的相同算法的速度大致相同，这表明 Haskell 确实有能力开发非常紧密的循环编码。可以说，一旦习惯了单子和非严格代码的 Haskell 范例，Haskell 代码比等效的“C/C++”代码更具可读性和安全性。埃拉托色尼筛子执行速度的进一步改进纯粹是对所用实现的调整，而不是在 Haskell 和“C/C++”之间选择语言。
总结： 当然，这还不是 Haskell 版本的埃拉托色尼筛法的终极速度，因为我们还没有调整内存访问以更有效地使用快速 CPU L1 缓存，也没有显着降低除了消除赔率处理之外，使用极端车轮分解所需的复合剔除操作的总数。然而，这足以回答这个问题，即可变数组是解决此类紧密循环类型问题的最有效方法，与使用列表或不可变数组相比，潜在的速度增益约为 100 倍。 END_EDIT_ADD2