Какова идеальная реализация решета Эратосфена между списками, массивами и изменяемыми массивами?

Этот

primes = sieve [2..] where
         sieve (p:x) = p : sieve [ n | n <- x, n `mod` p > 0 ]

это не сито. Это очень неэффективное пробное подразделение. Не используйте это!

Мне любопытно, как вы получили свое время, «решето» Тернера никак не может производить простые числа, не превышающие 2 000 000, за считанные секунды. Чтобы найти простые числа до 200 000, потребовалось

MUT     time    6.38s  (  6.39s elapsed)
GC      time    9.19s  (  9.20s elapsed)
EXIT    time    0.00s  (  0.00s elapsed)
Total   time   15.57s  ( 15.59s elapsed)

на моем компьютере (64-битный Linux, ghc-7.6.1, скомпилированный с -O2). Сложность этого алгоритма O(N² / log² N), почти квадратичная. Чтобы увеличить его до 2 000 000, потребуется около двадцати минут.

Ваши времена для версий массива тоже подозрительны, правда в другую сторону. Вы измеряли интерпретируемый код?

При просеивании до 2 000 000, скомпилированном с оптимизацией код изменяемого массива выполнялся за 0,35 секунды, а код неизменяемого массива — 0,12 секунды.

Теперь у него все еще есть изменяемый массив примерно в три раза медленнее, чем неизменный массив.

Но, это несправедливое сравнение. Для неизменяемого массива вы использовали Ints, а для изменяемого массива Integers. Изменение кода изменяемого массива для использования Ints — как и должно быть, поскольку под капотом массивы индексируются Int, поэтому использование Integer — это ненужная жертва производительности, которая ничего не покупает — заставило код изменяемого массива выполняться за 0,15 секунды. Близко к коду изменяемого массива, но не совсем там. Однако вы позволяете изменяемому массиву выполнять больше работы, поскольку в коде неизменяемого массива вы удаляете только нечетные числа, кратные нечетным простым числам, а в коде изменяемого массива вы помечаете все числа, кратные всем простым числам. Изменение кода изменяемого массива для специальной обработки 2 и устранения только нечетных кратных нечетных простых чисел также снижает это время до 0,12 секунды.

Но вы используете индексирование массива с проверкой диапазона, что медленно и, поскольку правильность индексов проверяется в самом коде, не нужно. Изменение этого на использование unsafeRead и unsafeWrite сокращает время для неизменяемого массива до 0,09 секунды.

Тогда у вас есть проблема, что использование

forM_ [x, y .. z]

использует упакованные Ints (к счастью, GHC исключает список). Написание цикла самостоятельно, чтобы использовались только неупакованные Int#s, время сократилось до 0,02 секунды.

{-# LANGUAGE MonoLocalBinds #-}
import Control.Monad (forM_, when)
import Control.Monad.ST
import Data.Array.ST
import Data.Array.Unboxed
import Data.Array.Base

primeSieve :: Int -> UArray Int Bool
primeSieve top = runSTUArray $ do
    a <- newArray (0,top) True
    unsafeWrite a 0 False
    unsafeWrite a 1 False
    let r = ceiling . sqrt $ fromIntegral top
        mark step idx
            | top < idx = return ()
            | otherwise = do
                unsafeWrite a idx False
                mark step (idx+step)
        sift p
            | r < p     = return a
            | otherwise = do
                prim <- unsafeRead a p
                when prim $ mark (2*p) (p*p)
                sift (p+2)
    mark 2 4
    sift 3

-- Return primes from sieve as list:
primesTo :: Int -> [Int]
primesTo top = [p | (p,True) <- assocs $ primeSieve top]

main :: IO ()
main = print .last $ primesTo 2000000

Итак, в заключение, для решета Эратосфена вы должны использовать массив — неудивительно, его эффективность зависит от возможности перехода от одного кратного к другому за короткое постоянное время.

Вы получаете очень простой и понятный код с неизменяемыми массивами, и этот код прилично работает при не слишком высоких ограничениях (при более высоких ограничениях ситуация становится относительно хуже, поскольку массивы все еще копируются и очищаются от мусора, но это не так уж плохо).

Когда вам нужна более высокая производительность, вам нужны изменяемые массивы. Написание эффективного кода изменяемого массива не совсем тривиально, нужно знать, как компилятор транслирует различные конструкции, чтобы выбрать правильный, и некоторые сочтут такой код однозначным. Но вы также можете использовать библиотеку (отказ от ответственности: я написал ее), которая обеспечивает довольно эффективную реализацию, а не писать ее самостоятельно.

Изменяемый массив всегда будет победителем с точки зрения производительности (и вам действительно нужно было скопировать версию, которая работает только с коэффициентами как минимум; он должен быть самым быстрым из трех — еще и потому, что он использует Int, а не Integer).

Для списков инкрементное сито в форме дерева должен работать лучше, чем тот, который вы показываете. Вы всегда можете использовать его с takeWhile (< limit), если это необходимо. Я утверждаю, что оно наиболее ясно передает истинную природу сита:

import Data.List (unfoldr)

primes :: [Int]         
primes = 2 : _Y ((3 :) . gaps 5 . _U . map (\p -> [p*p, p*p+2*p..]))

_Y g = g (_Y g)                                -- recursion 
_U ((x:xs):t) = (x :) . union xs . _U          -- ~= nub . sort . concat
                      . unfoldr (\(a:b:c) -> Just (union a b, c)) $ t

gaps k s@(x:xs) | k < x     = k : gaps (k+2) s   -- ~= [k,k+2..]\\s, when 
                | otherwise =     gaps (k+2) xs  --  k<=x && null(s\\[k,k+2..])

union a@(x:xs) b@(y:ys) = case compare x y of  -- ~= nub . sort .: (++)
         LT -> x : union  xs  b
         EQ -> x : union  xs ys
         GT -> y : union  a  ys

_U повторно реализует Data.List.Ordered.unionAll< /a>, а gaps 5 – это (minus [5,7..]), объединенное для эффективности с minus и union из одного и того же пакет.

Конечно, ничто не сравнится с краткостью Data.List.nubBy (((>1).).gcd) [2..] (но это очень медленно).

На ваш 1-й новый вопрос: нет. Он находит кратные путем подсчета, как и должно быть в любом истинном решете _{(хотя «минус» в списках, конечно, неэффективен; вышеприведенное улучшает это, перестраивая цепочку линейного вычитания ((((xs-a)-b)-c)- ... ) в вычитание древовидных дополнений, xs-(a+((b+c)+...))).}

Как уже было сказано, использование изменяемых массивов обеспечивает наилучшую производительность. Следующий код получен из этой версии "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

EDIT: Приведенный выше код был исправлен, а приведенное ниже время отредактировано, чтобы отразить исправление, а также были отредактированы комментарии, сравнивающие выгружаемую версию с невыгружаемой.

Время выполнения этого до указанных верхних диапазонов указано в таблице комментариев внизу приведенного выше исходного кода, измеренного ideone. .com и примерно в пять раз быстрее, чем ответ, опубликованный @WillNess, также измеренный на ideone.com. Этот код требует тривиального времени для отбраковки простых чисел до двух миллионов и всего 1,26 секунды для отбраковки до ста миллионов. Это время примерно в 2,86 раза быстрее при работе на моем i7 (3,5 ГГц) за 0,44 секунды до ста миллионов и за 6,81 секунды до одного миллиарда. Использование памяти составляет чуть более шести мегабайт для первого и шестидесяти мегабайт для второго, что является памятью, используемой огромным (битовым) массивом. Этот массив также объясняет нелинейную производительность тем, что по мере того, как размер массива превышает размер кэша ЦП, среднее время доступа к памяти ухудшается из-за отбраковки представления составного числа.

-- 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 ())))

На ideone.com этот код вычисляет количество простых чисел до ста миллионов примерно за 1,13 секунды и примерно за 12 секунд до единицы. миллиард (32-битная настройка). Вероятно, колесная факторизация и определенно многоядерная обработка сделают ее еще быстрее на многоядерном процессоре. На моем i7 (3,5 ГГц) просеивание до ста миллионов занимает 0,44 секунды, а до миллиарда — 4,7 секунды, при примерно линейной производительности с увеличением диапазона, как и ожидалось. Похоже, что в версии GHC, запущенной на ideone.com, есть какие-то нелинейные накладные расходы, которые имеют некоторые потери производительности для больших диапазонов, которых нет в i7, и которые, возможно, связаны с другой сборкой мусора, поскольку страница буферы создаются новые для каждой новой страницы. END_EDIT_ADD

EDIT_ADD2: кажется, что большая часть времени обработки приведенного выше сегментированного кода страницы используется в (ленивой) обработке списка, поэтому код соответственно переформулирован с несколькими улучшениями следующим образом:

1. Реализована функция подсчета простых чисел, которая не использует обработку списка и использует поиск в таблице «popCount» для подсчета количества «единичных» битов в 32-битном слове за раз. Таким образом, время на получение результатов незначительно по сравнению с фактическим временем отбраковки через сито.

2. Сохранение базовых простых чисел в виде списка битовых сегментов страницы, что намного более эффективно, чем хранение списка простых чисел, а время на преобразование сегментов страницы в простые числа по мере необходимости не требует больших вычислительных затрат.

3. Функция создания простого сегмента настроена таким образом, что для начального сегмента нулевой страницы он использует собственный битовый шаблон в качестве исходной страницы, что делает код отбраковки составных чисел короче и проще.

Затем код становится следующим:

{-# 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 с 0,02, 0,03, 0,05. , 0,30, 3,0 и 9,1 секунды, необходимые для работы на моем i7-2700K (3,5 ГГц) на один, два, десять, сто, тысячу (один миллиард) и три тысячи (три миллиарда) миллионов соответственно, с в значительной степени постоянный объем памяти, медленно увеличивающийся с количеством базовых простых чисел, меньшего, чем квадратный корень из диапазона, как требуется. При компиляции с помощью серверной части компилятора LLVM это время становится равным 0,01, 0,02, 0,02, 0,12, 1,35 и 4,15 секунды соответственно из-за более эффективного использования регистров и машинных инструкций; последнее довольно близко к той же скорости, что и при компиляции с помощью 64-битного компилятора, а не 32-битного компилятора, используемого, поскольку эффективное использование регистров означает, что наличие дополнительных регистров не имеет большого значения.

Как указано в коде, соотношение между производительностью на моей реальной машине и серверами ideone.com становится намного меньше, чем для гораздо более расточительных алгоритмов памяти, из-за того, что они не ограничиваются узкими местами доступа к памяти, поэтому ограничение скорости в основном просто соотношение Тактовые частоты ЦП, а также эффективность обработки ЦП за такт. Однако, как отмечается там, существует какая-то странная неэффективность генератора собственного кода GHC (NCG) при запуске под Windows (32-разрядный компилятор), поскольку время выполнения более чем на 50% медленнее, чем при запуске под Linux (как ideone. ком-сервер использует). Насколько я знаю, они оба имеют общую кодовую базу для одной и той же версии Haskell GHC, и единственное расхождение заключается в используемом компоновщике (который также используется с бэкэндом LLVM, который не затрагивается), поскольку GHC NCG не использует GCC, а только ассемблер mingw32. , что также должно быть одинаковым.

Обратите внимание, что этот код, скомпилированный с помощью серверной части компилятора LLVM, имеет примерно ту же скорость, что и тот же алгоритм, написанный для высокооптимизированных реализаций «C/C++», что указывает на то, что Haskell действительно имеет возможность разрабатывать код с очень плотным циклом. Можно сказать, что код на Haskell намного читабельнее и безопаснее, чем эквивалентный код «C/C++», если привыкнуть к парадигмам монадического и нестрогого кода на Haskell. Дальнейшие улучшения скорости выполнения для решета Эратосфена являются исключительно функцией настройки используемых реализаций, а не выбором языка между Haskell и C/C++.

Вывод: Конечно, это еще не максимальная скорость для Haskell-версии Решета Эратосфена, поскольку мы до сих пор не настроили доступ к памяти, чтобы более эффективно использовать быстрый кэш L1 процессора. , мы также не уменьшили значительно общее количество составных операций отбраковки, необходимых с использованием экстремальной факторизации колеса, за исключением исключения обработки шансов. Однако этого достаточно, чтобы ответить на вопрос, показывающий, что изменяемые массивы являются наиболее эффективным способом решения таких проблем с узким циклом, с потенциальным выигрышем в скорости примерно в 100 раз по сравнению с использованием списков или неизменяемых массивов. END_EDIT_ADD2