В комментариях к одному из прошлых постов о решете Эратосфена был упомянут этот короткий алгоритм из Википедии:

Алгоритм 1:

1: для i := 2, 3, 4, ..., до n: 
2:  если lp[i] = 0:
3:       lp[i] := i
4:       pr[] += {i}
5:   для p из pr пока p ≤ lp[i] и p*i ≤ n:
6:       lp[p*i] := p
Результат:
 lp - минимальный простой делитель для кажого числа до n
 pr - список всех простых до n.

Алгоритм простой, но не всем он показался очевидным. Главная же проблема в том, что на Википедии нет доказательства, а ссылка на первоисточник (pdf) содержит довольно сильно отличающийся от приведенного выше алгоритм.

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

Примечание о сложности алгоритма
Хоть этот алгоритм и асимптотически быстрее стандартного решета Эратосфена за O(n log log n), ему требуется гораздо больше памяти. Поэтому для по-настоящему больших n, где бы этот алгоритм засиял во всей красе, он не применим. Однако, он очень полезен даже для небольших n, если вам надо не только найти все простые числа, но и быстро раскладывать на множители все числа до n.

Определения


$\mathcal{P}$ — множество простых чисел.
$lp(x)$ — минимальный простой делитель числа: $lp(x) = min \{p \in \mathcal{P}\ \vert\ p\vert x \}$
$rd(x)$ — остальные множители: $rd(x) = x / lp(x)$

Обратите внимание, все определения выше даны для $x \ge 2$.

Некоторые очевидные свойства введенных выше функций, которые будут использоваться дальше:
$lp(a\times b) = min(lp(a), lp(b))$
$rd(x) < x$
$p \in \mathcal{P} => lp(p) = p$

Доказательство


Лемма 1: $lp(x) \le lp(rd(x)), \forall x \notin \mathcal{P}, x>1$
Доказательство: Т.к. любой делитель $rd(x)$ также является делителем $x$, а $lp(x)$ не превосходит любой делитель $x$, то $lp(x)$ не превосходит и любой делитель $rd(x)$, включая наименьший из них. Единственная проблема в этом утверждении, если $rd(x)$ не имеет простых делителей, но это возможно только в случае $x \in \mathcal{P}$, который исключен в условии леммы.
$\square$

Пусть $E = \{(lp(x), rd(x)) \vert \forall x \notin \mathcal{P}, 2\le x \le n\}$

Т.к. $lp(x)\times rd(x) = x$ (по определению $rd()$), если бы нам было дано множество $E$, то мы смогли бы вычислить $lp()$ для всех составных чисел до n. Это делает, например, следующий алгоритм:

Алгоритм 2:

1: Для всех (l,r) из E:
2:    lp[l*r] := l;

Заметим, что $\vert E \vert \le n$, поэтому алгоритм 2 выше работает за линейную сложность.

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

Пусть $E' = \{(p, r) \vert p \in \mathcal{P}, 2\le r < n, p \le lp(r), p \times r \le n\}$

Лемма 2: $E = E'$
Доказательство:

Пусть $(a,b)\in E$ => $\exists x \notin \mathcal{P},\ 2\le x\le n \mid a=lp(x), b=rd(x)$.

По определению $lp(), rd()$: $a \in \mathcal{P}$, $a\times b = x$. По лемме 1, $a \le lp(b)$.
т.к. $rd(x) < x$, то $b < n$.
Поскольку $x \notin \mathcal{P}$, $b \ge 2$.
Так же, по определению $E$, $x \le n$, следовательно, $a \times b \le n$.
Все 4 условия $E'$ выполнены, значит, $(a,b) \in E'$ => $E \subset E'$.

Пусть $(a,b)\in E'$. Пусть $x = a \times b$.
По определению $E'$, $a \in \mathcal{P}$. Значит, $a$ — простой делитеть $x$.
$lp(x) = lp(a\times b) = min(lp(a), lp(b)) = min(a, lp(b)).$
Т.к $a \le lp(b)$, то $lp(x) = a.$ Значит, $b = rd(x).$
По определению, $x = a \times b \le n.$ Также, очевидно, $x = a \times b \ge 2.$
Все условия для $E$ выполнены => $E' \subset E.$

Следовательно, $E = E'.$
$\square$

Теперь осталось перебрать правильные $r$ и $p$ из определения множества $E'$ и применить алгоритм 2. Именно это и делает Алгоритм 1 (только вместо r используется переменная i). Но есть проблема! Что бы перебрать все элементы множества $E'$ для вычисления $lp,$ нам надо знать $lp,$ ведь в определении есть условие $p \le lp(r)$.

К счастью, эта проблема обходится правильным порядком перебора. Надо перебирать $r$ во внешнем цикле, а $p$ — во внутреннем. Тогда $lp(r)$ уже будет вычислено. Этот факт доказывает следующая теорема.

Теорема 1:
Алгоритм 1 поддерживает следующий инвариант: После выполнения итерации внешнего цикла при i=k, все простые числа до k включительно будут выделены в список pr. Также будет подсчитано $lp()$ для всех $x \notin \mathcal{P} \mid x\le n,\ rd(x) \le k$ в массиве lp.

Доказательство:
Докажем по индукции. Для k=2 инвариант проверяется вручную. Единственное простое число 2 будет добавлено в список pr, потому что массив lp[] изначально заполнен нулями. Также, единственное составное число у которого $rd(x) \le 2$ — это 4. Несложно убедиться, что внутренний цикл выполнится ровно один раз (при n>3) и правильно выполнит lp[4] := 2.

Теперь допустим, что инвариант выполняется для итерации i=k-1. Докажем, что он будет выполнятся и для итерации i=k.

Для этого достаточно проверить что число i, если оно простое, будет добавлено в список pr и что $lp()$ будет подсчитано для всех составных чисел $x \le n,$ т.ч. $rd(x) = i.$ Именно эти числа из инварианта для k не покрыты уже инвариантом для k-1.

Если i простое, то lp[i] будет равно 0, ведь единственная операция записи в массив lp, которая теоретически могла бы подсчитать lp[i] (в строке 6), всегда пишет в составные индексы, ведь p*i (для i > 1) — всегда составное. Поэтому число i будет добавлено в список простых. Также, в строке 3 будет подсчитано lp[i].

Если же i составное, то на начало итерации lp[i] уже будет подсчитано по инварианту для i=k-1, ведь $rd(i) < i$ или $rd(i) \le i-1,$ следовательно i попадает под условие инварианта в предыдущей итерации. Поэтому составное число i не будет добавлено в список простых чисел.

Далее, имея корректное lp[i] и все простые числа до i включительно цикл в строках 5-6 переберет все элементы $(p, i) \in E'$, у которых вторая часть равна i:

  • $p \in P,$ потому что оно из списка pr
  • $p \le lp(i),$ по условию останова цикла
  • $p\times i \le n,$ по условию останова цикла
  • $i < n,$ следует из $p\times i \le n,\ p > 1$

Все нужные простые числа в списке pr есть, т.к. нужны только числа до $lp(i) \le i$. Поэтому будут подсчитаны значения lp[] для всех составных чисел $x$, у которых $rd(x)=i$. Это ровно все те числа, которых не хватало при переходе от инварианта для k-1 к инварианту для k.

Следовательно, инавриант выполняется для любых i = 2..n.

$\square$

По инварианту из Теоремы 1 для i=n получается, что все простые числа до n и все lp[] будут подсчитаны алгоритмом 1.

Более того, поскольку в строках 5-6 перебираются различные элементы множества $E$, то суммарно внутренний цикл выполнит не более $\vert E \vert < n$ операций. Операция проверки в цикле выполнится ровно $\vert E \vert + n - 1 < 2n$ раз ($\vert E \vert$ раз вернет true и n-1 раз, для каждого i, вернет false). Все оставшиеся операции вложены в один цикл по i от 2 до n.
Отсюда следует, что сложность алгоритма 1 — O(n).