Эффективное сито Эйлера в стиле потоковой обработки

Question

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

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

поэтому я пытаюсь написать самостоятельно:

{-sieve of Euler-}
primesEU = 2:sieve [3,5 ..] where
    sieve (p:i:xt) = p:(sieve (i:xt) `minus` (lps i)) where
        lps i = map (i*) (takeWhile' (\p->(p<i)&&(mod i p /= 0)) primesEU)

takeWhile'::(t->Bool)->[t]->[t]
takeWhile' p [] = []
takeWhile' p (x:xs) = case (p x) of
    True  -> x:takeWhile' p xs
    False -> [x]

minus::(Ord t) => [t]->[t]->[t]
minus [] _ = []
minus xs [] = xs 
minus xs@(x:xt) ys@(y:yt) 
                      |x<y  = x:(minus xt ys)
                      |x==y = minus xt yt
                      |x>y  = minus xs yt

minus похож на minus в Data.List.Ord,

takeWhile' похож на takeWhile, но имеет небольшую разницу: takeWhile отбрасывает первый элемент, который не удовлетворяет предикату; takeWhile' возьму это.

lsp i возвращает конечный поток произведений i и простых чисел не больше, чем наименьший простой коэффициент i.

К сожалению, моя реализация работает невероятно медленно... и я не могу найти виновника.

В любом случае, возможно ли реализовать эффективное сито Эйлера в стиле потоковой обработки? Или алгоритм имеет что-то присущее против природы потока?

Answer 1

primesEU = 2:sieve [3,5 ..] where
    sieve (p:i:xt) = p:(sieve (i:xt) `minus` (lps i)) where
        lps i = map (i*) (takeWhile' (\p->(p<i)&&(mod i p /= 0)) primesEU)

Прежде всего, реализация некорректна, она считает 9 простых. Я подозреваю, что вы имели в виду sieve (i:xt) `minus` lps p там.

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

Так что здесь происходит?

Суть этого (я еду на ужин через мгновение, когда я вернусь, расширится) в том, что каждый премьер p должен пройти через все minus звонки, созданные нечетными номерами (>= 3) меньше чем p Поэтому не совсем, несколько minus устранены до]. Это (p-3)/2 слои, каждый стоит один шаг.

Таким образом, постоянные факторы и композиты в стороне, производя простые числа <= n расходы O(∑ p)где сумма над простыми числами не превышает n, Это O(n²/log n),

Давайте проследим несколько шагов оценки sieve [3,5 .. ], Для краткости обозначу через s k список sieve [k, k+2 ..] а также minus от (\\), Я использую исправленное определение

sieve (k:ks) = k : (sieve ks `minus` lps k)

это не считает 9 премьер. Я сразу же расширяю список магазинов, что несколько сокращает объем работы. Мы получаем

 (:)
 / \
3  (\\)
   /  \
 s 5  [9]

сразу, а затем 3 могут быть использованы непосредственно. После этого minus вершина правильного поддерева требует, чтобы s 5 оцениваться достаточно, чтобы сказать, пусто ли оно.

   (\\)
   /  \
 (:)  [9]
 / \
5  (\\)
   /  \
 s 7  [15,25]

Это делается за один шаг. (Затем необходимо определить, lps 3 пусто, и аналогично позже, но для каждого члена lps k, это только один шаг, поэтому мы можем игнорировать это.) Тогда minus нужно сравнение, чтобы пузырь 5 на вершину, оставляя

   (:)
   / \
  5  (\\)
     /  \
  (\\)  [9]
  /  \
s 7  [15,25]

Теперь, после того, как 5 был потреблен, правый ребенок сверху (:) должен быть расширен

      (\\)
      /  \
   (\\)  [9]
   /  \
 (:)  [15,25]
 / \
7  (\\)
   /  \
 s 9  [21,35,49]

Одно сравнение, чтобы переместить 7 мимо кратных 5

     (\\)
     /  \
   (:)  [9]
   / \
  7  (\\)
     /  \
  (\\)  [15,25]
  /  \
s 9  [21,35,49]

и еще один, чтобы поднять его за кратное 3

      (:)
      / \
     7  (\\)
        /  \
     (\\)  [9]
     /  \
  (\\)  [15,25]
  /  \
s 9  [21,35,49]

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

         (\\)
         /  \
      (\\)  [9]
      /  \
   (\\)  [15,25]
   /  \
 (:)  [21,35,49]
 / \
9  (\\)
   /  \
 s 11 [27]

Этот кандидат должен быть поднят за два minusесли наконец достичь minus это устраняет

        (\\)
        /  \
     (\\)  [9]
     /  \
   (:)  [15,25]
   / \
  9  (\\)
     /  \
  (\\)  [21,35,49]
  /  \
s 11 [27]

        (\\)
        /  \
      (:)  [9]
      / \
     9  (\\)
        /  \
     (\\)  [15,25]
     /  \
  (\\)  [21,35,49]
  /  \
s 11 [27]

Теперь minus может сделать свою работу впервые, производя

           (\\)
           /  \
        (\\)  []
        /  \
     (\\)  [15,25]
     /  \
  (\\)  [21,35,49]
  /  \
s 11 [27]

(xs `minus` [] вверху удаляется пустой список только после того, как первый аргумент был определен непустым; поменять местами первые два уравнения minus изменит это, но разница в необходимых шагах будет небольшой). Затем оценка должна перейти на четыре уровня в дереве, чтобы получить следующего кандидата (11). Этот кандидат должен быть поднят за четыре minus пока он не достигнет вершины. Один minus при этом исключается из дерева, поэтому следующий кандидат (13) создается только на четыре уровня вниз по дереву, а не на пять ((13-3)/2), таким образом, количество шагов, чтобы принести p наверх, так что его можно употреблять не совсем (p-3)/2, но это не намного меньше.

Последний элемент в lps k всегда делится на квадрат наименьшего простого множителя kи все они странные, так что их максимум

1/2*(n/3² + n/5² + n/7² + n/11² + ...) ~ n/10

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

Таким образом, проблема в том, что каждое простое число производится глубоко в дереве, по крайней мере, p*0.4 уровни сверху. Это означает подъем p на вершину, чтобы сделать его потребляемым занимает по крайней мере p*0.4 шаги, дающие в основном квадратичную общую сложность, не считая работы, необходимой для составления списков кратных.

Тем не менее, фактическое поведение вначале хуже, оно приближается к квадратичному поведению, когда вычисление простых чисел в области между 100000 и 200000 - должно стать квадратичным по модулю логарифмическим фактором в пределе, но я не думаю, что у меня хватит терпения ждать миллион или два.

В любом случае, возможно ли реализовать эффективное сито Эйлера в стиле потоковой обработки? Или алгоритм имеет что-то присущее против природы потока?

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

Что делает Сито Эратосфена эффективным, так это то, что это тривиально, чтобы достичь следующего кратного простого числа (просто добавьте p или же 2*p или так к индексу), не заботясь о том, было ли оно уже вычеркнуто как кратное меньшего простого числа. Работа, необходимая для избежания многократного пересечения, намного больше, чем работа многократного пересечения.