如何优化这个Haskell代码总结次线性时间的素数？

Question

来自Project Euler的问题10 是找到给定n下面所有素数的总和.

我只是通过总结Eratosthenes筛子产生的素数来解决它.然后我通过Lucy_Hedgehog(次线性!)找到了更有效的解决方案.

对于n =2⋅10^ 9:

• Python代码(来自上面的引用)在Python 2.7.3中运行1.2秒.

• C++代码(我的)在大约0.3秒内运行(用g ++ 4.8.4编译).

我在Haskell中重新实现了相同的算法,因为我正在学习它:

import Data.List

import Data.Map (Map, (!))
import qualified Data.Map as Map

problem10 :: Integer -> Integer
problem10 n = (sieve (Map.fromList [(i, i * (i + 1) `div` 2 - 1) | i <- vs]) 2 r vs) ! n
              where vs = [n `div` i | i <- [1..r]] ++ reverse [1..n `div` r - 1]
                    r  = floor (sqrt (fromIntegral n))

sieve :: Map Integer Integer -> Integer -> Integer -> [Integer] -> Map Integer Integer
sieve m p r vs | p > r     = m
               | otherwise = sieve (if m ! p > m ! (p - 1) then update m vs p else m) (p + 1) r vs

update :: Map Integer Integer -> [Integer] -> Integer -> Map Integer Integer
update m vs p = foldl' decrease m (map (\v -> (v, sumOfSieved m v p)) (takeWhile (>= p*p) vs))

decrease :: Map Integer Integer -> (Integer, Integer) -> Map Integer Integer
decrease m (k, v) = Map.insertWith (flip (-)) k v m

sumOfSieved :: Map Integer Integer -> Integer -> Integer -> Integer
sumOfSieved m v p = p * (m ! (v `div` p) - m ! (p - 1))

main = print $ problem10 $ 2*10^9

我用它编译ghc -O2 10.hs并运行time ./10.

它给出了正确答案,但大约需要7秒钟.

我用它编译ghc -prof -fprof-auto -rtsopts 10并运行./10 +RTS -p -h.

10.prof显示decrease占用52.2%的时间和67.5%的分配.

运行后hp2ps 10.hp我得到了这样的堆配置文件:

看起来像是decrease占用了大部分堆.GHC版本7.6.3.

您将如何优化此Haskell代码的运行时间？

---

更新13.06.17:

我试着更换不变Data.Map与可变Data.HashTable.IO.BasicHashTable从hashtables包,但我可能做坏事,因为对于微小的N = 30,它已经花费太长时间,大约10秒.怎么了？

更新18.06.17:

对HashTable性能问题感到好奇是一个很好的阅读.我使用了可变的Sherh的代码Data.HashTable.ST.Linear,但却放弃Data.Judy了.它在1.1秒内运行,仍然相对较慢.

Answer 1

我做了一些小的改进,所以它3.4-3.5在我的机器上运行几秒钟.使用IntMap.Strict帮助了很多.除此之外,我ghc只是为了确保手动执行一些优化.并使Haskell代码更接近链接中的Python代码.作为下一步,您可以尝试使用一些可变的HashMap.但我不确定...... IntMap不能比一些可变容器快得多,因为它是一个不可变的容器.虽然我对它的效率感到惊讶.我希望这可以更快地实施.

这是代码:

import Data.List (foldl')
import Data.IntMap.Strict (IntMap, (!))
import qualified Data.IntMap.Strict as IntMap

p :: Int -> Int
p n = (sieve (IntMap.fromList [(i, i * (i + 1) `div` 2 - 1) | i <- vs]) 2 r vs) ! n
               where vs = [n `div` i | i <- [1..r]] ++ [n', n' - 1 .. 1]
                     r  = floor (sqrt (fromIntegral n) :: Double)
                     n' = n `div` r - 1

sieve :: IntMap Int -> Int -> Int -> [Int] -> IntMap Int
sieve m' p' r vs = go m' p'
  where
    go m p | p > r               = m
           | m ! p > m ! (p - 1) = go (update m vs p) (p + 1)
           | otherwise           = go m (p + 1)

update :: IntMap Int -> [Int] -> Int -> IntMap Int
update s vs p = foldl' decrease s (takeWhile (>= p2) vs)
  where
    sp = s ! (p - 1)
    p2 = p * p
    sumOfSieved v = p * (s ! (v `div` p) - sp)
    decrease m  v = IntMap.adjust (subtract $ sumOfSieved v) v m

main :: IO ()
main = print $ p $ 2*10^(9 :: Int) 

更新:

使用mutable hashtables我已经设法~5.5sec通过这个实现在Haskell上实现性能.

此外,我在几个地方使用了未装箱的向量而不是列表.Linear哈希似乎是最快的.我认为这可以做得更快.我注意到sse42选项的hasthables包.不确定我是否设法正确设置它,但即使没有它运行得那么快.

更新2(19.06.2017)

我已经设法3x通过删除judy hashmap 来使它更快,然后从@Krom(使用我的代码+他的地图)获得最佳解决方案.相反,只使用普通数组.你可以拿出同样的想法,如果你发现钥匙SHashMap的是无论从顺序1到n'或n div i为i从1到r.因此,我们可以将这样的HashMap表示为两个数组,这些数组根据搜索关键字在数组中进行查找.

我的代码+ Judy HashMap

$ time ./judy
95673602693282040

real    0m0.590s
user    0m0.588s
sys     0m0.000s

我的代码+我的稀疏地图

$ time ./sparse
95673602693282040

real    0m0.203s
user    0m0.196s
sys     0m0.004s

如果不使用IOUArray已生成的向量和Vector库并readArray替换为,则可以更快地完成此操作unsafeRead.但是,如果只是你真的没有兴趣尽可能地优化它,我认为不应该这样做.

与此解决方案的比较是作弊并且不公平.我希望在Python和C++中实现的相同想法会更快.但是具有封闭散列图的@Krom解决方案已经在作弊,因为它使用自定义数据结构而不是标准数据结构.至少你可以看到Haskell中的标准和最流行的哈希映射并不那么快.使用更好的算法和更好的ad-hoc数据结构可以更好地解决这些问题.

这是结果代码.