wiki/pandoc-pages
1
0
Fork 0
Code Issues Pull Requests Projects Releases Wiki Activity
master
pandoc-pages/pages/linear_nt/main.md

13 KiB
Raw Permalink Blame History

Линейные алгоритмы в теории чисел

Зачастую в теории чисел сложно сказать что-то про одно число, а про несколько значительно проще.
Так, например, до недавнего времени человечеству не было известно, как определить, является ли число простым за полиномиальное время, а найти все простые числа от 1 до $N$ можно решетом Эратосфена, и даже за линейное время.

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

Результаты аналитической теории чисел

Которые, мы, однако, доказывать не будем, но будем активно использовать(здесь \log -- натуральный логарифм)

Вторая теорема Мертенса. $$\sum\limits_{p \leqslant n, p \text{ -- простое}} \frac{1}{p} = \log \log n + C + o(1)$$ Доказательство более слабого утверждения можно найти здесь

Некоторые из её следствий, которыми мы будем пользоваться:

  • \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, и "просеиваем" им, то есть помечаем все числа, которые на него делятся как составные. Его псевдо-реализация:

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.

Его псевдо-реализация:

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++

// 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]
        }
    }
}

Полезные ссылки и источники