Быстрое преобразование Фурье

Быстрое преобразование Фурье

Рассмотрим такую распространённую операцию как умножение двух целых чисел. Квадратичный алгоритм — умножения в столбик — все знают со школы. Долгое время предполагалось, что ничего быстрее придумать нельзя.

Первым эту гипотезу опроверг Анатолий Карацуба. Его алгоритм сводит умножение двух \(n\)-значиных чисел к трём умножениям \(\frac{n}{2}\)-значных чисел, что даёт оценку времени работы

\[ T(n)=3T\left(\dfrac n 2\right)+O(n)=O\left(n^{\log_2 3}\right)\approx O(n^{1.58}) \]

Чтобы перейти к алгоритму с лучшей оценкой, нам нужно сначала установить несколько фактов о многочленах.

Умножение многочленов

Обратим внимание на то, что любое число можно представить многочленом:

\[ \begin{aligned} A(x) &= a_0 + a_1\cdot x + a_2 \cdot x^2 + \dots + a_n \cdot x^n \\ &= a_0 + a_1\cdot 2 + a_2 \cdot 2^2 + \dots + a_n \cdot 2^n \end{aligned} \]

Основание x при этом может быть выбрано произвольно.

Чтобы перемножить два числа, мы можем перемножить соответствующие им многочлены, а затем произвести каррирование: пройтись от нижних разрядов получившегося многочлена и «сдвинуть» переполнившиеся разряды:

const int base = 10;

vector<int> normalize(vector<int> a) {
    int carry = 0;
    for (int &x : a) {
        x += carry;
        carry = x / base;
        x %= base;
    }
    while (carry > 0) {
        a.push_back(carry % base);
        carry /= base;
    }
    return a;
}

vector<int> multiply(vector<int> a, vector<int> b) {
    return normalize(poly_multiply(a, b));
}

Прямая формула для произведения многочленов имеет вид

\[ \left(\sum_{i=0}^n a_i x^i\right)\cdot\left(\sum_{j=0}^m b_j x^j\right)=\sum_{k=0}^{n+m}x^k\sum_{i+j=k}a_i b_j \]

Её подсчёт требует \(O(n^2)\) операций, что нас не устраивает. Подойдём к этой задаче с другой стороны.

Интерполяция

Теорема. Пусть есть набор различных точек \(x_0, x_1, \dots, x_{n}\). Многочлен степени \(n\) однозначно задаётся своими значениями в этих точках. (Коэффициентов у этого многочлена столько же, сколько и точек — прим. К. О.)

Доказательство будет конструктивным — можно явным образом задать многочлен, который принимает заданные значения \(y_0, y_1, \ldots, y_n\) в этих точках:

\[ y(x)=\sum\limits_{i=0}^{n}y_i\prod\limits_{j\neq i}\dfrac{x-x_j}{x_i-x_j} \]

Корректность. Проверим, что в точке \(x_i\) значение действительно будет равно \(y\):

  1. Для \(i\)-го слагаемого внутреннее произведение будет равно единице, если вместо \(x\) подставить \(x_i\): в этом случае просто перемножается \((n-1)\) единица. Эта единица помножится на \(y_i\) и войдёт в сумму.

  2. Для всех остальных слагаемых произведение занулится: один из множетелей будет равен \((x_i - x_i)\).

Уникальность. Предположим, есть два подходящих многочлена степени \(n\)\(A(x)\) и \(B(x)\). Рассмотрим их разность. В точках \(x_i\) значение получившегося многочлена \(A(x) - B(x)\) будет равняться нулю. Если так, то точки \(x_i\) должны являться его корнями, и тогда разность можно записать так:

\[ A(x) - B(x) = \alpha \prod_{i=0}^n (x-x_i) \]

для какого-то числа \(\alpha\). Тут мы получаем противоречие: если раскрыть это произведение, то получится многочлен степени \(n+1\), который нельзя получить разностью двух многочленов степени \(n\).

Этот многочлен называется интерполяционным многочленом Лагранжа, а сама задача проведения многочлена через точки — интерполяцией.

Примечание. На практике интерполяцию решают методом Гаусса: её можно свести к решению линейного уравнения \(aX = y\), где \(X\) это матрица следующего вида:

\[ \begin{pmatrix} 1 & x_0 & x_0^2 & \ldots & x_0^n \\ 1 & x_1 & x_1^2 & \ldots & x_1^n \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 1 & x_n & x_n^2 & \ldots & x_n^n \\ \end{pmatrix} \]

Важный факт: многочлен можно однозначно задать не только своими коэффициентами, но также корнями и значениями хотя бы в \((n+1)\)-ой точке.

Умножение через интерполяцию

Что происходит со значениями многочлена-произведения \(A(x) B(x)\) в конкретной точке \(x_i\)? Оно просто становится равным \(A(x_i) B(x_i)\).

Основная идея алгоритма: если мы знаем значения в каких-то различных \(n + m\) точках для обоих многочленов \(A\) и \(B\), то, попарно перемножив их, мы за \(O(n + m)\) операций можем получить значения в тех же точках для многочлена \(A(x) B(x)\) — а с их помощью можно интерполяцией получить исходный многочлен и решить задачу.

vector<int> poly_multiply(vector<int> a, vector<int> b) {
    vector<int> A = evaluate(a);
    vector<int> B = evaluate(b);
    for (int i = 0; i < A.size(); i++)
        A[i] *= B[i];
    return interpolate(A);
}

Если притвориться, что evaluate и interpolate работают за линейное время, то умножение тоже будет работать за линейное время.

К сожалению, непосредственное вычисление значений требует \(O(n^2)\) операций, а интерполяция — как методом Гаусса, так и через символьное вычисление многочлена Лагранжа — и того больше, \(O(n^3)\).

Но что, если бы мы могли вычислять значения в точках и делать интерполяцию быстрее?

Комплексные числа

Определение. Комплексные числа — это числа вида \(a + bi\), где \(a\) и \(b\) это обычные вещественные числа, а \(i\) это так называемая мнимая единица: это число, для которого выполняется равенство \(i^2 = -1\).

Комплексные числа ввели в алгебре, чтобы работать с корнями из отрицательных чисел: \(i\) в каком-то смысле равно \(\sqrt{-1}\). Так же, как и отрицательные числа, они как бы «не существуют» в реальном мире, а только в сознании математиков.

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

Комплексная плоскость

Комплексные числа удобно изображать на плоскости в виде вектора \((a, b)\) и считать через них всякую геометрию.

Модулем комплексного числа называется действительное число \(r = \sqrt{a^2 + b^2}\) . Геометрически, это длина вектора \((a, b)\).

Аргументом комплексного числа называется действительное число \(\phi \in (-\pi, \pi]\), для которого выполнено \(\tan \phi = \frac{b}{a}\). Геометрически, это значение угла между \((a, 0)\) и \((a, b)\). Для нуля — вектора \((0, 0)\) — аргумент не определён.

Таким образом комплексное число можно представить в полярных координатах:

\[ a+bi = r ( \cos \phi + i \sin \phi ) \]

Подобное представление удобно по следующей причине: чтобы перемножить два комплексных числа, нужно перемножить их модули и сложить аргументы.

Упражнение. Докажите это.

Формула эйлера

Определим число Эйлера \(e\) как число со следующим свойством:

\[ e^{i\phi} = \cos \phi + i \sin \phi \]

Просто введём такую нотацию для выражения \(\cos \phi + i \sin \phi\). Не надо думать, почему это так.

Геометрически, все такие точки живут на единичном круге:

Такая нотация удобна, потому что можно обращаться с \(e^{i\phi}\) как с обычной экспонентой. Пусть мы, например, хотим перемножить два числа на единичном круге с аргументами \(a\) и \(b\). Тогда это можно записать так:

\[ (\cos a + i \sin a) \cdot (\cos b + i \sin b) = e^{i (a+b)} \]

Упражнение. Проверьте это: раскройте скобки и проделайте немного алгебры.

Корни из единицы

У комплексных чисел есть много других замечательных свойств, но нам для алгоритма на самом деле потребуется только следующее:

Утверждение. Для любого натурального \(n\) есть ровно \(n\) комплексных «корней из единицы», то есть чисел \(w_k\), для которых выполнено:

\[ w_k^n = 1 \]

А именно, это будут числа вида:

\[ w_k = e^{i \tau \frac{k}{n}} \]

где \(\tau\) обозначает \(2 \pi\), «целый круг». Это довольно новая нотация.

На комплексной плоскости эти числа располагаются на единичном круге на равном расстоянии друг от друга:

Все 9 комплексных корней степени 9 из единицы
Все 9 комплексных корней степени 9 из единицы

Первый корень \(w_1\) (точнее второй — единицу считаем нулевым корнем) называют образующим корнем степени \(n\) из единицы. Возведение его в нулевую, первую, вторую и так далее степени порождает последовательность нужных корней единицы, при этом на \(n\)-ном элементе последовательность зацикливается:

\[ w_n = e^{i \tau \frac{n}{n}} = e^{i \tau} = e^{i \cdot 0} = w_0 = 1 \]

Будем обозначать \(w_1\) как просто \(w\).

Упражнение. Докажите, что других корней быть не может.

Дискретное преобразование Фурье

Дискретным преобразованием Фурье называется вычисление значений многочлена в комплексных корнях из единицы:

\[ y_j = \sum_{k=0}^{n-1} x_n e^{i\tau \frac{kj}{n}} = \sum_{k=0}^{n-1} x_n w_1^{kj} \]

Обратным дискретным преобразованием Фурье назвается, как можно догадаться, обратная операция — интерполяция коэффициентов \(x_i\) по значениям \(X_i\).

\[ x_j = \frac{1}{n} \sum_{k=0}^{n-1} y_n e^{-i\tau \frac{kj}{n}} = \frac{1}{n} \sum_{k=0}^{n-1} y_n w_{n-1}^{kj} \]

Почему эта формула верна? При вычислении ПФ мы практически применяем матрицу к вектору:

\[ \begin{pmatrix} w^0 & w^0 & w^0 & w^0 & \dots & w^0 \\ w^0 & w^1 & w^2 & w^3 & \dots & w^{-1} \\ w^0 & w^2 & w^4 & w^6 & \dots & w^{-2} \\ w^0 & w^3 & w^6 & w^9 & \dots & w^{-3} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ w^0 & w^{-1} & w^{-2} & w^{-3} & \dots & w^1 \end{pmatrix}\begin{pmatrix} a_0 \\ a_1 \\ a_2 \\ a_3 \\ \vdots \\ a_{n-1} \end{pmatrix} = \begin{pmatrix} y_0 \\ y_1 \\ y_2 \\ y_3 \\ \vdots \\ y_{n-1} \end{pmatrix} \]

То есть преобразование Фурье — это просто линейная операция над вектором: \(W a = y\). Значит, обратное преобразование можно записать так: \(a = W^{-1}y\).

Как будет выглядеть эта \(W^{-1}\)? Автор не будет пытаться изображать логичный способ рассуждений о её получении и сразу её приведёт:

\[ W^{-1} = \dfrac 1 n\begin{pmatrix} w^0 & w^0 & w^0 & w^0 & \dots & w^0 \\ w^0 & w^{-1} & w^{-2} & w^{-3} & \dots & w^{1} \\ w^0 & w^{-2} & w^{-4} & w^{-6} & \dots & w^{2} \\ w^0 & w^{-3} & w^{-6} & w^{-9} & \dots & w^{3} \\ \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ w^0 & w^{1} & w^{2} & w^{3} & \dots & w^{-1} \end{pmatrix} \]

Проверим, что при перемножении \(W\) и \(W^{-1}\) действительно получается единичная матрица:

  1. Значение \(i\)-того диагонального элемента будет равно \(\frac{1}{n} \sum_k w^{ki} w^{-ki} = \frac{1}{n} n = 1\).

  2. Значение любого недиагонального (\(i \neq j\)) элемента \((i, j)\) будет равно \(\frac{1}{n} \sum_k w^{ik} w^{-jk} = \frac{1}{n} \sum_k w^k w^{i-j} = \frac{w^{i-j}}{n} \sum_k w^k = 0\), потому что все комплексные корни суммируются в ноль, то есть \(\sum w^k = 0\) (см. картинку — там всё симметрично).

Внимательный читатель заметит симметричность форм \(W\) и \(W^{-1}\), а также формул для прямого и обратного преобразования. На самом деле, эта симметрия нам сильно упростит жизнь: для обратного преобразования Фурье можно использовать тот же алгоритм, только вместо \(w^k\) использовать \(w^{-k}\), а в конце результат поделить на \(n\).

Зачем это надо?

Напомним, что мы изначально хотели перемножать многочлены следующим алгоритмом:

  1. Посчитаем значения в \(n+m\) каких-нибудь точках обоих многочленов

  2. Перемножим эти значения за \(O(n+m)\).

  3. Интерполяцией получим многочлен-произведение.

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

Соответствующий алгоритм называется быстрым преобразованием Фурье (англ. fast Fourier transform). Он использует парадигму «разделяй-и-властвуй» и работает за \(O(n \log n)\).

Схема Кули-Тьюки

Обычно, алгоритмы «разделяй-и-властвуй» делят задачу на две половины: на первые \(\frac{n}{2}\) элементов и вторые \(\frac{n}{2}\) элементов. Здесь мы поступим по-другому: поделим все элементы на чётные и нечётные.

Представим многочлен в виде \(P(x)=A(x^2)+xB(x^2)\), где \(A(x)\) состоит из коэффициентов при чётных степенях \(x\), а \(B(x)\) — из коэффициентов при нечётных.

Пусть \(n = 2k\). Тогда заметим, что

\[ w^{2t}=w^{2t \bmod 2k}=w^{2(t \bmod k)} \]

Зная это, исходную формулу для значения многочлена в точке \(w^t\) можно записать так:

\[ P(w^t)=A\left(w^{2(t\bmod k)}\right)+w^tB\left(w^{2(t\bmod k)}\right) \]

Сам алгоритм заключается в следующем: рекурсивно посчитаем БПФ для многочленов \(A\) и \(B\) и объединим ответы с помощью формулы выше. При этом в рекурсии нам нужно считать значения на корнях степени не \(n\), а \(k = \frac{n}{2}\), то есть на всех «чётных» корнях степени \(n\) (вида \(w^{2t}\)).

Заметим, что если \(w\) это образующий корень степени \(n = 2k\) из единицы, то \(w^2\) будет образующим корнем степени \(k\), то есть в рекурсию мы можем просто передать другое значение образующего корня.

Таким образом, мы свели преобразование размера \(n\) к двум преобразованиям размера \(\dfrac n 2\). Следовательно, общее время вычислений составит

\[ T(n)=2T\left(\dfrac n 2\right)+O(n)=O(n\log n) \]

Заметим, что предположение о делимости \(n\) на \(2\) имело существенную роль. Значит, \(n\) должно быть чётным на каждом уровне, кроме последнего, из чего следует, что \(n\) должно быть степенью двойки.

Реализация

Приведём код, считающий БПФ по схеме Кули-Тьюки:

typedef complex<double> ftype;
const double pi = acos(-1);

template<typename T>
vector<ftype> fft(vector<T> p, ftype w) {
    int n = p.size();
    if(n == 1)else { {
        return {p[0]};
    else {
        vector<T> AB[2];
        for(int i = 0; i < n; i++)
            AB[i % 2].push_back(p[i]);
        auto A = fft(AB[0], w * w);
        auto B = fft(AB[1], w * w);
        vector<ftype> res(n);
        ftype wt = 1;
        int k = n / 2;
        for(int i = 0; i < n; i++) {
            res[i] = A[i % k] + wt * B[i % k];
            wt *= w;
        }
        return res;
    }
}

vector<ftype> evaluate(vector<int> p) {
    while(__builtin_popcount(p.size()) != 1)
        p.push_back(0);
    return fft(p, polar(1., 2 * pi / p.size()));
}

Как обсуждалось выше, обратное преобразование Фурье удобно выразить через прямое:

vector<int> interpolate(vector<ftype> p) {
    int n = p.size();
    auto inv = fft(p, polar(1., -2 * pi / n));
    vector<int> res(n);
    for(int i = 0; i < n; i++)
        res[i] = round(real(inv[i]) / n);
    return res;
}

Теперь мы умеем перемножать два многочлена за \(O(n \log n)\):

vector<int> poly_multiply(vector<int> a, vector<int> b) {
    vector<int> A = fft(a);
    vector<int> B = fft(b);
    for (int i = 0; i < A.size(); i++)
        A[i] *= B[i];
    return interpolate(A);
}

Примечание. Приведённый выше код, являясь корректным и имея асимптотику \(O(n\log n)\), едва ли пригоден для использования на реальных контестах. Он имеет большую константу и далеко не так численно устойчивый, чем оптимальные варианты написания быстрого преобразования Фурье. Мы его приводим, потому что он относительно простой.

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

Здесь приведена одна из условно пригодных реализаций.

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

Number-theoretic transform

Нам от комплексных чисел на самом деле нужно было только одно свойство: что у единицы есть \(n\) комплексных корней. На самом деле, это не единственные алгебраические объекты, обладающие таким свойством.

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

Значит, для любого простого p в поле остатков от деления на него есть корень \(g\) степени \(p-1\) из единицы. Если при этом \((p-1)=c\cdot 2^k\), то \(g^{c}\) будет корнем степени \(2^k\), что позволяет применять метод Кули-Тьюки. Отсюда следует, что \(p=c\cdot2^k+1\). Практика показывает, что чисел такого вида очень много.

Применения

На самом деле, быстрое преобразование Фурье было придумано совсем не для работы с многочленами. Преобразование Фурье нужно само по себе.

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

Каждый раз, когда вы открываете джипег или проигрываете .mp3, где-то в компьютере считается БПФ.

Почти во всех алгоритмах обработки звука и эффективных алгоритмах обработки изображений оно используется.

Свёртки и корреляции

Пусть есть \(\{a_i\}_{i=0}^{n}\) и \(\{b_j\}_{j=0}^{m}\). Тогда свёрткой называют \(\{c_k\}_{k=0}^{m+n}\):

\[c_k=\sum\limits_{i=0}^{n} a_i b_{k-i}\]

Как мы видим, это просто \(k\)-ый коэффициент из произведения. Корреляцией же называют \(\{d_k\}_{-n}^{m}\):

\[d_k=\sum\limits_{i=0}^{n} a_i b_{k+i}\]

В обоих случаях мы предполагаем, что вне допустимых индексов последовательности равны нулям. Корреляцию можно интерпретировать двумя способами. С одной стороны, это коэффициент в произведении \(A(x) \cdot B(x^{-1})\), т.е. сдвинутая свёртка первой последовательности и развёрнутой второй. С другой стороны, \(d_k\) – в точности скалярное произведение последовательности \(a_i\) и отрезка последовательности \(b_j\), начинающегося в позиции \(k\). Именно свёртка и корреляция являются теми величинами, которые чаще всего нужно считать в задачах на преобразование Фурье.

Chirp Z-transform

Пусть нам дано некоторое число \(z\) и мы хотим вычислить значение многочлена в числах вида \(\{z^i\}_{i=0}^{n-1}\), т.е. множество чисел \(y_k=\sum\limits_{i=0}^{n-1} a_i z^{ik}\). Для этого сделаем замену \(ik=\dfrac{i^2+k^2-(i-k)^2}{2}\), после которой получим, что нам нужно вычислить

\[y_k=z^{\tfrac{k^2}{2}}\sum\limits_{i=0}^{n-1}\left(a_i z^{\tfrac{i^2}{2}}\right)z^{-\tfrac{(i-k)^2}{2}}\]

Что с точностью до множителя \(z^{\tfrac{k^2}{2}}\) является свёрткой двух последовательностей

\[u_i=a_i z^{\tfrac{i^2}{2}},~v_i=z^{-\tfrac{i^2}{2}}\]

Которая считается через произведение многочленов с такими коэффициентами. Но следует учесть, что здесь \(v_i\) определена также для отрицательных номеров. Данный метод среди прочего позволяет за \(O(n \log n)\) посчитать преобразование Фурье произвольной длины.

Одновременное преобразование вещественных многочленов

Пусть есть два многочлена

\[A(x)=\sum\limits_{i=0}^{n-1} a_i x^i,~B(x)=\sum\limits_{i=0}^{n-1} b_i x^i\]

С вещественными коэффициентами. Рассмотрим \(P(x)=A(x)+iB(x)\) и сопряжённый к нему.

\[\overline{P(w^k)}=A(\overline{w^k})-iB(\overline{w^k})=A(w^{n-k})-iB(w^{n-k})\]

Отсюда следует выражение для преобразования Фурье \(A(x)\) и \(B(x)\):

\[\begin{dcases} A(w^k)=\dfrac{P(w^k)+P(w^{n-k})}{2},\\ B(w^k)=\dfrac{P(w^k)-P(w^{n-k})}{2i} \end{dcases}\]

Одновременное преобразование можно произвести и в обратную сторону, рассматривая последовательность \(P(w^k)=A(w^k)+iB(w^k)\). После обратного преобразования мы получим \(P(x)=A(x)+iB(x)\).

Умножение по произвольному модулю

Нам нужно перемножить два многочлена, а затем вывести коэффициенты результата по модулю \(M\), не являющимся подходящим для быстрого преобразования Фурье. При этом достаточно большому, чтобы обычному преобразованию не хватало точности. Для разрешения данной ситуации представим многочлены в виде \[A(x)=A_1(x)+A_2(x)\cdot 2^k\]\[B(x)=B_1(x)+B_2(x) \cdot 2^k\]

где \(2^k \approx \sqrt M\). Тогда все коэффициенты будет \(O(\sqrt M)\), а произведение разложится как

\[A \cdot B = A_1 B_1 + (A_1 B_2+A_2 B_1) \cdot 2^k + A_2 B_2 \cdot 2^{2k}\]

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

Многомерное преобразование Фурье

Ранее мы работали с многочленами от одной переменной. Но аналогичные конструкции работают для многочлена от нескольких переменных. Считать значения многочлена теперь нужно в точках \((w_1^{k_1},w_2^{k_2}, \dots, w_m^{k_m})\). Оказывается, для такого преобразования достаточно поочерёдно сделать одномерное преобразование Фурье вдоль каждой координаты. В двумерном случае, например, нужно сначала сделать одномерное преобразование каждой строки, а затем каждого столбца.

Докажем это для двумерного случая. Мы хотим получить набор чисел

\[P_{uv}=P(w_1^u,w_2^v)=\sum_{i=0}^{n-1}\sum_{j=0}^{m-1} a_{ij} w_1^uw_2^v\]

Изначально мы имеем таблицу \(A_{uv}=a_{uv}\), после преобразования строк, мы получим

\[A'_{uv}=P_u(w_2^v)=\sum_{j=0}^{m-1}A_{uj}w_2^v=\sum_{j=0}^{m-1}a_{uj}w_2^v\]

После последующего преобразования столбцов же мы получим

\[A''_{uv}=P'_v(w_1^u)=\sum_{i=0}^{n-1} A'_{iv}w_1^u=\sum_{i=0}^{n-1}P_i(w_2^v) w_1^u=\sum_{i=0}^{n-1}\sum_{j=0}^{m-1} w_1^u w_2^v\]

Такое преобразование позволяет быстро вычислять двумерные свёртки \(C(x,y)=A(x,y)\cdot B(x,y)\) вида

\[c_{uv}=\mathop{\sum\sum}_{\substack{i_1+j_1=u\\i_2+j_2=v}} a_{i_1 i_2} b_{j_1 j_2}\]

Преобразование Уолша-Адамара и другие свёртки

Вычисляя значения многомерного многочлена в некоторых особых точках, мы можем научиться считать свёртки с другими условиями суммирования:

\[c_k=\sum\limits_{i | j=k} a_i b_j,~~~c_k=\sum\limits_{i \oplus j=k} a_i b_j,~~~c_k=\sum\limits_{i \& j=k} a_i b_j\]

Здесь \(|\), \(\&\) и \(\oplus\) соответствуют операциям побитового \(or\), \(and\) и \(xor\) соответственно.

  1. \(xor\). Рассмотрим значения многочлена в точках гиперкуба

    \(x \in \{-1,1\}^k\). Для таких точек верно соотношение

    \(x_i^a x_i^b=x_i^{a~xor~b}\), поэтому произведения значений

    многочленов в этих точках будут равны значениям многочлена, в

    котором мономы умножаются с учётом данного условия.

    Иначе говоря, если рассматривать степень \(x_i\) в мономе, как \(i\)-ый

    бит номера данного коэффициента, мы можем считать, что при

    произведении двух мономов мы получаем моном, чьему номеру

    соответствует \(xor\) номеров исходных мономов.

    Заметим, что такое вычисление есть ни что иное как вычисление

    многомерного преобразования Фурье в корнях степени \(2\) из единицы.

    Оно также называется преобразованием Уолша-Адамара. Здесь есть

    некоторое упрощение по сравнению с обычным преобразованием Фурье:

    во-первых, все вычисления можно производить в целых числах,

    во-вторых, \(w^{-1}=w=-1\), поэтому для обратного преобразования можно

    просто применить прямое и разделить всё на \(n\).

  2. \(or\). Теперь рассмотрим значения в точках \(x \in \{0,1\}^k\). Для них

    имеет место \(x_i^a x_i^b=x_i^{a~or~b}\), из чего следует, что при

    произведении мономов можно трактовать результат как моном с номером,

    равным побитовому \(or\) их номеров. Отдельно заметим, что посчитанное

    значение многочлена в точке это сумма его коэффициентов по всем

    подмаскам номера данной точки.

  3. \(and\). Чтобы посчитать свёртку по данной операции, нужно либо

    поменять все маски на их дополнения, посчитать свёртку по \(or\), а

    потом вернуться, либо воспользоваться идеей из прошлого пункта и

    провести суммирование по всем надмаскам. Это будет соответствовать

    значению многочлена в тех же точках, но с неявной перенумерацией,

    соответствующей переходу к дополнениям.

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

Метод Ньютона для функций над многочленами

Хотим решить уравнение \(f(x)=0\). \(f(x)\) можно представить в виде \(f(x)=f(x_0)+f'(x_0)\Delta x+O(\Delta x^2)\). Будем последовательно искать её нули, приближая линейной \(g(x_{n+1})=f(x_n)+f'(x_n)(x_{n+1}-x_n)\) на каждом шаге. Решая \(g(x_{n+1})=0\), приходим к \[x_{n+1}=x_n-\dfrac{f(x_n)}{f'(x_n)}\]

При этом \(f(x_{n+1})=O((x_{n+1}-x_n)^2)=O\left(\dfrac{f(x_n)^2}{f'(x_n)^2}\right)\). В случае обратимой производной это \(O(f(x_n)^2)\). Если \(x\) – многочлен, это значит, что используя метод Ньютона, мы будем на каждом шаге удваивать число точно известных коэффициентов. Наиболее распространённые функции от многочленов:

  1. Обратный ряд. Надо решить \(PQ=1 \Rightarrow f(P)=Q-P^{-1}\) и

    \(P_{n+1}=P_n-\dfrac{Q-P_n^{-1}}{P_n^{-2}}=P_n(2-QP_n)\).

  2. Экспонента. \(Q = \ln P \Rightarrow f(P)=Q-\ln P\) и

    \(P_{n+1}=P_n-\dfrac{Q-\ln P_n}{-P_n^{-1}}=P_n(1+Q-\ln P_n)\).

  3. Корень. \(Q=P^k \Rightarrow f(P)=Q-P^k\) и

    \(P_{n+1}=P_n+\dfrac{Q-P_n^k}{kP_n^{k-1}}=P_n\left(\dfrac{k-1}{k}+\dfrac{Q}{kP_n^k}\right)\).

В выражении для экспоненты есть логарифм, для его вычисления следует воспользоваться тем, что \(\left(\log P\right)'=P'P^{-1}\), что позволит восстановить коэффициенты при положительных степенях, а коэффициент при нулевой степени можно посчитать встроенными методами. Отметим, что все указанные алгоритмы работают за \(O(n \log n)\).

Разделяй и властвуй

Дано уравнение \(AX=B\) и мы хотим посчитать первые \(n\) коэффициентов \(X\). При этом \(B\) нам дан не весь, а подаётся по мере того, как мы узнаём коэффициенты \(X\) (например, \(B\) может зависеть от \(X\), см. пример ниже). Пусть мы нашли первые \(m \approx n/2\) коэффициентов и хотим найти следующие \(m\), тогда

\[\begin{cases} X=X_1+x^m X_2,\\ A=A_1+x^m A_2,\\ B=B_1+x^m B_2,\\ AX_1=B_1+x^m B_2' \end{cases} \implies AX\bmod x^{2m} = B_1+x^m(A_1X_2+B_2')=B_1+x^m B_2\] \[A_1X_2=B_2-B_2'\mod x^{m}\]

Значит, если запомнить \(B_2'\), который потом нужно будет вычесть из \(B_2\), можно свести вычисление оставшихся коэффициентов к этой же задаче размера \(\approx n/2\). Что даст алгоритм за \(O(n \log^2 n)\).

В качестве примера, этим методом также можно считать экспоненту, так как \[e^A=X \implies X'=A'X\]

На практике это проще, чем по схеме Ньютона и скрытая константа должна быть меньше, т.к. не нужно вызывать такие затратные функции как логарифм, который требует подсчёта обратного ряда. Заметим, что аналогичная схема может быть применена для подсчёта \(X=AB\), где \(A\) задан заранее, а \(B\) зависит от \(X\).

\[AB_1 = X_1+x^m X_2'\implies AB \bmod x^{2m}=X_1+x^m(A_1B_2+X_2')\]

\[X_2 = A_1 B_2 + X_2' \mod x^{m}\]

Деление и интерполяция

В завершение научимся делить многочлены с остатком, а также делать то, с чего всё началось – интерполировать многочлен и вычислять его на произвольных точках.

  1. Деление с остатком. Нам нужно представить

    \(A(x)=B(x)D(x)+R(x),~\deg R(x) < \deg B(x)\). Пусть

    \(\deg A = n,~\deg B=m\). Тогда \(\deg D = n-m\). При этом с учётом

    \(\deg R < m\) приходим к выводу, что коэффициенты при

    \(\{x^k\}_{k=m}^n\) не зависят от \(R(x)\). Получается, мы имеем систему

    из \(n-m+1\) линейных уравнений на \(n-m+1\) неизвестных (коэффициенты

    \(D\)).

    Рассмотрим

    \(A^r(x)=x^nA\left(x^{-1}\right),~B^r(x)=x^mB\left(x^{-1}\right),~~D^r(x)=x^{n-m}D\left(x^{-1}\right)\)

    – многочлены, в которых коэффициенты идут в обратном порядке. Для

    \(n-m+1\) старших коэффициентов исходных многочленов с учётом того,

    что \(P(x) \mod z^k\) – первые \(k\) коэффициентов \(P(x)\) имеем систему

    \[A^r(x)= B^r(x) D^r(x)\mod z^{n-m+1}\]

    Её решением будет \(D^r(x)=A^r(x)[B^r(x)]^{-1}\bmod z^{n-m+1}\), что

    позволяет найти \(D(x)\) и из него \(R(x)\).

  2. Многоточечное вычисление. Нужно вычислить \(P(x_i)\) для

    \(\{x_i\}_{i=1}^n\). Учитывая \(P(x_i)=P \mod (x - x_i)\), вычислим

    \(P \mod \prod\limits_{i=1}^{n/2-1}(x-x_i)\) и

    \(P \mod \prod\limits_{i=n/2}^{n} (x-x_i)\) и запустимся рекурсивно.

    Получим \(O(n \log^2 n)\).

  3. Интерполяция. Дан набор \(\{(x_i, y_i)\}_{i=0}^{n-1}\), нужно найти

    \(P: P(x_i)=y_i\). Пусть мы нашли многочлен \(P_1\) для первых \(n/2\)

    точек. Тогда \(P=P_1+P_2 \prod\limits_{i=0}^{n/2-1}(x-x_i)=P_1+P_2Q\).

    Нахождение \(P_2\) сведём к интерполяции и многоточечному вычислению:

    \(P_2(x_i)=\dfrac{y_i-P_1(x_i)}{Q(x_i)}\) для \(i>n/2\). Получим

    \(O(n \log^3 n)\).

Упражнения

Рюкзак

Есть \(n\) типов предметов. Предмет \(i\)-го типа имеет стоимость \(s_i\). Пусть \(s=\sum\limits_{i=1}^n s_i\). Предложите алгоритм, который для каждого \(w\leq s\) находит число способов выбрать подмножество предметов ровно с таким весом за \(O(s \log s \log n)\).

Степенной ряд

Даны числа \(k\) и \(n\). Найдите \(\sum\limits_{m=0}^n m^k\) за \(O(k \log k)\).

Общая схема Кули-Тьюки

Пусть \(n=pq\). Придумайте алгоритм, сводящий преобразование Фурье размера \(n\) к \(p\) преобразованиям Фурье размера \(q\) за \(O(n)\) дополнительных операций.

Арифметические прогрессии

Дано множество из \(n\) чисел от \(0\) до \(m\). Найдите число арифметических прогрессий длины \(3\) в этом множестве за \(O(m \log m)\).

Расстояние между точками

Даны \(n\) точек в прямоугольнике \(A \times B\). Для каждой возможной пары \((\Delta x,\Delta y)\) посчитайте сколько есть пар точек, таких что разность по \(x\)-координате между ними равна \(\Delta x\), а по \(y\)-координате соответственно \(\Delta y\) за \(O(AB \log AB)\).

Сопоставление шаблонов

Даны две строки \(s\) и \(t\). В них могут встречаться символы из множества \(\Sigma\), а также знаки вопроса. Найдите все позиции \(i\) такие, что если приложить строку \(t\) к строке \(s\) начиная с \(i\), то в любой позиции соответствующие символы в \(s\) и \(t\) должны либо совпадать, либо хотя бы один из них должен быть знаком вопроса за \(O(\Sigma n \log n)\) и за \(O(n \log n)\)*.

Линейные рекурренты*

Последовательность \(F_n\) задана как \(F_n = \sum\limits_{i=1}^k a_{k-i} F_{n-i}\). Даны коэффициенты \(\{a_i\}_{i=0}^{k-1}\) и начальные величины \(\{F_i\}_{i=0}^{k-1}\). Предложите алгоритм, вычисляющий \(F_n\) за \(O(k \log k \log n)\).

Степень многочлена*

Дан \(P(x)=\sum\limits_{i=0}^n a_i x^i\). Нужно найти первые \(n\) коэффициентов \(P^k(x)\) за \(O(n \log n)\).