202 lines
13 KiB
Markdown
202 lines
13 KiB
Markdown
# Линейные алгоритмы в теории чисел
|
||
|
||
Зачастую в теории чисел сложно сказать что-то про одно число, а про несколько значительно проще.
|
||
Так, например, до [недавнего времени](https://en.wikipedia.org/wiki/AKS_primality_test) человечеству не было
|
||
известно, как определить, является ли число простым за полиномиальное время, а найти все простые числа от $1$ до $N$
|
||
можно решетом Эратосфена, и даже за линейное время.
|
||
|
||
**Важное замечание:** далее во всех асимптотиках алгоритмов мы считаем количество арифметических операций. Иначе говоря, мы не берем во внимание, сколько выполняется операция умножения или деления, взятия по модулю.
|
||
|
||
## Результаты аналитической теории чисел
|
||
|
||
Которые, мы, однако, доказывать не будем, но будем активно использовать(здесь $\log$ -- натуральный логарифм)
|
||
|
||
**Вторая теорема Мертенса**. $$\sum\limits_{p \leqslant n, p \text{ -- простое}} \frac{1}{p} = \log \log n + C + o(1)$$
|
||
Доказательство более слабого утверждения можно найти
|
||
[здесь](https://math.stackexchange.com/questions/2678885/easy-proof-of-a-weak-form-of-mertenss-second-theorem)
|
||
|
||
Некоторые из её следствий, которыми мы будем пользоваться:
|
||
|
||
* $\sum\limits_{x \leqslant n} \omega(x) = n \log \log n + C_1n + o(n)$, где $\omega(n)$ -- количество различных простых делителей
|
||
* $\sum\limits_{x \leqslant n} \Omega(x) = n \log \log n + C_2n + o(n)$, где $\Omega(n)$ -- количество простых делителей в разложении
|
||
|
||
На сумму количества делителей известна следующая оценка: $\sum\limits_{x \leqslant n} \sigma_0(x) = n \log n + Cn + O\left(\sqrt{n}\right)$.
|
||
|
||
## Решето Эратосфена
|
||
|
||
Самый простой вариант решета Эратосфена заключается в том, что мы берём каждое число от $2$ до $N$, и "просеиваем" им, то есть помечаем
|
||
все числа, которые на него делятся как составные. Его псевдо-реализация:
|
||
``` {.c}
|
||
for (int i = 2; i <= N; i++) {
|
||
for (int k = 2; i * k <= N; k++) {
|
||
primes[i * k] = false;
|
||
}
|
||
}
|
||
```
|
||
|
||
Нетрудно оценить его время работы. Когда мы просеваем числом $i$, мы делаем $\frac{N}{i}$ операций, то есть суммарно
|
||
$\frac{N}{1} + \frac{N}{2} + \ldots + \frac{N}{N} = N \left(1 + \frac{1}{2} + \ldots + \frac{1}{N}\right)$,
|
||
что легко оценивается сверху как $O(N \log N)$ суммой гармонического ряда.
|
||
|
||
Перед нами встаёт вопрос, как улучшить наш алгоритм. Первое, что приходит в голову, просеивать только числами до $\sqrt{N}$, потому что
|
||
если число $x \leqslant N$ делится на $k \geqslant \sqrt{N}$, то оно делится на $\frac{d}{k} \leqslant \sqrt{N}$.
|
||
|
||
Однако, это никак не улучшает асимптотику, поскольку $\log \left(\sqrt{N}\right) = \frac{1}{2} \log N$.
|
||
|
||
Ещё одна идея -- просеивать только простыми числами, то есть теми, которыми алгоритм ещё не пометил. И, оказывается, что эта идея
|
||
улучшает асимптотику до $O(N \log \log N)$. Действительно, теперь мы делаем
|
||
$\sum\limits_{p \leqslant N} \frac{N}{p}$ операций, где $p$ -- простые
|
||
числа. Из второй теоремы Мертенса это и есть $O(N \log \log N)$.
|
||
|
||
### Линейное решето Эратосфена
|
||
|
||
Решим более общую задачу. Найдём для каждого числа от $1$ до $N$ его наименьший делитель, а также будем поддерживать список
|
||
уже найденных простых чисел.
|
||
|
||
Рассмотрим следующий алгоритм, на $i$-ом шаге которого мы:
|
||
|
||
1. Если наименьший делитель числа $i$ не был найден до этого, добавляем $i$ в список простых.
|
||
2. Для всех простых чисел $p$, не больших $q$, наименьшего делителя $i$, объявляем наименьшим делителем числа $pi$ число $p$.
|
||
|
||
Его псевдо-реализация:
|
||
``` {.c}
|
||
for (int i = 2; i <= N; i++) {
|
||
if (lowest_divisor[i] == 0) {
|
||
primes.push_back(i);
|
||
lowest_divisor[i] = i;
|
||
}
|
||
for (int prime : primes) {
|
||
if (prime * i > N || prime > lowest_divisor[i]) {
|
||
break;
|
||
}
|
||
lowest_divisor[i * prime] = prime;
|
||
}
|
||
}
|
||
```
|
||
|
||
Докажем корректность алгоритма, а именно, что на $i$-ом шаге корректно посчитаны все наименьшие делители чисел, не превосходящих $i$
|
||
, а также найдены все простые числа не превосходящие $i$ по индукции.
|
||
|
||
База очевидна.
|
||
|
||
Переход. Пусть до числа $i$ утверждение выполнено.
|
||
|
||
Если число $i$ -- простое, то ни на каком прошлом шаге $j$ не нашлось бы число $p$, что $jp = i$,
|
||
значит мы добавим его в список простых, и объявим его наименьшим делителем $i$, что и требовалось.
|
||
|
||
Если число $i$ -- составное, тогда его наименьшим делителем является простое число, обозначим его за $p$.
|
||
|
||
Наименьший делитель числа $\frac{i}{p}$ не меньше $p$ из выбора $p$. Значит, на шаге с номером $\frac{i}{p}$ мы объявили наименьшим делителем
|
||
$i$ число $p$.
|
||
|
||
Пусть наименьший делитель $i$ обновлялся ещё и на шаге $j \neq \frac{i}{p}$. Тогда $jq = i$, причём $q \neq p$ -- простое. Пусть
|
||
|
||
* $q < p$. Тогда у числа $i$ есть делитель, меньший $p$, а именно $q$. Противоречие.
|
||
* $q > p$. Тогда у числа $j = \frac{i}{q}$ наименьший делитель равен $p < q$, а значит мы бы не обновляли ответ для $i$ на шаге $j$. Противоречие.
|
||
|
||
|
||
Значит, алгоритм корректно находит наименьший общий делитель всех чисел от $1$ до $N$, причём в процессе каждому числу присваивался
|
||
наименьший делитель ровно один раз.
|
||
Значит, наш алгоритм работает за линейное время.
|
||
|
||
## Вычисление значений мультипликативной функции $f$ для всех чисел от $1$ до $N$ за $O(N)$
|
||
|
||
Будем считать, что наименьшие делители всех чисел уже посчитаны до этого за $O(N)$.
|
||
|
||
Единственное требование к функции $f$ -- умение по простому числу $p$ и числу $k$
|
||
вычислять $f(p^k)$ за $O(k)$.
|
||
|
||
Так, например при $f = \varphi$, функции Эйлера, $f(p^k) = p^{k - 1} \cdot (p - 1)$.
|
||
|
||
Заметим, что $f(n) = f(\frac{n}{p^k}) \cdot f(p^k)$, где $p$ -- наименьший делитель $n$, а $d(n) = k$ -- его степень
|
||
вхождения
|
||
|
||
Из мультипликативности $f(ab) = f(a) \cdot f(b)$ если $\gcd(a, b) = 1$, а $\gcd(p^k, \frac{n}{p^k}) = 1$ из опредения $p^k$.
|
||
|
||
Откуда вытекает следующий алгоритм, на $i$-ом шаге которого мы:
|
||
|
||
1. Считаем для данного $i$ его $p^k$ за $O(d(i))$: делим на наименьший делитель, пока он не меняется.
|
||
2. Если $p^k = i$, считаем $f(p^k)$ за $O(d(i))$ или за $O(1)$, что и будет ответом для числа $i$.
|
||
3. Иначе объявляем $f(i)$ равным $f(p^k) \cdot f(\frac{n}{p^k})$, вычисляя это выражение за $O(1)$, поскольку оно
|
||
использует только значения $f$ для аргументов меньших $i$, для которых всё уже посчитано.
|
||
|
||
Итоговая сложность алгоритма составляет $O\left( \sum\limits_{i \leqslant N} d(i) \right)$ времени и $O(N)$ памяти.
|
||
|
||
**Основная лемма.** $\sum\limits_{i \leqslant N} d(i) = \Theta(N)$.
|
||
|
||
Пусть $n = p_1^{\alpha_1} \cdot \ldots p_k^{\alpha_k}$. Тогда
|
||
$d(n) \leqslant 1 + (\alpha_1 - 1) + \ldots + (\alpha_k - 1) = 1 + \Omega(n) - \omega(n)$.
|
||
|
||
Откуда из теорем в начале для исходной суммы выполнена оценка:
|
||
|
||
$$N + \sum\limits_{x \leqslant N} \Omega(x) - \sum\limits_{x \leqslant N} \omega(x) = N + (C_2 - C_1)N + o(N) = O(N)$$
|
||
|
||
Что и требовалось. То есть на самом деле наш алгоритм линейный, и даже с очень неплохой константой. Также
|
||
он требует всего $2N$ памяти, не считая списка простых(его можно не хранить после вычисления наименьших делителей).
|
||
|
||
Его можно применять для нахождения функции Мёбиуса, Эйлера, количества делителей, суммы делителей, и так далее.
|
||
|
||
Также стоит отметить, что можно было просто предподсчитать $d(i)$ для всех чисел за $O(N)$
|
||
потратив ещё $N$ памяти.
|
||
|
||
Ценность данной техники в том, что она позволяет сэкономить память, а также решать задачу почти в общем случае, даже
|
||
когда считать $f$ для степени простого очень сложно.
|
||
|
||
### Реализация
|
||
|
||
Ниже приведена реализация на языке C++
|
||
|
||
``` {.c}
|
||
// ld -- наименьший делитель числа, предпосчитан линейным решетом
|
||
// fs -- список значений функции f для чисел от $1$ до $N$
|
||
|
||
// Также можно возвращать, например не $k$, а $p^k$ или пару $(k, p^k)$, чтобы не считать дважды
|
||
int d(int n) {
|
||
if (n == 1) {
|
||
return 1;
|
||
}
|
||
int p = ld[n];
|
||
int ans = 1;
|
||
while (ld[n / p] == p) {
|
||
ans += 1;
|
||
n /= p;
|
||
}
|
||
return ans;
|
||
}
|
||
|
||
int f(int p, int k) {
|
||
// Результат функции для p^k
|
||
// В качестве примера возьмём функцию Эйлера
|
||
int ans = (p - 1);
|
||
for (int i = 0; i < k - 1; i++) {
|
||
ans *= p;
|
||
}
|
||
}
|
||
|
||
void calculate(int N) {
|
||
fs[1] = 1;
|
||
for (int i = 2; i <= N; i++) {
|
||
int k = d(i), p = ld[i];
|
||
int pk = 1; // $p^k$
|
||
for (int j = 0; j < k; j++) {
|
||
pk *= p;
|
||
}
|
||
if (pk == i) {
|
||
fs[i] = f(p, k);
|
||
} else {
|
||
fs[i] = fs[i / pk] * fs[pk]
|
||
}
|
||
}
|
||
}
|
||
|
||
```
|
||
|
||
## Полезные ссылки и источники
|
||
|
||
* [Оригинальное доказательство основной леммы](https://math.stackexchange.com/questions/3886140/upper-bound-of-sum-of-powers-of-lowest-divisors)
|
||
* [Hardy-Wright: An introduction to the theory of numbers pdf](https://blngcc.files.wordpress.com/2008/11/hardy-wright-theory_of_numbers.pdf)
|
||
, книга со всеми доказательствами.
|
||
* [Шаблоны алгоритмики](https://github.com/algorithmica-org/algorithmica), то, откуда были взяты css-файлы и шрифты для отображения этой
|
||
статьи
|
||
|