Разделяй и властвуй. Повышение эффективности алгоритмов. Часть 3

В прошлой части мы рассмотрели общий подход к расчету эффективности алгоритмов с принципом "разделяй и властвуй", а также применение принципа к различным базовым алгоритмам.

Сегодня поговорим о следующем приеме. Как известно, составная часть принципа, это поделить задачу на подзадачи. Мы ни разу не касались, как именно делить. Просто делили на равные части. Но тут вот и есть нюанс, если поделить не абы как, а используя какую-то стратегию, то можно добиться применения принципа там, где это даже не очевидно. И именно эта тема станет предметом данной статьи на примере задачи умножения полиномов. Как и в предыдущих частях, я упрощаю математическую часть, опуская различные нюансы и частные случаи, с целью сохранить научно-популярный характер публикации. При этом я пытаюсь сохранить основные элементы алгоритмов и математических основ. Я хочу подать информацию в кратком доступном виде, в виде математического пересказа, и если у кого-либо возникнет интерес, тот может легко найти полные и строгие математические выкладки по данной теме.

Сигналы

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

Сигнал в общем случае, это некоторый показатель изменяющийся во времени. Это может быть огибающая звука, диаграмма распределения, статистические данные, яркость видеосигнала и др. Для того чтобы обрабатывать сигнал компьютерными алгоритмами сигнал "оцифровывают", т.е. заменяют непрерывный сигнал a(t) дискретной величиной, в разные моменты времени.

$a(t); t\in[0\dots T]\rightarrow \{a(t_0),a(t_1),\dots a(t_k)\};t_0,t_1,\dots t_k\in[0\dots T]$

Также принято следующее компактное представление оцифрованного сигнала, облегчающего математику с сигналом

$a(t)=\sum_{i=0}^{T-1}a(i)\delta (t-i) \quad(1)$

Данный сигнал оцифрован в набор из дискретных значений. Функция $\delta (t)$ это, так называемая дельта-функция, функция, которая везде равна 0, кроме точки 0, где она равна 1.
Далее компьютерная программа занимается преобразованием F(a) данного входного сигнала. И получает в результате некоторый выходной сигнал c(t)

Большой интерес в работе с сигналами представляют, так называемые, линейные преобразования. Это такие преобразования c(t)=[F(a)](t) исходного сигнала a(t) , что обладают следующими свойствами

не зависят от сдвига по времени. То есть, если сигнал $a'(t)=a(t+\lambda )$ сдвинут по времени, то результат его преобразования также сдвигается по времени на ту же величину $[F(a')](t)=c(t+\lambda )$
сохраняют линейное преобразование сигнала. То есть если сигнал $a'(t)=\lambda a(t)$ , то результат его преобразования будет $[F(a')](t)=\lambda c(t)$

Используя второе свойство линейности, можем записать преобразование используя представление сигнала из формулы (1)

$c(k)=[F(a)](k)=[F(\sum_{i=0}^{T-1}a(i)\delta(t-i))](k)=\sum_{i=0}^{T-1}a(i)[F(\delta(t-i))](k)=$

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

$=\sum_{i=0}^{T-1}a(i)[F(\delta(t))](k-i) \quad(2)$

И получается, чтобы выполнить линейное преобразование любого сигнала,нам нужно лишь посчитать одно преобразование от "дельта-сигнала" $F(\delta )$ , а преобразования любых других сигналов рассчитываются операциями сумм и произведения. Допустим, мы посчитали это одно преобразование $b(t)=F(\delta )$ и получили в результате сигнал b(t) . Давайте этот полученный сигнал, также разложим по формуле (1)

$b(t)=[F(\delta)](t)=\sum_{j=0}^{T-1}b(j)\delta(t-j)$

Тогда, продолжим преобразование в формуле (2)

$bc(k)=\sum_{i=0}^{T-1}a(i)[F(\delta)](k-i)=\sum_{i=0}^{T-1}a(i)\sum_{j=0}^{T-1}b(j)\delta(k-i-j)=\sum_{i=0}^{k}a(i)b(k-i)\quad (3)$

Наивное вычисление по этой формуле всех составляющих сигнала c(k) имеет сложность O(n^2) , но метод "Разделяй и властвуй" позволил получить сложность $O(n\cdot \log n)$ и это было очень существенное изменение. Это было революционное усовершенствование эффективности обработки сигналов. Благодаря этому изменению у нас работает цифровое видео, сотовая связь и другие вещи, без которых не представить жизнь в современности.
Резюме:

Любой сигнал для обработки оцифровывается и представляется в виде некоторого набора дискретных значений $a(i); \quad 0\leq i\leq T$
Линейное преобразование сигнала также возможно задать набором дискретных значений
Формула (3) позволяет вычислить результат линейного преобразования сигнала

Модель полиномов

Итак выше мы получили задачу, которую нам надо оптимально вычислять, чтобы обрабатывать сигналы. Сформулируем ее еще раз:
Даны k+1 значения a_i , b_i , где $0\leq i\leq k$ . Расширим эти наборы значений до размера 2k+1 , считая a_i=b_i=0 , если $2k\geq i \gt k$ . Требуется вычислить 2k+1 значений

$c_j=\sum_{i=0}^ja_ib_{j-i}\quad 0≤j≤2k \quad (4)$

Эффективное решение данной задачи позволит нам эффективно обрабатывать оцифрованные сигналы, как мы видели в разделе выше.
Как будем решать такую задачу? Более того, как свести эту задачу к методу "Разделяй и властвуй"? Если подойти к решению наивным методом, то мы для каждого c_j будем выполнять операций умножения. И всего у нас 2k+1 штук таких c_j , каждый из которых требует вычисления по формуле (4). В итоге, наивное вычисление дает эффективность вида O(n^2) и надо что-то улучшить, и где то надо "разделить и властвовать", чтобы улучшить эту эффективность.
Формула (4) выглядит довольно простой. Применима она не только для линейного преобразования сигналов. Можно легко вспомнить еще одну простую математическую задачу, для решения которой потребуется выполнить вычисления по формуле (4)
Пусть у нас есть два полинома:

$p_1=a_kx^k+a_{k-1}x^{k-1}+\dots +a_1x+a_0$ $p_2=b_kx^k+b_{k-1}x^{k-1}+\dots +b_1x+b_0$

И нам нужно вычислить новый полином, который является произведением этих двух данных полиномов

$p_1p_2=c_{2k}x^{2k}+c_{2k-1}x^{2k-1}+\dots c_1x+c_0=[a_kb_k]x^{2k}+[a_kb_{k-1}+a_{k-1}b_k]x^{2k-1}+$ $+\dots +[a_1b_0+a_0b_1]x+a_0b_0$

Легко проверяется, что в общем случае, чтобы вычислить c_j придется воспользоваться формулой

$c_j=\sum_{i=0}^ja_ib_{j-i}$

Что в точности соответствует нашей исходной задаче (4)
Таким образом мы свели исходную задачу к задаче нахождения произведения двух полиномов за оптимальное время, и теперь мы будем решать эту задачу.

Альтернативное представление полинома

Любой полином степени однозначно задается его коэффициентами a_j при степени x_j^j , где $0\leq j\leq k$ . Но есть еще один способ однозначно задать полином. Любой полином степени однозначно задается его k+1 значениями в любых различных точках $x_0,x_1,\dots x_k: A(x_0), A(x_1), \dots , A(x_k)$ Мы знаем, что прямая линия однозначно задается двумя различными точками, константа задается одной точкой. И это также верно для других степеней полиномов. Но что очень удобно, в применении к нашей задаче, если два перемножаемых полинома будут заданы не коэффициентами, а заданы в этой альтернативной "системе счисления" точками и значениями в одинаковых точках

$x_0, x_1, \dots x_k : A(x_0), A(x_1), \dots , A(x_k): B(x_0), B(x_1), \dots , B(x_k)$

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

$C(x_0)=A(x_0)\cdot B(x_0); C(x_1)=A(x_1)\cdot B(x_1);\dots C(x_{2k})=A(x_{2k})\cdot B(x_{2k})$

Что потребует всего 2k+1 операций умножения
Если мы найдем алгоритм, которым сможем эффективно переводить полиномы, задаваемые коэффициентами a_i ; b_i при степенях x_i^i в полиномы задаваемыми 2k+1 значениями в точках

$x_0, x_1, \dots x_{2k} : A(x_0), A(x_1), \dots ,A(x_{2k}); B(x_0), B(x_1), \dots , B(x_{2k})$

, то далее у нас будет задача линейной сложности по вычислению $C(x_0), C(x_1),\dots ,C(x_{2k})$ Под конец у нас останется еще одна задача, перевести полином, заданный 2k+1 значениями в точках $x_0,x_1,\dots x_{2k}$ : $C(x_0), C(x_1),\dots ,C(x_{2k})$ обратно в форму коэффициентов c_j при степени x_j^j , где $0\leq j\leq 2k$ . Тогда исходная задача (4) будет решена.
Как эффективно вычислить полином

$p=a_{2k}x^{2k-1}+a_{2k-1}x^{2k-2}+\dots +a_2x+a_1\quad (5)$

в точках? Допустим, мы за линейное время выберем случайные различные точек и будем наивно вычислять значения полинома в этих точках. Во первых, нам нужно вычислить все степени x^i , от 1 до 2k-1 , что потребует 2k-2 операций умножения. Потом нам еще потребуется 2k-1 операций умножения каждой степени на свой коэффициент a_i и, наконец, сложить слагаемых. Везде только линейные операции, но мы вычислим за линейное время только одну точку. А у нас таких точек и мы опять выходим на сложность O(n^2) , что перечеркивает все наши усилия.
Но, что если точки, в количестве , будут заданы не совсем произвольно? Например, это будет положительно-отрицательных пар

$\pm x_1, \pm x_2, \dots , \pm x_k$

Полином из (5) можно разрезать на две части, с четными степенями и с нечетными степенями .

$p=x[a_{2k}x^{2k-2}+a_{2k-2}x^{2k-4}+\dots a_2]+[a_{2k-1}x^{2k-2}+a_{2k-3}x^{2k-4}+\dots a_1]$

В этом случае данный полином можно записать как

$p=xP_{нечет}(x^2)+P_{чет}(x^2)\quad (6)$

Где два полинома $P_{нечет}$ и $P_{чет}$ имеют степени, как минимум в два раза меньше, чем в нашем исходном полиноме (5)
Но самое главное, что для четно-нечетных пар нам потребуется вычислить значения этих полиномов уже не в точках, а в точках, так как исходные точки - это четно-нечетные пары, а нам нужны лишь их квадраты, которые одинаковы для каждой четно-нечетной пары. После чего за линейное время мы сделаем вычисления значений по формуле (6) для нашего исходного полинома.
Ничего не напоминает? Мы начали вычислять значения полинома степени 2k-1 в точках, но свели эту задачу к вычислению двух полиномов степени k-1 , каждого в точках. Это же наш "разделяй и властвуй". Если мы рекурсивно продолжим и дальше делить степени полиномов, то на следующем уровне рекурсии у нас будет четыре полинома степени k/2-1 , которые нам надо будет посчитать в k/2 точках.
Рекурсивное соотношение

Согласно мастер-теореме дает результат $O(n\cdot \log n)$ . Ровно то, что доктор прописал. Гораздо эффективнее исходных O(n^2) , что давал исходный наивный алгоритм.
Но, позвольте, ведь, когда мы спускались на нижний уровень рекурсии, мы возвели в квадрат наши положительно-отрицательные пары и получили вместо положительно-отрицательных пар положительные числа, с которыми уже так просто не спустишься на следующий уровень рекурсии.
Казалось бы да, но на помощь приходят комплексные числа. Допустим у нас есть набор из 4-ёх положительно-отрицательных пар:

$\pm 1;\pm(\frac1{\sqrt 2}+\frac{i}{\sqrt 2}); \pm i; \pm(\frac1{\sqrt 2}-\frac{i}{\sqrt 2})$

После того как мы каждую из этих пар возведем в квадрат, исчезнет плюс/минус, и получится 4 числа.

Заметьте, у нас опять вышли две положительно отрицательные пары:

$\pm 1;\pm i;$

Если мы возведем их в квадрат, то опять исчезнет плюс/минус, первое число даст +1, второе даст -1 и у нас опять будет одна положительно отрицательная пара.
То есть, для вычисления значения полинома можно выбрать такие точки, которые будут положительно отрицательными парами, и при возведении в квадрат останутся также набором положительно-отрицательных пар, что позволит совершать рекурсии. Тогда значения полинома можно вычислить по формуле (6) для квадратов этих положительно отрицательных пар. Затем, так как эти квадраты все еще остались положительно-отрицательными парами, мы можем повторить рекурсивно разбиение по формуле (6) и выйти на полиномы еще меньшей степени.
В общем случае, чтобы выбрать нужные точки, в которых провести вычисления, в количестве штук нам нужно найти все решения уравнения:

$\omega^n= 1$

Тогда в качестве точек, образующих положительно-отрицательные пары будет для n=2 , пара (+1;-1) , Для n=4 пара (+1;-1;+i;-i) и так далее.
Уравнение в комплексных числах решается довольно просто и имеет следующие корни

$\Omega = \{1;e^{i2\pi/n};(e^{i2\pi/n})^2;(e^{i2\pi/n})^3\dots (e^{i2\pi/n})^{n-1}\} \quad (7)$

Для чётного у нас в этом множестве образуются необходимые положительно-отрицательные пары:

для 1 это $(e^{i2\pi/n})^{n/2}=-1$
для $e^{i2\pi/n}$ это $(e^{i2\pi/n})^{n/2+1}=-e^{i2\pi/n}$ и т.д.

При возведении этих чисел в квадрат, каждая положительно-отрицательная пара дает число, которое также является одним из корней, то есть также находится в множестве $\Omega$ . Количество чисел при возведении в квадрат сокращается в два раза. На предпоследней итерации алгоритма останутся два числа +1 и -1, и наконец одно число 1.
Подведем промежуточный итог:

Мы получили алгоритм, который позволяет вычислить значения полиномов степени в точке за время $O(n\cdot \log n)$ и, используя этот алгоритм, перевели перемножаемые полиномы, задаваемые коэффициентами ; при степенях в полиномы задаваемыми значениями в точках

$x_0, x_1, \dots x_{n} : A(x_0), A(x_1), \dots ,A(x_{n}); B(x_0), B(x_1), \dots , B(x_{n})$

Далее мы за линейное время вычислили произведения этих полиномов

$C(x_0)=A(x_0)\cdot B(x_0); C(x_1)=A(x_1)\cdot B(x_1);\dots C(x_{n})=A(x_{n})\cdot B(x_{n})$

Осталась одна задача. Перевести полином, заданный значениями в точке в привычный вид

$c_nx^n+c_{n-1}x^{n-1}+\dots +c_1x+c_0$

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

Это известная задача - интерполяция полинома степени n-1 по точкам.

В точках $\{x_1, x_2, \dots , x_n\}$ известны значение полинома $\{C(x_1), C(x_2), \dots , C(x_n)\}$

$C(x_1)=c_nx_1^{n-1}+c_{n-1}x_1^{n-2}+\dots +c_1\quad (8)$ $C(x_2)=c_nx_2^{n-1}+c_{n-2}x_2^{n-1}+\dots +c_1$ $\dots$ $C(x_n)=c_nx_n^{n-1}+c_{n-1}x_n^{n-2}+\dots +c_1$

По этим значениям надо вычислить коэффициенты $\{c_1, c_2, \dots , c_n\}$

В общем случае (8) это линейная система уравнений, которую надо решить относительно неизвестных $\{c_1, c_2, \dots , c_n\}$ Чтобы увидеть решение этой задачи в общем случае, запишем систему (8) в привычном матричном виде

$\left( \begin{array}{cccc} C(x_1)\\ C(x_2)\\ \vdots \\ C(x_n) \end{array} \right)=\left( \begin{array}{cccc} 1 & x_1 & \ldots & x_1^{n-1}\\ 1 & x_2 & \ldots & x_2^{n-1}\\ \vdots & \vdots & \ddots & \vdots\\ 1 & x_n & \ldots & x_n^{n-1} \end{array} \right)\left( \begin{array}{cccc} c_1\\ c_2\\ \vdots \\ c_n \end{array} \right)=M\left( \begin{array}{cccc} c_1\\ c_2\\ \vdots \\ c_n \end{array} \right)$

Надо просто найти матрицу, обратную к квадратной $n \times n$ матрице , и в этом случае получим интерполяцию, которая и будет решением нашей задачи

$\left( \begin{array}{cccc} c_1\\ c_2\\ \vdots \\ c_n \end{array} \right)=M^{-1}\left( \begin{array}{cccc} C(x_0)\\ C(x_1)\\ \vdots \\ C(x_{n-1}) \end{array} \right) \quad (9)$

Проблема в том, что в общем случае нахождение обратной матрицы имеет сложность O(n^3) , что нас не устраивает и нам необходимо как-то решить эту задачу более эффективным способом.

На помощь приходит то, что множество точек для интерполяции $\{x_1, x_2, \dots , x_n\}$ у нас не произвольное, а жестко задано в (7).

Тогда, с учетом (7) матрица M выглядит следующим образом

$M=\left( \begin{array}{cccc} 1 & 1 & \ldots & 1\\1 & \omega & \ldots & \omega^{n-1}\\ 1 & \omega^2 & \ldots & \omega^{2(n-1)}\\ \vdots & \vdots & \ddots & \vdots\\ 1 & \omega^{n-1} & \ldots & \omega^{(n-1)(n-1)} \end{array} \right)\quad, где\quad \omega=e^{i2\pi/n}$

Если задать формулу для $M=[m_{ij}]$ элемента этой матрицы (нумеруя индексы для упрощения формулы с 0 до n-1 ), то это будет

$m_{ij}=\omega^{i \cdot j}$

Данная матрица называется оператором дискретного преобразования Фурье и обладает рядом замечательных свойств, применяемые в различных областях. Заметим, что эта матрица является постоянной, совершенно не зависящей от того, какие полиномы мы перемножаем и зависящей только от размерности . Чтобы получить эту матрицу для заданной размерности надо вычислить $\omega = e^{i2\pi/n}$ , после чего заполнить матрицу $n\times n$ по схеме выше, где $m_{ij}=\omega^{i \cdot j}$ . Вычислительная система может вычислить заранее эту матрицу и держать ее в памяти, чтобы дальше заниматься только умножением полиномов степени и не тратить время на вычисление элементов матрицы. Обозначим матрицу дискретного преобразования Фурье размерности n, как $DFT^{(n)}$

Возьмем матрицу $M^*=[m_{ij}^*]$ , где каждый элемент матрицы получен из исходной матрицы путем комплексного сопряжения каждого элемента. И давайте вычислим, чему равняется матрица $X=[x_{ij}]=M\times M^*$ равная произведению исходной матрицы, и матрицы где все элементы получены из исходной путем сопряжения.

По известной формуле умножения двух матриц

$x_{ij}=\sum_{k=0}^{n-1}m_{ik}m^*_{kj}=\sum_{k=0}^{n-1}\omega^{i \cdot k}\cdot(\omega^{k \cdot j})^*$

Заметим, что для сопряженного значения справедливо следующее:

$(\omega^k)^*=(e^{i2\pi k/n})^*=e^{-i2\pi k/n}=\omega^{-k}\quad (10)$

Тогда

$x_{ij}=\sum_{k=0}^{n-1}\omega^{i \cdot k}\cdot(\omega^{k \cdot j})^*=\sum_{k=0}^{n-1}\omega^{(i-j) \cdot k}=n\delta_{ij}$

Действительно, если i=j то под суммой будут единицы. Если $i\neq j$ , то под сумой будет набор комплексных чисел, которые в своем множестве на комплексной плоскости похожи на спицы колеса, эти комплексные вектора смотрят в разные стороны и "гасят" друг друга.

Интуитивно понятно, но если хотите более математического обоснования, то оно следующее:
Пусть $\sum_{k=0}^{n-1}\omega^{m \cdot k}=C$ для некоторого $\omega^m\neq1$ , $\omega^n=1$
Умножим эту сумму на $\omega^m$ , тогда

$\omega^mC=\sum_{k=0}^{n-1}\omega^{m \cdot (k+1)}=\omega^n+\sum_{k=1}^{n-1}\omega^{m \cdot k}=1+\sum_{k=1}^{n-1}\omega^{m \cdot k}=\sum_{k=0}^{n-1}\omega^{m \cdot k}=C$

Так как $\omega^m\neq1$ , значит такое равенство возможно только, если C=0

Следовательно, матрица которую мы вычислили, это просто единичная матрица умноженная на . Получается, что $M\times M^*=nI$ . А значит, матрица $(1/n)M^*=[m_{ij}^*/n]$ это и будет та самая обратная матрица, которую нам надо найти. Используя (10) мы легко можем выписать, как выглядит эта матрица

$M^{-1}=\left( \begin{array}{cccc} 1 & 1 & \ldots & 1\\1 & \omega^{-1} & \ldots & \omega^{-(n-1)}\\ 1 & \omega^{-2} & \ldots & \omega^{-2(n-1)}\\ \vdots & \vdots & \ddots & \vdots\\ 1 & \omega^{-(n-1)} & \ldots & \omega^{-(n-1)(n-1)} \end{array} \right)\quad, где\quad \omega=e^{i2\pi/n}$

То есть $M^{-1}=[m_{ij}^*]=\omega ^{-ij}$ (10)
Заметим, что имея вычисленную матрицу $DFT^{(n)}$ нет необходимости вычислять элементы обратной матрицы, достаточно отразить матрицу по вертикали. Действительно, учитывая, что $\omega ^n=1$ :

$\omega ^{-ij}=\omega^{in}\times \omega ^{-ij}=\omega ^{i(n-j)}$

Аналогично, можно отразить элементы и по горизонтали. Обозначим эту матрицу как $TFD^{(n)}$ , то есть по смысле это обратная матрица к $DFT^{(n)}$ , или просто отраженная по вертикали, как мы видели выше, что я отразил, переставив буквы DFT в TFD.
Чтобы наивно провести вычисления по формуле (9) потребуется вычисления сложности O(n^2) , коэффициентов c_i , по каждому надо провести операций умножения и сложения. Но не волнуйтесь, немного терпения, мы почти на финише. Осталось еще раз "разделить и повластвовать", чтобы эффективно провести вычисления по формуле (9)
Прием следующий. Давайте в матрице (10) переставим колонки. Сначала будут идти все четные колонки (0,2,4...), потом все нечетные (1,3,5...)
Чтобы результат вычислений по формуле (9) после такой перестановки не изменился, то мы одновременно переставим коэффициенты C() , сверху будут все те, что были с четными индексами (C(x_0), C(x_2), C(x_4), ...), потом те, что с нечетными индексами (C(x_1), C(x_3), C(x_4), ...). В итоге формула (9) после раскрытия станет выглядеть так:

$\left( \begin{array}{cccc} c_1\\ c_2\\ c_3\\ \vdots \\c_n \end{array} \right)=\left( \begin{array}{cccc}1 & 1 & 1 & \ldots & 1 & 1 & \vdots \\1 & \omega^{-2} & \omega^{-4} & \ldots & \omega^{-1} & \omega^{-3} & \vdots \\1 & \omega^{-4} & \omega^{-8} & \ldots & \omega^{-2} & \omega^{-6} & \vdots \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \ddots \\1 & \omega^{-2(n-2)} & \omega^{-4(n-2)} & \ldots & \omega^{-(n-2)} & \omega^{-3(n-2)} & \vdots \\1 & \omega^{-2(n-1)} & \omega^{-4(n-1)} & \ldots & \omega^{-(n-1)} & \omega^{-3(n-1)} & \vdots \end{array} \right)\left( \begin{array}{cccc} C(x_0)\\ C(x_2)\\ C(x_4)\\ \vdots \\C(x_1)\\ C(x_3) \\ \vdots \end{array} \right)\quad, где\quad \omega=e^{i2\pi/n}$

Разобьем эту матрицу на 4 равных матричных блока A,B,C,D , размером $(n/2\times n/2)$ каждый

$\left( \begin{array}{cccc} c_1\\ c_2\\ c_3\\ \vdots \\c_n \end{array} \right)=\left( \begin{array}{cccc} A & B \\ C & D\end{array} \right)\left( \begin{array}{cccc} C(x_0)\\ C(x_2)\\ C(x_4)\\ \vdots \\C(x_1)\\ C(x_3) \\ \vdots \end{array} \right)$

Элементы матрицы $A=[a_{ij}]=[(\omega^2)^{-ij}]=[TFD^{n/2}]$
Элементы матрицы $C=[c_{ij}]=[(\omega^2)^{-(i-n/2)j}]=[\omega^n(\omega^2)^{-ij}]=[(\omega^2)^{-ij}]=[TFD^{n/2}]$
Элементы матрицы $B=[b_{ij}]=[(\omega)^{-j(2i+1)}]=[\omega^{-j}(\omega^2)^{-ij}]=[\omega^{-j}TFD^{n/2}]$
Элементы матрицы $D=[d_{ij}]=[(\omega)^{-(j+n/2)(2i+1)}]=[\omega^{-j}(\omega)^{ni}(\omega^{-n/2}]=[-\omega^{-j}(\omega^2)^{-ij}]=[-\omega^{-j}TFD^{n/2}]$

Получается, что вместо наивного вычисления по формуле (9), мы можем провести следующие рекурсивные вычисления:
Для j<n/2

$c_j=A \left( \begin{array}{cccc} C(x_0)\\ C(x_2)\\ C(x_4)\\ \vdots \\ C(x_n) \end{array} \right)+B \left( \begin{array}{cccc} C(x_1)\\ C(x_3)\\ C(x_5)\\ \vdots \\ C(x_n-1) \end{array} \right)=TFD^{(n/2)} \left( \begin{array}{cccc} C(x_0)\\ C(x_2)\\ C(x_4)\\ \vdots \\ C(x_n) \end{array} \right)+\omega^{-j}TFD^{n/2} \left( \begin{array}{cccc} C(x_1)\\ C(x_3)\\ C(x_5)\\ \vdots \\ C(x_n-1) \end{array} \right)$

Для $j\geq n/2$

$c_j=C \left( \begin{array}{cccc} C(x_0)\\ C(x_2)\\ C(x_4)\\ \vdots \\ C(x_n) \end{array} \right)+D \left( \begin{array}{cccc} C(x_1)\\ C(x_3)\\ C(x_5)\\ \vdots \\ C(x_n-1) \end{array} \right)=TFD^{(n/2)} \left( \begin{array}{cccc} C(x_0)\\ C(x_2)\\ C(x_4)\\ \vdots \\ C(x_n) \end{array} \right)-\omega^{-j}TFD^{n/2} \left( \begin{array}{cccc} C(x_1)\\ C(x_3)\\ C(x_5)\\ \vdots \\ C(x_n-1) \end{array} \right)$

И мы опять свели задачу вычисления формулы (9) размера к двум подзадачам размерности n/2 . Таким образом, мы опять приходим к рекурсивному соотношению:

что по мастер-теореме дает нам результат $O(n\times log(n))$

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

Заключение

Мы рассмотрели решение вычислительной задачи

$c_j=\sum_{i=0}^ja_ib_{j-i}\quad 0≤j≤2k \quad (4)$

Мы разделили ее на следующие этапы, где почти на каждом этапе применяли метод "Разделяй и властвуй"

Свели задачу к умножению двух полиномов
Привели полиномы к альтернативному представлению за время $O(n\times log(n))$
За линейное время нашли произведениб полиномов этом альтернативном представлении
Привели результирующий полином из альтернативного к традиционному представлению за время $O(n\times log(n))$

Источник: https://habr.com/ru/articles/745360/

Вернуться к списку

Интересные статьи

Оказание первой помощи. Практическая часть

Судя по предыдущей части, оказание первой помощи оказалась трепетной темой. Поэтому я решилась написать вторую часть, посвященную практике. Можете сохранить эту статью себе в заметки, в экстренно...

Ностальгические игры: Fallout 2 (часть 1)

В 1997 году в свет вышла Fallout – игра, ставшая эталоном жанра. Она покорила умы и сердца миллионов геймеров по всему миру, которые не спали ночами, говоря себе «Ну еще один квест и точно на боковую!...

Безопасное хранение данных IoT в частном блокчейне Ethereum. Часть 3

Основная задача:Изучить, как хранить данные IoT на комбинации on-chain (Ethereum Blockchain) и off-chain хранилищ (IPFS и Ethereum Swarm) в зашифрованном виде и использовать их в модели публикации-под...

Чистый код для TypeScript — Часть 2

Первая часть, судя по комментариям вызвала неоднозначное мнение, особенно что касалось части enum. Где-то я так же могу не соглашаться, как с автором оригинала, так и с некоторыми комментариями. ...

Часть 4: Всё-таки запускаем Linux на RocketChip RISC-V

На картинке Linux kernel шлёт вам привет через GPIO. В этой части истории с портированием RISC-V RocketChip на китайскую плату с Cyclone IV мы всё-таки запустим Linux, а также научимся сами кон...