LINUX.ORG.RU

Как быстрее найти сумму простых чисел, меньших 2e6, на хаскелле?

 , ,


0

3

На форуме projecteuler наибыстрейшие реализации на других языках по сути решето эратосфена используют, и предлагаемый там(хотя может на других страницах есть другие) код на хаскелле тоже его использует

module Main where

seive :: [Integer] -> [Integer] -> [Integer]
seive ps []     = reverse ps
seive ps (p:cs) = 
    seive (p:ps) (filter ((/= 0) . (`mod` p)) cs)

main = print . sum $ seive [] [2..1999999]
.

Я пока пришёл к выводу, что вариант

isPrime :: Integral a => a -> Bool
check n factor = if factor <= floor(sqrt(fromIntegral n)) then if n `mod` factor == 0 then False;
	else check n (factor+2);
	else True;
isPrime x | x == 2 = True
					| x `mod` 2 == 0 = False
					| otherwise = check x 3
primeSum i n sum = if i >= n then sum;
	else if isPrime(i) then primeSum(i+1) n (sum+i);
		else primeSum (i+1) n sum
, где primeSum 2 1999999 0 даёт результат(должно), быстрее, по крайней мере в ghci(для 64000). На С практически такой же такой же алгоритм(цикл вместо рекурсии, конечно) секунд 10 считает на i5 2.67 GHz.

Может дело просто в том. что не в ghci нужно запускать, и тогда хотя бы секунд за 30 посчитает?

sudo cast dave, Macil

★★★★

Последнее исправление: aptyp (всего исправлений: 1)

Не совсем понял, в чем вопрос. Найти более быстрый алгоритм, или как использовать компилятор ghc вместо интерпретатора ghci?

Если второе, то можно собрать программу так (предполагаем, что основной файл с main называется Main.hs):

ghc --make -O2 Main.hs

Кстати, компиляция, профилирование и прочие технические вещи хорошо описаны в общедоступной книге «Real World Haskell».

Если все же вопрос касается оптимального алгоритма, то тут мне что-то совсем голову не хочется нагружать :) У меня еще есть рабочие задачи на сегодня.

dave ★★★★★
()
Ответ на: комментарий от aptyp

А понятно. Могу сказать, что задача не так проста. По-моему верхний код очень сильно напоминает то, что написано в SICP в главе про потоки. В хаскеле список можно трактовать как бесконечный поток, но делать это надо умеючи. На мой взгляд сложная тема, которая только кажется простой. Так вот, через поток и реализуется это решето. У тебя же алгоритм по-моему работает в-лоб с небольшой оптимизацией четных делителей.

dave ★★★★★
()
Ответ на: комментарий от dave

Такой поток по-английски называется stream.

dave ★★★★★
()

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

KblCb ★★★★★
()
Ответ на: комментарий от dave

Да нет) Просто захотелось почитать насчёт lazy evaluation, бесконечные потоки восхитили, хотя и не эффективны.
Дальше по треду на сайте projecteuler были примеры, когда уже найденные простые числа сохранялись, вроде ещё быстрее.

aptyp ★★★★
() автор топика
Ответ на: комментарий от aptyp

бесконечные потоки восхитили, хотя и не эффективны.

нужно просто уметь ими пользоваться, и собирать с оптимизациями, в этом случае список не конструируется, а разбирается inplace.

qnikst ★★★★★
()
Ответ на: комментарий от qnikst

Я сознательно пока к монадам не хочу лезть.

aptyp ★★★★
() автор топика

Во-первых «factor <= floor(sqrt(fromIntegral n))» замени на «factor*factor < n» Умножение быстрее.

Во-вторых тут слишком много делений. Правильнее вычёркивать, но тут засада. Так что думай.

AlexVR ★★★★★
()
Ответ на: комментарий от AlexVR

ммм, ну да. Люди и вычёркивали:-) Буду думать.

aptyp ★★★★
() автор топика

Для затравки могу дать ещё такое плохое решение:

primes :: [Integer]
primes = 2: 3: sieve 0 (tail primes) 3
sieve k (p:ps) x = [n | n <- [x+2,x+4..p*p-2], and [n`rem`p/=0 | p <- fs]]
    ++ sieve (k+1) ps (p*p) 
    where fs = take k (tail primes)
AlexVR ★★★★★
()
Ответ на: комментарий от aptyp

Если мне память не изменяет, решето Эратосфена требует изменяемых структур данных, которых за пределами ST (и C-шной кучи) просто нет.

KblCb ★★★★★
()
Ответ на: комментарий от KblCb

А вот это что такое?

seive :: [Integer] -> [Integer] -> [Integer]
seive ps []     = reverse ps
seive ps (p:cs) = 
    seive (p:ps) (filter ((/= 0) . (`mod` p)) cs)
Или я ошибаюсь?

aptyp ★★★★
() автор топика
Ответ на: комментарий от aptyp

на immutable данных же. считай у тебя тут создаётся куча списков, в итоге это работает, но не эффективно.

qnikst ★★★★★
()
Ответ на: комментарий от aptyp

чтобы в вычислениях не использовались предыдущие элементы списка, но это не для этой задачи. Напр. mapM_ print [ very_difficult_function x| x<-[1..1e100]] не будет создавать списка.

qnikst ★★★★★
()
Ответ на: комментарий от unC0Rr

иммутабельный массив, использующийся не только ради поиска по индексу это вообще ужас.

qnikst ★★★★★
()

Довольно скучная задача. Более интересна задача: найти, например, десятимиллиардное простое число.

buddhist ★★★★★
()
Ответ на: комментарий от buddhist

Там была задача 10001 найти.

aptyp ★★★★
() автор топика
Ответ на: комментарий от buddhist

Она скучная, но долгая, я минуты три ждал ответа.

aptyp ★★★★
() автор топика

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

buddhist ★★★★★
()
Ответ на: комментарий от buddhist

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

qnikst ★★★★★
()
Ответ на: комментарий от buddhist

Ты многого требуешь. Человек недавно приступил к изучению языка. И на мой взгляд такую задачу еще рано рассматривать.

dave ★★★★★
()
Ответ на: комментарий от aptyp

ну менять алгоритм на более адекватный, учитывающий эту особенность. Ну или использовать ST, который чистый снаружи, главное не обращать внимание на слово «монада», т.к. ничего страшного в этом нет.

qnikst ★★★★★
()
Ответ на: комментарий от dave

Человек недавно приступил к изучению языка. И на мой взгляд такую задачу еще рано рассматривать.

Как раз самое то, решение небольших задач с http://projecteuler.net/, http://spoj.pl и т.д. очень хорошо помогают разобраться в новом языке, понять его особенности достоинства и недостатки, где-то поставить задачи для дальнейшего изучения.

AlexVR ★★★★★
()
Ответ на: комментарий от AlexVR

Да уж))) День уйдёт разобраться.

aptyp ★★★★
() автор топика
Ответ на: комментарий от AlexVR

Да, порядка 10 задачек решил, просто вот тут наткнулся на ограничения, ушёл от решета, и за 2.5 минуты посчиталось.
правда дальше хуже, с треугольными числами) стека не хватает))

aptyp ★★★★
() автор топика
Ответ на: комментарий от aptyp

это не ограничения, ты точно делаешь, что-то не то

qnikst ★★★★★
()
Ответ на: комментарий от aptyp

Что-то долго на моём целероне это меньше занимает

AlexVR ★★★★★
()
{-# LANGUAGE BangPatterns #-}

import Data.Bits
import Data.Array.Base
import Control.Applicative
import Control.Monad
import Control.Monad.ST
import System.Environment

newSieveArray :: Int -> Int -> ST s (STUArray s Int Bool)
newSieveArray n m = newArray (n, m) True

sumPrimes :: Int -> Int
sumPrimes m = runST $ newSieveArray 3 m >>= go 3 0 where
  max :: Int
  max = 1 + truncate (sqrt $ fromIntegral m)
  go :: Int -> Int -> STUArray s Int Bool -> ST s Int
  go !n !c !a = if n >= m then return c else do
    e <- unsafeRead a n
    if e then
      if n < max then let
          loop !j = if j >= m then go (n + 2) (c + n) a else do
            x <- unsafeRead a j
            when x $ unsafeWrite a j False
            loop (j + n)
          in loop (n ^ 2)
        else go (n + 2) (c + n) a
      else go (n + 2) c a

main :: IO ()
main = print . (2 +) . sumPrimes . (1 +) =<< read . head <$> getArgs
$ ghc -O3 -o p p.hs
[1 of 1] Compiling Main             ( p.hs, p.o )
Linking p ...
$ time ./p 1000000   
37550402023
./p 1000000  0.01s user 0.00s system 95% cpu 0.017 total
$ time ./p 2000000
142913828922
./p 2000000  0.03s user 0.00s system 89% cpu 0.035 total
$ time ./p 1000000000
24739512092254535
./p 1000000000  24.66s user 0.10s system 99% cpu 24.767 total
$ time ./p 2000000000
95673602693282040
./p 2000000000  51.45s user 0.21s system 99% cpu 51.701 total

Тут IOUArray никаких преимуществ не даст, тогда как STUArray позволяет вытащить результат из ST (runST). На 32 битах может иметь смысл взять Integer вместо Int для параметра c.

quasimoto ★★★★
()
Ответ на: комментарий от quasimoto

[code][denis@delete[nixos]:~]$ ghc -O3 -o p p.hs
[1 of 1] Compiling Main ( p.hs, p.o )
Linking p ...

[denis@delete[nixos]:~]$ time ./p 1000000
-1104303641

real 0m0.030s
user 0m0.019s
sys 0m0.004s

[denis@delete[nixos]:~]$ time ./p 1000000
-1104303641

real 0m0.029s
user 0m0.018s
sys 0m0.004s

[denis@delete[nixos]:~]$ time ./p 2000000
1179908154

real 0m0.048s
user 0m0.036s
sys 0m0.002s

[denis@delete[nixos]:~]$ time ./p 1000000000
-2043879097

real 0m47.153s
user 0m38.164s
sys 0m0.168s

[denis@delete[nixos]:~]$
[/code]

Простите, но что за цифры выплевывает ваша программа и какое они имеют отношение к задаче ТС?

delete83 ★★
()
Ответ на: комментарий от delete83

Простите, но что за цифры выплевывает ваша программа и какое они имеют отношение к задаче ТС?

У меня 64-битная система, поэтому можно обходиться арифметикой фиксированной точности:

maxBound :: Int = 9223372036854775807 = 2 ^ 63 - 1

на 32-битной системе, как я сказал, Int может быть недостаточно:

maxBound :: Int = 2147483647 = 2 ^ 31 - 1

поэтому нужно использовать длинную арифметику Integer:

sumPrimes :: Int -> Integer
sumPrimes m = runST $ newSieveArray 3 m >>= go 3 0 where
  max :: Int
  max = 1 + truncate (sqrt $ fromIntegral m)
  go :: Int -> Integer -> STUArray s Int Bool -> ST s Integer
  go !n !c !a = if n >= m then return c else do
    e <- unsafeRead a n
    if e then
      if n < max then let
          loop !j = if j >= m then go (n + 2) (c + fromIntegral n) a else do
            x <- unsafeRead a j
            when x $ unsafeWrite a j False
            loop (j + n)
          in loop (n ^ 2)
        else go (n + 2) (c + fromIntegral n) a
      else go (n + 2) c a

При переходе Int -> Integer для чисел < maxBound :: Int разница в скорости мало заметна, то есть на 64-битной системе такой вариант будет работать примерно также, на 32-битной - будут задействованы функции gmp.

Ну и ещё на 64-битной системе имеет смысл вызывать:

./p 9223372036854775807

на 32-битной - только

./p 2147483647

иначе нужно переходить на STUArray s Integer Bool.

Ещё нужно заметить, что newSieveArray аллоцирует массив в (m - n) / 8 байт (./p 1000000000 аллоцирует (1000000000 / 8) / 1024 / 1024 ~= 120MB, у меня в task manager даже точно такое значение и показывается), а не (m - n) * sizeof(int) (для ./p 1000000000 это было бы (1000000000 * 8) / 1024 / 1024 / 1024 ~= 7GB). То есть массив Bool-ов это битовый массив. Ну и этот достаточно большой массив существует в единственном экземпляре в памяти, рекурсивные функции работают как циклы, изменяющие этот массив (тут строгость с BangPatterns помогает, иначе бы стек забился).

quasimoto ★★★★
()
Ответ на: комментарий от quasimoto

на 32-битной - будут задействованы функции gmp

Проверил в chroot-е, работает:

$ time ./p 1000000000    
24739512092254535
./p 1000000000  34,13s user 0,14s system 99% cpu 34,280 total
$ time ./p 2000000000    
95673602693282040
./p 2000000000  68,66s user 0,26s system 99% cpu 1:08,94 total

всего на 10-20% медленнее.

quasimoto ★★★★
()
Ответ на: комментарий от quasimoto

./p 9223372036854775807

Только это миллиарды гигабайт :)

./p 2147483647

Это 1/2 гигабайта.

quasimoto ★★★★
()
Ответ на: комментарий от aptyp

Это, фактически, адаптация сишного императивного кода - мутабельный массив в качестве плоской памяти, хвостово-рекурсивные функции в качестве циклов и их аккумуляторы в качестве переменных. Вот ещё немного короче и быстрее (без sqrt, как выше сказали):

sumPrimes :: Int -> Int
sumPrimes m = runST $ newArray (3, m) True >>= go 3 0 where
  go :: Int -> Int -> STUArray s Int Bool -> ST s Int
  go !n !c !a = if n >= m then return c else do
    e <- unsafeRead a n
    if e then
      if n * n < m then let
        loop !j = if j >= m then go (n + 2) (c + n) a else do
          x <- unsafeRead a j
          when x $ unsafeWrite a j False
          loop (j + n)
        in loop (n ^ 2)
        else go (n + 2) (c + n) a
      else go (n + 2) c a
quasimoto ★★★★
()
Вы не можете добавлять комментарии в эту тему. Тема перемещена в архив.