Лінійне решето

Дано число $n$; знайти всі прості числа на відрізку $[2;n]$.

Стандартний спосіб розв'язати цю задачу — скористатися решетом Ератосфена. Цей алгоритм дуже простий, але має час роботи $O(n \log \log n)$.

Хоча відомо чимало алгоритмів із сублінійним часом роботи (тобто $o(n)$), описаний нижче алгоритм цікавий своєю простотою: він не складніший за класичне решето Ератосфена.

До того ж наведений тут алгоритм як побічний ефект обчислює розклади на множники всіх чисел на відрізку $[2; n]$, що може стати в пригоді в багатьох практичних застосуваннях.

Слабке місце цього алгоритму — він використовує більше пам'яті, ніж класичне решето Ератосфена: йому потрібен масив із $n$ чисел, тоді як класичному решету Ератосфена достатньо $n$ бітів пам'яті (що у 32 рази менше).

Тому має сенс застосовувати описаний алгоритм лише для чисел порядку $10^7$ і не більших.

Алгоритм належить Полу Прітчарду (Paul Pritchard). Це варіант Algorithm 3.3 із (Pritchard, 1987: див. список літератури в кінці статті).

Коли підходить цей алгоритм?

• Чи потрібен як побічний продукт найменший простий дільник $lp[i]$ кожного числа (для швидкої факторизації)? Якщо потрібен лише список простих — вистачить простішого решета Ератосфена.
• Чи $n \lesssim 10^7$? Алгоритм тримає масив із $n$ int, а не $n$ бітів, тож для більших $n$ пам'яті не вистачить — беріть решето Ератосфена.
• Чи важлива гарантія, що кожне число викреслюється рівно один раз (істинно лінійний час $O(n)$), а не $O(n \log \log n)$?

Алгоритм

Наша мета — обчислити найменший простий дільник $lp [i]$ для кожного числа $i$ на відрізку $[2; n]$.

Окрім того, нам потрібно зберігати список усіх знайдених простих чисел — назвімо його $pr []$.

Значення $lp [i]$ ми ініціалізуємо нулями, що означає припущення, ніби всі числа прості. Під час виконання алгоритму цей масив поступово заповнюватиметься.

Тепер ми проходитимемо числа від 2 до $n$. Для поточного числа $i$ маємо два випадки:

• $lp[i] = 0$ — це означає, що $i$ просте, тобто ми не знайшли для нього жодного меншого дільника.
  Тож ми присвоюємо $lp [i] = i$ і додаємо $i$ в кінець списку $pr[]$.

• $lp[i] \neq 0$ — це означає, що $i$ складене, а його найменший простий дільник — $lp [i]$.

В обох випадках ми оновлюємо значення $lp []$ для чисел, які діляться на $i$. Однак наша мета — навчитися робити це так, щоб задавати значення $lp []$ щонайбільше один раз для кожного числа. Зробити це можна так:

Розгляньмо числа $x_j = i \cdot p_j$, де $p_j$ — це всі прості числа, менші або рівні $lp [i]$ (саме тому нам потрібно зберігати список усіх простих чисел).

Для всіх чисел такого вигляду ми задамо нове значення $lp [x_j] = p_j$.

Доведення коректності цього алгоритму та його час роботи можна знайти після реалізації.

Реалізація

C++

const int N = 10000000;
vector<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; i * pr[j] <= N; ++j) {
		lp[i * pr[j]] = pr[j];
		if (pr[j] == lp[i]) {
			break;
		}
	}
}

Python

N = 10_000_000
lp = [0] * (N + 1)  # найменший простий дільник кожного числа
pr = []             # список знайдених простих чисел

for i in range(2, N + 1):
    if lp[i] == 0:          # i просте
        lp[i] = i
        pr.append(i)
    for p in pr:            # перебираємо прості p_j <= lp[i]
        if i * p > N:
            break
        lp[i * p] = p
        if p == lp[i]:      # дійшли до lp[i] — далі не йдемо
            break

TypeScript

const N = 10_000_000;
const lp: number[] = new Array(N + 1).fill(0); // найменший простий дільник
const pr: number[] = [];                       // список простих чисел

for (let i = 2; i <= N; ++i) {
  if (lp[i] === 0) {        // i просте
    lp[i] = i;
    pr.push(i);
  }
  for (let j = 0; i * pr[j] <= N; ++j) { // прості pr[j] <= lp[i]
    lp[i * pr[j]] = pr[j];
    if (pr[j] === lp[i]) break;          // дійшли до lp[i]
  }
}

Go

const N = 10_000_000
lp := make([]int, N+1) // найменший простий дільник кожного числа
pr := []int{}          // список знайдених простих чисел

for i := 2; i <= N; i++ {
	if lp[i] == 0 { // i просте
		lp[i] = i
		pr = append(pr, i)
	}
	for j := 0; i*pr[j] <= N; j++ { // прості pr[j] <= lp[i]
		lp[i*pr[j]] = pr[j]
		if pr[j] == lp[i] { // дійшли до lp[i]
			break
		}
	}
}

Доведення коректності

Нам потрібно довести, що алгоритм задає всі значення $lp []$ правильно і що кожне значення буде задано рівно один раз. Тоді алгоритм матиме лінійний час роботи, бо всі інші дії алгоритму, очевидно, виконуються за $O (n)$.

Зауважимо, що кожне число $i$ має рівно одне подання у вигляді:

$$i = lp [i] \cdot x,$$

де $lp [i]$ — найменший простий дільник $i$, а число $x$ не має жодних простих дільників, менших за $lp [i]$, тобто

$$lp [i] \le lp [x].$$

Тепер порівняймо це з діями нашого алгоритму: насправді для кожного $x$ він проходить усі прості числа, на які його можна було б помножити, тобто всі прості числа аж до $lp [x]$ включно, щоб отримати числа у наведеному вище вигляді.

Отже, алгоритм пройде кожне складене число рівно один раз, задаючи там правильні значення $lp []$. Що й треба було довести.

Час роботи та пам'ять

Хоча час роботи $O(n)$ кращий за $O(n \log \log n)$ класичного решета Ератосфена, різниця між ними не така вже й велика. На практиці лінійне решето працює приблизно так само швидко, як і типова реалізація решета Ератосфена.

Порівняно з оптимізованими версіями решета Ератосфена, наприклад із сегментованим решетом, воно значно повільніше.

Якщо взяти до уваги вимоги цього алгоритму до пам'яті — масив $lp []$ довжини $n$ та масив $pr []$ довжини $\frac n {\ln n}$, — то цей алгоритм видається гіршим за класичне решето в усьому.

Однак його рятівна риса в тому, що цей алгоритм обчислює масив $lp []$, який дозволяє нам знаходити розклад будь-якого числа на відрізку $[2; n]$ за час порядку розміру цього розкладу. Більше того, лише за допомогою одного додаткового масиву ми можемо уникнути ділень під час пошуку розкладу.

Знання розкладів усіх чисел на множники дуже корисне для деяких задач, і цей алгоритм — один із небагатьох, що дозволяють знаходити їх за лінійний час.

Список літератури

• Paul Pritchard, Linear Prime-Number Sieves: a Family Tree, Science of Computer Programming, vol. 9 (1987), pp.17-35.