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

e-maxx algorithm оригинал: C/C++ #algorithm #emaxx #math #number-theory

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

Дано число

. Требуется найти все простые в отрезке

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

за время

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

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

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

решето Эратосфена: требуется заводить массив из

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

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

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

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

эти два алгоритма.)

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

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

от в отрезке

его минимальный простой делитель

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

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

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

Будем теперь перебирать текущее число

от

до

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

— это означает, что число

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

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

и добавить

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

.

— это означает, что текущее число

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

Реализация

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

.

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

}

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

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

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

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

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

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

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

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

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

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

.

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

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

где

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

, а число

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

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

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

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

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

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

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

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

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

Хотя асимптотика

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

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

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

длины

и массив всех простых

длины примерно

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

Однако спасает его то, что массив

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

в отрезке

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

Литература

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

C# решение

auto-draft, проверить перед отправкой
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++ решение

сопоставлено/оригинал
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 решение

auto-draft, проверить перед отправкой
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];
    }
}

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

Вакансии для этой задачи

Показаны активные вакансии с пересечением по тегам задачи.

Все вакансии
Активных вакансий пока нет.