Решето Эратосфена 

Статья написана для Сборника Олпрогера.

Достаточно часто в задачах бывает полезно - уметь факторизировать числа. В большинстве таких задач как раз используется решето Эратосфена.

Вспомним несколько определений.

Простое число - число, которое имеет ровно два различных натуральных делителя — единицу и самого себя.

Если число не простое, то оно составное.

Единица не является ни простым, ни составным числом.

Давайте сначала решим задачу капельку полегче, чем исходная: для каждого числа от 1 до n поймем является ли оно простым или нет.

Список простых чисел: 2, 3, 5, 7, 11, 13, 17, ...

Заметим, что если число x делится на какое-либо из чисел от 2 до x - 1, то оно составное, иначе простое. Сделаем какие-то замечания:

• Все числа, которые делятся на 2, кроме самой двойки точно не простые: 4, 6, 8, 10, 12, ...
• Все числа, которые делятся на 3, кроме самой тройки точно не простые: 6, 9, 12, 15, 18, ...
• Все числа, которые делятся на 4, кроме самой четверки точно не простые: 8, 12, 16, 20, ...
• Все числа, которые делятся на 5, кроме самой пятерки точно не простые: 10, 15, 20, 25, ...
• И так далее

Из этого замечания вырисовывается интересный алгоритм: давайте перебирать число x от 2 до n, и помечать, что числа 2x, 3x, 4x, ... - составные. Если мы ни разу не пометили число, то оно является простым. Это и есть решето Эратосфена.

// isPrime[i] - хранит bool в зависимости от того, что i простое или нет
// Изначального IsPrime - хранит n значений true
isPrime[1] = false; // 1 - не простое число
for (int x = 2; x <= n; x++) {
    // перебираем x
    for (int y = 2 * x; y <= n; y += x) {
        // перебираем числа вида 2x, 3x, 4x
        isPrime[y] = false;
    }
}
// после этого массив isPrime хранит корректные значения

Сперва может показаться, что этот алгоритм работает долго, но это не так. Оценка временной сложности алгоритма:

Это так, потому что сумма гармонического ряда сходится к O(log n).

Давайте сделаем небольшую оптимазацию нашего алгоритма, а именно достаточно перебирать только простые x, и для них помечать 2x, 3x, 4x, ...

// isPrime[i] - хранит bool в зависимости от того, что i простое или нет
// Изначального IsPrime - хранит n значений true
isPrime[1] = false; // 1 - не простое число
for (int x = 2; x <= n; x++) {
    // перебираем x
    if (isPrime[x]) {
        // берем только простые x
        for (int y = 2 * x; y <= n; y += x) {
            // перебираем числа вида 2x, 3x, 4x
            isPrime[y] = false;
        }
    }
}
// после этого массив isPrime хранит корректные значения

Используя математику можно доказать, что это работает за O(n log (logn)). При желании с доказательсвом можно ознакомиться в статьях, отмеченных внизу.

Вернемся к исходной задаче. Мне нужно факторизовать все числа от 1 до n.

Пусть я для каждого числа x знаю его минимальный простой делитель. То есть для 4 - это 2, 21 - это 3, 17 - это 17. Единицу проигнорируем. Обозначим минимальный простой делитель x - minPrimeDivisor[x].

Тогда, чтобы разложить на простые множители число x, сделаем следующее: добавим в список его простых множителей minPrimeDivisor[x], и разложим число x / minPrimeDivisor[x]. Будем повторять этот процесс, пока не получим единицу. Достаточно понятно, что алгоритм корректен. Код:

for (int x = 2; x <= n; x++) {
    // перебираем число x, которое хотим разложить
    vector<int> factors; // список простых множителей
    int current = x; // текущее значение x
    while (current != 1) {
         factors.push_back(minPrimeDivisor[current]);
         current /= minPrimeDivisor[current];
    }
    // factor - корректный список множителей
}

Тогда, чтобы решить исходную задачу, нужно научиться строить массив minPrimeDivisor.

Заметим, что этот массив можно поддерживать с помощью решета Эратосфена. Действительно, давайте изначально заполним массив minPrimeDivisor[i] = i. Теперь, когда я перебираю простое x и перебираю y = 2x, 3x, 4x, ..., прорелаксируем minPrimeDivisor[y] с x. Достаточно понятно, что для числа я переберу все его простые делители, а значит минимум будет посчитан корректно. Код:

// isPrime[i] - хранит bool в зависимости от того, что i простое или нет
// Изначального IsPrime - хранит n значений true
// minPrimeDivisor[i] - хранит минимальный простой делитель i
// Изначально minPrimeDivisor[i] = i
isPrime[1] = false; // 1 - не простое число
for (int x = 2; x <= n; x++) {
    // перебираем x
    if (isPrime[x]) {
        // берем только простые x
        for (int y = 2 * x; y <= n; y += x) {
            // перебираем числа вида 2x, 3x, 4x
            isPrime[y] = false;
            // релаксируем minPrimeDivisor[y] с x
            minPrimeDivisor[y] = min(minPrimeDivisor[y], x);
        }
    }
}
// после этого массивы isPrime, minPrimeDivisor хранят корректные значения

Также можно сократить код:

for (int x = 2; x <= n; x++) {
    if (minPrimeDivisor[x] == x) {
        for (int y = 2 * x; y <= n; y += x) {
            minPrimeDivisor[y] = min(minPrimeDivisor[y], x);
        }
    }
}

Итого мы научились решать задачу, которую хотели за O(nlog(logn)).

Другие статьи:

E-maxx - описание и реализация

Алгоритмика - в статье также содержится оптимазация алгоритма O(n).