E004. Sieve of Eratosthenes

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

Sieve of Eratosthenes — это Algorithme, позволяющий find все простые числа в отрезке

за

операций.

Идея проста — запишем ряд чисел

, и будем вычеркивать сначала все числа, делящиеся на

, кроме самого

числа

, затем деляющиеся на

, кроме самого числа

, затем на

, затем на

,

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

.

Implémentation

Сразу приведём реализацию Algorithmeа:

int n;

vector<char> prime (n+1, true);

prime[0] = prime[1] = false;

for (int i=2; i<=n; ++i)

if (prime[i])

if (i * 1ll * i <= n)

for (int j=i*i; j<=n; j+=i)

prime[j] = false;

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

до

, и, если текущее number

простое,

то помечаем все числа, кратные ему, как составные.

При этом мы начинаем идти от

, поскольку все меньшие числа, кратные

, обязательно имеют простой делитель

меньше

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

легко может переполнить тип

, в коде

перед вторым вложенным циклом делается дополнительная проверка с использованием типа

.)

При такой реализации Algorithme потребляет

памяти (что очевидно) и выполняет

действий

(это доказывается в следующем разделе).

Asymptotic complexity

Докажем, что Asymptotic complexity Algorithmeа равна

.

Итак, для каждого простого

будет выполняться внутренний цикл, который совершит

действий. Следовательно, нам нужно оценить следующую величину: Вспомним здесь два известных факта: что number простых, меньше либо равных

, приблизительно равно

, и что

-

ое простое number приблизительно равно

(это следует из первого утверждения). Тогда сумму можно записать таким образом:

Здесь мы выделили первое простое из суммы, поскольку при

согласно приближению

получится

,

что приведёт к делению на нуль. Теперь оценим такую сумму с помощью интеграла от той же функции по

от

до

(мы можем производить

такое приближение, поскольку, фактически, сумма относится к интегралу как его приближение по формуле прямоугольников):

Первообразная для подынтегральной функции есть

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

порядка, получаем: Теперь, возвращаясь к первоначальной сумме, получаем её приближённую оценку: что и требовалось доказать. Более строгое Preuve (и дающее более точную оценку с точностью до константных множителей) можно find в книге Hardy и Wright "An Introduction to the Theory of Numbers" (стр. 349).

Различные оптимизации решета Эратосфена

Самый большой недостаток Algorithmeа — то, что он "гуляет" по памяти, постоянно Sortieя за пределы кэш-памяти, из-

за чего константа, скрытая в

, сравнительно велика.

Кроме того, для достаточно больших

узким местом становится объём потребляемой памяти. Ниже рассмотрены методы, позволяющие как уменьшить number выполняемых операций, так и значительно сократить потребление памяти.

Просеивание простыми до корня

Самый очевидный момент — что для того, чтобы find все простые до

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

только простыми, не превосходящими корня из

. Таким образом, изменится внешний цикл Algorithmeа:

for (int i=2; i*i<=n; ++i)

На асимптотику такая оптимизация не влияет (действительно, повторив приведённое выше Preuve, мы

получим оценку

, что, по свойствам логарифма, асимптотически есть то же самое), хотя number операций заметно уменьшится.

Решето только по нечётным числам

Поскольку все чётные числа, кроме

, — составные, то можно вообще не обрабатывать никак чётные числа, а оперировать только нечётными числами. Во-первых, это позволит вдвое сократить объём требуемой памяти. Во-вторых, это уменьшит number делаемых Algorithmeом операций Exempleно вдвое.

Уменьшение объёма потребляемой памяти

Заметим, что Algorithme Эратосфена фактически оперирует с

битами памяти. Следовательно, можно

существенно сэкономить потребление памяти, храня не

байт — переменных булевского типа, а

бит, т.е.

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

Таким образом, этот подход оправдан, только если

настолько большое, что

байт памяти выделить уже

нельзя. Сэкономив память (в

раз), мы заплатим за это существенным замедлением Algorithmeа. В завершение стоит отметить, что в языке C++ уже реализованы контейнеры, автоматически осуществляющие битовое сжатие: vector<bool> и bitset<>. Впрочем, если скорость работы очень важна, то лучше реализовать битовое сжатие вручную, с помощью битовых операций — на сегодняшний день компиляторы всё же не в состоянии генерировать достаточно быстрый код.

Блочное решето

Из оптимизации "просеивание простыми до корня" следует, что нет необходимости хранить всё время весь

tableau

. Для выполнения просеивания достаточно хранить только простые до корня из , т. е.

, а остальную часть tableauа

строить поблочно, храня в текущий момент времени

только один блок.

Пусть

— константа, определяющая размер блока, тогда всего будет

блоков,

-ый блок

(

) содержит числа в отрезке

. Будем обрабатывать блоки по очереди, т.е.

для каждого

-го блока будем перебирать все простые (от

до

) и выполнять ими просеивание только

внутри текущего блока. Аккуратно стоит обрабатывать первый блок — во-первых, простые из

не должны

удалить сами себя, а во-вторых, числа

и

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

Приведём реализацию блочного решета. Программа считывает number

и находит количество простых от

до

:

const int SQRT_MAXN = 100000; // корень из максимального значения N

const int S = 10000;

bool nprime[SQRT_MAXN], bl[S];

int primes[SQRT_MAXN], cnt;

int main() {

int n;

cin >> n;

int nsqrt = (int) sqrt (n + .0);

for (int i=2; i<=nsqrt; ++i)

if (!nprime[i]) {

primes[cnt++] = i;

if (i * 1ll * i <= nsqrt)

for (int j=i*i; j<=nsqrt; j+=i)

nprime[j] = true;

}

int result = 0;

for (int k=0, maxk=n/S; k<=maxk; ++k) {

memset (bl, 0, sizeof bl);

int start = k * S;

for (int i=0; i<cnt; ++i) {

int start_idx = (start + primes[i] - 1) / primes[i];

int j = max(start_idx,2) * primes[i] - start;

for (; j<S; j+=primes[i])

bl[j] = true;

}

if (k == 0)

bl[0] = bl[1] = true;

for (int i=0; i<S && start+i<=n; ++i)

if (!bl[i])

++result;

}

cout << result;

} Asymptotic complexity блочного решета такая же, как и обычного решета Эратосфена (если, конечно, размер

блоков не

будет совсем маленьким), зато объём используемой памяти сократится до

и уменьшится "блуждание"

по памяти. Но, с другой стороны, для каждого блока для каждого простого из

будет выполняться деление,

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

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

имеет значение приблизительно от

до

.

Улучшение до линейного времени работы

Algorithme Эратосфена можно преобразовать в другой Algorithme, который уже будет работать за линейное время — см. статью "Sieve of Eratosthenes с линейным временем работы". (Впрочем, этот Algorithme имеет и недостатки.)

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.
    int n;
    List<char> prime (n+1, true);
    prime[0] = prime[1] = false;
    for (int i=2; i<=n; ++i)
            if (prime[i])
                    if (i * 1ll * i <= n)
                            for (int j=i*i; j<=n; j+=i)
                                    prime[j] = false;
    for (int i=2; i*i<=n; ++i)
    const int SQRT_MAXN = 100000; // корень из максимального значения N
    const int S = 10000;
    bool[] nprime = new bool[SQRT_MAXN], bl[S];
    int primes[SQRT_MAXN], cnt;
    int main() {
            int n;
            cin >> n;
            int nsqrt = (int) sqrt (n + .0);
            for (int i=2; i<=nsqrt; ++i)
                    if (!nprime[i]) {
                            primes[cnt++] = i;
                            if (i * 1ll * i <= nsqrt)
                                    for (int j=i*i; j<=nsqrt; j+=i)
                                            nprime[j] = true;
                    }
            int result = 0;
            for (int k=0, maxk=n/S; k<=maxk; ++k) {
                    memset (bl, 0, sizeof bl);
                    int start = k * S;
                    for (int i=0; i<cnt; ++i) {
                            int start_idx = (start + primes[i] - 1) / primes[i];
                            int j = max(start_idx,2) * primes[i] - start;
                            for (; j<S; j+=primes[i])
                                    bl[j] = true;
                    }
                    if (k == 0)
                            bl[0] = bl[1] = true;
                    for (int i=0; i<S && start+i<=n; ++i)
                            if (!bl[i])
                                    ++result;
            }
            Console.WriteLine( result;
    }
}

int n;
vector<char> prime (n+1, true);
prime[0] = prime[1] = false;
for (int i=2; i<=n; ++i)
        if (prime[i])
                if (i * 1ll * i <= n)
                        for (int j=i*i; j<=n; j+=i)
                                prime[j] = false;
for (int i=2; i*i<=n; ++i)
const int SQRT_MAXN = 100000; // корень из максимального значения N
const int S = 10000;
bool nprime[SQRT_MAXN], bl[S];
int primes[SQRT_MAXN], cnt;
int main() {
        int n;
        cin >> n;
        int nsqrt = (int) sqrt (n + .0);
        for (int i=2; i<=nsqrt; ++i)
                if (!nprime[i]) {
                        primes[cnt++] = i;
                        if (i * 1ll * i <= nsqrt)
                                for (int j=i*i; j<=nsqrt; j+=i)
                                        nprime[j] = true;
                }
        int result = 0;
        for (int k=0, maxk=n/S; k<=maxk; ++k) {
                memset (bl, 0, sizeof bl);
                int start = k * S;
                for (int i=0; i<cnt; ++i) {
                        int start_idx = (start + primes[i] - 1) / primes[i];
                        int j = max(start_idx,2) * primes[i] - start;
                        for (; j<S; j+=primes[i])
                                bl[j] = true;
                }
                if (k == 0)
                        bl[0] = bl[1] = true;
                for (int i=0; i<S && start+i<=n; ++i)
                        if (!bl[i])
                                ++result;
        }
        cout << result;
}

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.
    int n;
    ArrayList<Character> prime (n+1, true);
    prime[0] = prime[1] = false;
    for (int i=2; i<=n; ++i)
            if (prime[i])
                    if (i * 1ll * i <= n)
                            for (int j=i*i; j<=n; j+=i)
                                    prime[j] = false;
    for (int i=2; i*i<=n; ++i)
    const int SQRT_MAXN = 100000; // корень из максимального значения N
    const int S = 10000;
    boolean nprime[SQRT_MAXN], bl[S];
    int primes[SQRT_MAXN], cnt;
    int main() {
            int n;
            cin >> n;
            int nsqrt = (int) sqrt (n + .0);
            for (int i=2; i<=nsqrt; ++i)
                    if (!nprime[i]) {
                            primes[cnt++] = i;
                            if (i * 1ll * i <= nsqrt)
                                    for (int j=i*i; j<=nsqrt; j+=i)
                                            nprime[j] = true;
                    }
            int result = 0;
            for (int k=0, maxk=n/S; k<=maxk; ++k) {
                    memset (bl, 0, sizeof bl);
                    int start = k * S;
                    for (int i=0; i<cnt; ++i) {
                            int start_idx = (start + primes[i] - 1) / primes[i];
                            int j = max(start_idx,2) * primes[i] - start;
                            for (; j<S; j+=primes[i])
                                    bl[j] = true;
                    }
                    if (k == 0)
                            bl[0] = bl[1] = true;
                    for (int i=0; i<S && start+i<=n; ++i)
                            if (!bl[i])
                                    ++result;
            }
            System.out.println( result;
    }
}