MAXimal

добавлено: 6 Sep 2011 1:03
редактировано: 23 Mar 2012 3:58

Решето Эратосфена с линейным временем работы

Дано число n. Требуется найти все простые в отрезке [2; n].

Классический способ решения этой задачи — решето Эратосфена. Этот алгоритм очень прост, но работает за время O (n \log \log n).

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

Кроме того, приводимый здесь алгоритм в качестве "побочного эффекта" фактически вычисляет факторизацию всех чисел в отрезке [2; n], что может быть полезно во многих практических применениях.

Недостатком приводимого алгоритма является то, что он использует больше памяти, чем классическое решето Эратосфена: требуется заводить массив из n чисел, в то время как классическому решету Эратосфена достаточно лишь n бит памяти (что получается в 32 раза меньше).

Таким образом, описываемый алгоритм имеет смысл применять только до чисел порядка 10^7, не более.

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

Описание алгоритма

Наша цель — посчитать для каждого числа i от в отрезке [2; n] его минимальный простой делитель lp[i].

Кроме того, нам потребуется хранить список всех найденных простых чисел — назовём его массивом pr[].

Изначально все величины lp[i] заполним нулями, что означает, что мы пока предполагаем все числа простыми. В ходе работы алгоритма этот массив будет постепенно заполняться.

Будем теперь перебирать текущее число i от 2 до n. У нас может быть два случая:

  • lp[i] = 0 — это означает, что число i — простое, т.к. для него так и не обнаружилось других делителей.

    Следовательно, надо присвоить lp[i] = i и добавить i в конец списка pr[].

  • lp[i] \ne 0 — это означает, что текущее число i — составное, и его минимальным простым делителем является lp[i].

В обоих случаях дальше начинается процесс расстановки значений в массиве lp[]: мы будем брать числа, кратные i, и обновлять у них значение lp[]. Однако наша цель — научиться делать это таким образом, чтобы в итоге у каждого числа значение lp[] было бы установлено не более одного раза.

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

 x_j = i \cdot p_j,

где последовательность p_j — это все простые, не превосходящие lp[i] (как раз для этого нам понадобилось хранить список всех простых чисел).

У всех чисел такого вида проставим новое значение lp[x_j] — очевидно, оно будет равно p_j.

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

Реализация

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

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];
}

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

Доказательство корректности

Докажем корректность алгоритма, т.е. что он корректно расставляет все значения lp[], причём каждое из них будет установлено ровно один раз. Отсюда будет следовать, что алгоритм работает за линейное время — поскольку все остальные действия алгоритма, очевидно, работают за O (n).

Для этого заметим, что у любого числа i единственно представление такого вида:

 i = lp[i] \cdot x,

где lp[i] — (как и раньше) минимальный простой делитель числа i, а число x не имеет делителей, меньших lp[i], т.е.:

 lp[i] \le lp[x].

Теперь сравним это с тем, что делает наш алгоритм — он фактически для каждого x перебирает все простые, на которые его можно домножить, т.е. простые до lp[x] включительно, чтобы получить числа в указанном выше представлении.

Следовательно, алгоритм действительно пройдёт по каждому составному числу ровно один раз, поставив у него правильное значение lp[].

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

Время работы и требуемая память

Хотя асимптотика O (n) лучше асимптотики O (n \log \log n) классического решета Эратосфена, разница между ними невелика. На практике это означает лишь двукратную разницу в скорости, а оптимизированные варианты решета Эратосфена и вовсе не проигрывают приведённому здесь алгоритму.

Учитывая затраты памяти, которые требует этот алгоритм — массив чисел lp[] длины n и массив всех простых pr[] длины примерно n / \ln n — этот алгоритм кажется уступающим классическому решету по всем статьям.

Однако спасает его то, что массив lp[], вычисляемый этим алгоритмом, позволяет искать факторизацию любого числа в отрезке [2; n] за время порядка размера этой факторизации. Более того, ценой ещё одного дополнительного массива можно сделать, чтобы в этой факторизации не требовались операции деления.

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

Литература

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