E020. Sieve of Eratosthenes с линейным временем работы

e-maxx algorithm original: C/C++ #algorithm #emaxx #math #number-theory
Task text is translated from Russian for the selected interface language. Code is left unchanged.

Источник: e-maxx.ru/algo, страница PDF 54.

given number

. it is required find все простые в отрезке

. Классический способ решения этой задачи — Sieve of Eratosthenes. Этот Algorithm очень прост, но работает

за время

. Хотя в настоящий момент известно достаточно много Algorithmов, работающих за сублинейное время (т.е. за ), описываемый ниже Algorithm интересен своей простотой — он практически не сложнее классического решета Эратосфена. Кроме того, приводимый здесь Algorithm в качестве "побочного эффекта" фактически вычисляет

факторизацию всех чисел в отрезке

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

Sieve of Eratosthenes: it is required заводить array из

чисел, в то время как классическому решету Эратосфена

достаточно лишь

бит памяти (что получается в

раза меньше). Таким образом, описываемый Algorithm имеет смысл применять только до чисел порядка , не более. Авторство Algorithmа, по всей видимости, принадлежит Грайсу и Мисра (Gries, Misra, 1978 г. — см. список литературы в конце). (И, собственно говоря, называть данный Algorithm "решетом Эратосфена" некорректно: слишком отличаются

эти два Algorithmа.)

Описание Algorithmа

Наша цель — посчитать для каждого числа

от в отрезке

его minimum простой делитель

. Кроме того, нам поit is required хранить список всех найденных простых чисел — назовём его arrayом .

Изначально все величины

заполним нулями, что означает, что мы пока предполагаем все числа простыми. В ходе работы Algorithmа этот array будет постепенно заполняться.

Будем теперь перебирать текущее number

от

до

. У нас может быть два случая:

— это означает, что number

— простое, т.к. для него так и не обнаружилось других делителей.

Следовательно, надо присвоить

и добавить

в конец списка

.

— это означает, что текущее number

— составное, и его минимальным простым делителем является

. В обоих случаях дальше начинается процесс расстановки значений в arrayе

: мы будем брать

числа, кратные

, и обновлять у них значение

. Однако наша цель — научиться делать это таким образом, чтобы

в итоге у каждого числа значение

было бы установлено не более одного раза. Утверждается, что для этого можно поступить таким образом. Рассмотрим числа вида:

где последовательность

— это все простые, не превосходящие

(как раз для этого нам понадобилось

хранить список всех простых чисел).

У всех чисел такого вида проставим новое значение

— очевидно, оно будет равно

. Почему такой Algorithm корректен, и почему он работает за линейное время — см. ниже, пока же приведём его реализацию.

Implementation

Решето выполняется до указанного в константе числа

.

const int N = 10000000;

int lp[N+1];

vector<int> pr;

for (int i=2; i<=N; ++i) {
if (lp[i] == 0) {

lp[i] = i;

pr.push_back (i);

}

for (int j=0; j<(int)pr.size() && pr[j]<=lp[i] && i*pr[j]<=N; ++j)

lp[i * pr[j]] = pr[j];

}

Эту реализацию можно немного ускорить, избавившись от вектора

(заменив его на обычный array со счётчиком),

а также избавившись от дублирующегося умножения во вложенном цикле

(для чего результат произведения

надо просто запомнить в какой-либо переменной).

Correctness proof

Докажем корректность Algorithmа, т.е. что он корректно расставляет все значения

, причём каждое из них

будет установлено ровно один раз. Отсюда будет следовать, что Algorithm работает за линейное время — поскольку

все остальные действия Algorithmа, очевидно, работают за

.

Для этого заметим, что у любого числа

единственно представление такого вида:

где

— (как и раньше) minimum простой делитель числа

, а number

не имеет делителей, меньших

, т.е.: Теперь сравним это с тем, что делает наш Algorithm — он фактически для каждого

перебирает все простые, на

которые его можно домножить, т.е. простые до

включительно, чтобы получить числа в указанном

выше представлении. Следовательно, Algorithm действительно пройдёт по каждому составному числу ровно один раз, поставив у

него правильное значение

. Это означает корректность Algorithmа и то, что он работает за линейное время.

Running time и требуемая память

Хотя Asymptotic complexity

лучше асимптотики

классического решета Эратосфена, разница между

ними невелика. На практике это означает лишь двукратную разницу в скорости, а оптимизированные варианты решета Эратосфена и вовсе не проигрывают приведённому здесь Algorithmу. given затраты памяти, которые требует этот Algorithm — array чисел

длины

и array всех простых

длины Exampleно

— этот Algorithm кажется уступающим классическому решету по всем статьям.

Однако спасает его то, что array

, вычисляемый этим Algorithmом, позволяет искать факторизацию любого числа

в отрезке

за время порядка размера этой факторизации. Более того, ценой ещё одного дополнительного arrayа можно сделать, чтобы в этой факторизации не требовались операции деления. Знание факторизации всех чисел — очень полезная информация для некоторых задач, и этот Algorithm является одним из немногих, которые позволяют искать её за линейное время.

References

● David Gries, Jayadev Misra. A Linear Sieve Algorithm for Finding Prime Numbers [1978]

C# solution

auto-draft, review before submit
using System;
using System.Collections.Generic;
using System.Linq;

public static class AlgorithmDraft
{
    // Auto-generated C# draft from the original e-maxx C/C++ listing. Review before production use.
    const int N = 10000000;
    int lp[N+1];
    List<int> pr;
    for (int i=2; i<=N; ++i) {
            if (lp[i] == 0) {
                    lp[i] = i;
                    pr.push_back (i);
            }
            for (int j=0; j<(int)pr.size() && pr[j]<=lp[i] && i*pr[j]<=N; ++j)
                    lp[i * pr[j]] = pr[j];
    }
}

C++ solution

matched/original
const int N = 10000000;
int lp[N+1];
vector<int> pr;
for (int i=2; i<=N; ++i) {
        if (lp[i] == 0) {
                lp[i] = i;
                pr.push_back (i);
        }
        for (int j=0; j<(int)pr.size() && pr[j]<=lp[i] && i*pr[j]<=N; ++j)
                lp[i * pr[j]] = pr[j];
}

Java solution

auto-draft, review before submit
import java.util.*;
import java.math.*;

public class AlgorithmDraft {
    // Auto-generated Java draft from the original e-maxx C/C++ listing. Review before production use.
    const int N = 10000000;
    int lp[N+1];
    ArrayList<Integer> pr;
    for (int i=2; i<=N; ++i) {
            if (lp[i] == 0) {
                    lp[i] = i;
                    pr.push_back (i);
            }
            for (int j=0; j<(int)pr.size() && pr[j]<=lp[i] && i*pr[j]<=N; ++j)
                    lp[i * pr[j]] = pr[j];
    }
}

Материал разбит как Algorithmическая Task: изучить постановку, понять асимптотику и реализовать Algorithm на выбранном языке.

Vacancies for this task

Active vacancies with overlapping task tags are shown.

All vacancies
There are no active vacancies yet.