MAXimal

добавлено: 10 Jun 2008 17:57
редактировано: 23 Mar 2012 3:53

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

Решето Эратосфена — это алгоритм, позволяющий найти все простые числа в отрезке $[1; n]$ за $O (n \log \log n)$ операций.

Идея проста — запишем ряд чисел $1 \ldots n$ , и будем вычеркивать сначала все числа, делящиеся на $2$ , кроме самого числа $2$ , затем деляющиеся на $3$ , кроме самого числа $3$ , затем на $5$ , затем на $7$ , $11$ , и все остальные простые до $n$ .

Реализация

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

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;

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

При этом мы начинаем идти от $i^2$ , поскольку все меньшие числа, кратные $i$ , обязательно имеют простой делитель меньше $i$ , а значит, все они уже были отсеяны раньше. (Но поскольку $i^2$ легко может переполнить тип $int$ , в коде перед вторым вложенным циклом делается дополнительная проверка с использованием типа $long~long$ .)

При такой реализации алгоритм потребляет $O (n)$ памяти (что очевидно) и выполняет $O (n \log \log n)$ действий (это доказывается в следующем разделе).

Асимптотика

Докажем, что асимптотика алгоритма равна $O (n \log \log n)$ .

Итак, для каждого простого $p \le n$ будет выполняться внутренний цикл, который совершит $\frac{n}{p}$ действий. Следовательно, нам нужно оценить следующую величину:

$\sum_{ \scriptstyle p \le n, \atop\scriptstyle p~[...]$

Вспомним здесь два известных факта: что число простых, меньше либо равных $n$ , приблизительно равно $\frac{n}{\ln n}$ , и что $k$ -ое простое число приблизительно равно $k \ln k$ (это следует из первого утверждения). Тогда сумму можно записать таким образом:

$\sum_{ \scriptstyle p \le n, \atop\scriptstyle p~[...]$

Здесь мы выделили первое простое из суммы, поскольку при $k = 1$ согласно приближению $k \ln k$ получится $0$ , что приведёт к делению на нуль.

Теперь оценим такую сумму с помощью интеграла от той же функции по $k$ от $2$ до $\frac{n}{\ln n}$ (мы можем производить такое приближение, поскольку, фактически, сумма относится к интегралу как его приближение по формуле прямоугольников):

$\sum_{k=2}^{\frac{n}{\ln n}} \frac{1}{k \ln k} \a[...]$

Первообразная для подынтегральной функции есть $\ln \ln k$ . Выполняя подстановку и убирая члены меньшего порядка, получаем:

$\int\limits_2^{\frac{n}{\ln n}} \frac{1}{k \ln k}[...]$

Теперь, возвращаясь к первоначальной сумме, получаем её приближённую оценку:

$\sum_{ \scriptstyle p \le n, \atop\scriptstyle p~[...]$

что и требовалось доказать.

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

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

Самый большой недостаток алгоритма — то, что он "гуляет" по памяти, постоянно выходя за пределы кэш-памяти, из-за чего константа, скрытая в $O (n \log \log n)$ , сравнительно велика.

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

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

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

Самый очевидный момент — что для того, чтобы найти все простые до $n$ , достаточно выполнить просеивание только простыми, не превосходящими корня из $n$ .

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

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

На асимптотику такая оптимизация не влияет (действительно, повторив приведённое выше доказательство, мы получим оценку $n \ln \ln \sqrt{n} + o(n)$ , что, по свойствам логарифма, асимптотически есть то же самое), хотя число операций заметно уменьшится.

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

Поскольку все чётные числа, кроме $2$ , — составные, то можно вообще не обрабатывать никак чётные числа, а оперировать только нечётными числами.

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

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

Заметим, что алгоритм Эратосфена фактически оперирует с $n$ битами памяти. Следовательно, можно существенно сэкономить потребление памяти, храня не $n$ байт — переменных булевского типа, а $n$ бит, т.е. $n/8$ байт памяти.

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

Таким образом, этот подход оправдан, только если $n$ настолько большое, что $n$ байт памяти выделить уже нельзя. Сэкономив память (в $8$ раз), мы заплатим за это существенным замедлением алгоритма.

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

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

Из оптимизации "просеивание простыми до корня" следует, что нет необходимости хранить всё время весь массив $prime[1 \ldots n]$ . Для выполнения просеивания достаточно хранить только простые до корня из $n$ , т.е. $prime[1 \ldots \sqrt{n}]$ , а остальную часть массива $prime$ строить поблочно, храня в текущий момент времени только один блок.

Пусть $s$ — константа, определяющая размер блока, тогда всего будет $\left\lceil \frac{n}{s} \right\rceil$ блоков, $k$ -ый блок ( $k = 0 \ldots \left\lfloor \frac{n}{s} \right\rfloor$ ) содержит числа в отрезке $[ks; ks+s-1]$ . Будем обрабатывать блоки по очереди, т.е. для каждого $k$ -го блока будем перебирать все простые (от $1$ до $\sqrt{n}$ ) и выполнять ими просеивание только внутри текущего блока. Аккуратно стоит обрабатывать первый блок — во-первых, простые из $[1; \sqrt{n}]$ не должны удалить сами себя, а во-вторых, числа $0$ и $1$ должны особо помечаться как не простые. При обработке последнего блока также следует не забывать о том, что последнее нужное число $n$ не обязательно находится в конце блока.

Приведём реализацию блочного решета. Программа считывает число $n$ и находит количество простых от $1$ до $n$ :

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

Асимптотика блочного решета такая же, как и обычного решета Эратосфена (если, конечно, размер $s$ блоков не будет совсем маленьким), зато объём используемой памяти сократится до $O(\sqrt{n} + s)$ и уменьшится "блуждание" по памяти. Но, с другой стороны, для каждого блока для каждого простого из $[1; \sqrt{n}]$ будет выполняться деление, что будет сильно сказываться при меньших размерах блока. Следовательно, при выборе константы $s$ необходимо соблюсти баланс.

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

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

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