УказательРазделыОбозначенияАвторО проекте


Статус материалов: черновик.


Умножение полиномов

Идея

§

Для понимания материалов настоящего пункта крайне желательно просмотреть начальные пункты раздела ИНТЕРПОЛЯЦИЯ.

Задача. Заданы два полинома в канонической форме:

f(x)=a_0+a_1x+a_2x^2+\dots+a_{n}x^n,\ g(x)= b_0+b_1x+b_2x^2+\dots+b_{n}x^n ;

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

h(x)=f(x)g(x)=c_0+c_1x+\dots+c_{2n}x^{2n} \ .

Предположим сначала, что коэффициенты полинома и известны приближенно (т.е. не являются целыми числами). Нам нужно организовать нахождение набора коэффициентов (c_0,c_1,\dots,c_{2n}), (также приближенно, с известной точностью). Предлагается следующий вариант решения этой задачи, альтернативный традиционной схеме умножения «в столбик».

Вычисляем две таблицы значений для полиномов f_{}(x) и g_{}(x) при 2n+1 различных значениях переменной:

\begin{array}{c|ccccc} x & x_1 & x_2 & \dots & x_{2n+1} \\ \hline f(x) & f(x_1) & f(x_2) &\dots & f(x_{2n+1}) \end{array} \qquad u \qquad \begin{array}{c|ccccc} x & x_1 & x_2 & \dots & x_{2n+1} \\ \hline g(x) & g(x_1) & g(x_2) &\dots & g(x_{2n+1}) \end{array}

Очевидно, если h(x)=f(x)g(x), то h(x_j)=f(x_j)g(x_j) при j\in \{1,\dots,2n+1\}. Таким образом, для неизвестного полинома h(x) мы строим интерполяционную таблицу

\begin{array}{c|ccccc} x & x_1 & x_2 & \dots & x_{2n+1} \\ \hline h(x) & f(x_1)g(x_1) & f(x_2)g(x_2) &\dots & f(x_{2n+1}) g(x_{2n+1}) \end{array}

из которой вполне определяется единственный полином степени \le 2n. Теперь задача свелась к извлечению из этой таблицы коэффициентов полинома c_{j}. Ну, это можно организовать разными способами, например, изложенным ☞ЗДЕСЬ. Формализуем предложенный алгоритм:

В чем же преимущество этого метода перед традиционным — «в столбик»? — Да пока что ни в чем. А преимущества начинаются когда мы специально подберем значения переменной x_{}, т.е. узлы интерполяции. Ведь их-то можно выбирать произвольными — воспользуемся же этим произволом!

Умножение полиномов по методу Карацубы

§

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

Рассмотрим два линейных полинома f(x)=a_0+a_1x и g(x)= b_0+b_1x. Их умножение по стандартной схеме, т.е. вычисление коэффициентов произведения

h(x)= a_0 b_0+(a_1b_0+a_0b_1)x+ a_1 b_1 x^2

«стоит» 4_{} операции умножения и одну операцию сложения. Попробуем организовать вычисление этого же произведения по иной схеме — имея целью сэкономить на операции умножения. Для этого сформулируем задачу интерполяции: найти линейный полином L(x) такой, что

h(x)\equiv L(x) + a_1 b_1 x^2 \quad \iff \quad L(x)\equiv f(x)g(x)- a_1 b_1 x^2 \ .

Узлами интерполяции возьмем \{ 0,\, 1 \}; интерполяционная таблица —

\begin{array}{c|cc} x & 0 & 1 \\ \hline L(x) & a_0b_0 & (a_0+a_1)(b_0+b_1) - a_1 b_1 \end{array}

Интерполяционный полином в записываем в форме Лагранжа:

L(x)=-a_0b_0 (x-1)+[ (a_0+a_1)(b_0+b_1) - a_1 b_1 ]x \equiv a_0b_0+[ (a_0+a_1)(b_0+b_1)-a_0b_0 - a_1 b_1 ]x \ .

В результате схема вычисления произведения f(x)g(x) становится следующей:

h(x)=a_0b_0+[ (a_0+a_1)(b_0+b_1)-a_0b_0 - a_1 b_1 ]x+ a_1 b_1 x^2 \ ,

и эта схема «стоит» 3_{} операции умножения и 4_{} операции сложения. В программировании операция умножения «стоит» существенно дороже операции сложения, так что предлагаемая схема умножения — известная как метод Карацубы — реализуется на процессорном уровне.

Она допускает очевидное развитие. Так, к примеру, умножение полиномов степени 7_{} можно свести к умножению полиномов степени 3_{}:

f(x)\equiv\overbrace{(a_0+a_1x+a_2x^2+a_3x^3)}^{f_0(x)}+x^4\overbrace{(a_4+a_5x+a_6x^2+a_7x^3)}^{f_1(x)},
g(x)\equiv \underbrace{(b_0+b_1x+b_2x^2+b_3x^3)}_{g_0(x)}+x^4\underbrace{(b_4+b_5x+b_6x^2+b_7x^3)}_{g_1(x)},
\Rightarrow h(x)=f_0g_0+[ (f_0+f_1)(g_0+g_1)-f_0g_0 - f_1 g_1 ]x^4+ f_1 g_1 x^8 \ .

В свою очередь, каждый полином \{f_j(x),g_j(x)\}_{j=0,1} можно представить с помощью линейных полиномов f_0(x)\equiv a_0+a_1x+x(a_2+a_3x),\dots, что позволяет свести умножение к линейному случаю, уже рассмотренному выше.

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

Можно, например, выбирать узлы интерполяции симметричными относительно 0_{} — см. пример, разбираемый ЗДЕСЬ. Такой выбор позволит разбить задачу нахождения коэффициентов полинома h(x) отдельно для его четных коэффициентов и отдельно — для его нечетных.

А что если пойти дальше и рассмотреть «сверхсимметричное» множество — множество комплексных корней (2n+1)-й степени из единицы? На первый взгляд, это кажется неразумным — поскольку, как правило, перемножаются полиномы с вещественными коэффициентами, то и результат должен быть вещественным! Оказывается, однако, что такой выбор узлов интерполяции здорово упрощает вычисления. — Почему? — Например, потому, что матрицу Вандермонда для такого выбора узлов довольно просто обратить.

П

Пример. Перемножить полиномы f(x)=2+3\,x+5\,x^2 и g(x)=-1-2\,x+2\,x^2.

Решение. В данном примере n_{}=2. Рассмотрим множество корней 5_{}-й степени из 1_{}:

\varepsilon_k = \cos \frac{2 \pi k}{n} + \mathbf i \sin \frac{2 \pi k}{n} \quad npu \ k \in\{0,1,2,3,4,5\} \ .

Их значения мы берем узлами интерполяции — только перенумеруем их в обратном порядке:

x_1=\varepsilon_0=1,
x_2=\varepsilon_{-1}= \frac{1}{4} \left( \scriptstyle{(\sqrt{5}-1)} - \displaystyle{\mathbf i} \scriptstyle{\sqrt{2 (5+\sqrt{5})}} \right),\ x_3=\varepsilon_{-2}=\frac{1}{4} \left( -\scriptstyle{(\sqrt{5}+1)} - \displaystyle{\mathbf i} \scriptstyle{\sqrt{2 (5-\sqrt{5})}} \right),
x_4=\varepsilon_{-3}=\frac{1}{4} \left( -\scriptstyle{(\sqrt{5}+1)} +\displaystyle{\mathbf i} \scriptstyle{\sqrt{2 (5-\sqrt{5})}} \right), \ x_5=\varepsilon_{-4}=\frac{1}{4} \left( \scriptstyle{(\sqrt{5}-1)} +\displaystyle{\mathbf i} \scriptstyle{\sqrt{2 (5+\sqrt{5})}} \right) \ .

Действуя по предложенной в предыдущем пункте схеме, вычисляем значения интерполяционных таблиц.

\begin{array}{c|c|c|c|c|c} k & 0 & -1 & -2 & -3 & - 4 \\ \hline \varepsilon_k & 1 & \frac{1}{4} \left( \scriptstyle{(\sqrt{5}-1)} - \displaystyle{\mathbf i} \scriptstyle{\sqrt{2 (5+\sqrt{5})}} \right) & \frac{1}{4} \left( \scriptstyle{(-\sqrt{5}-1)} - \displaystyle{\mathbf i} \scriptstyle{\sqrt{2 (5-\sqrt{5})}} \right) & \\ \hline f(\varepsilon_k) & 10 & -\scriptstyle{\frac{\sqrt{5}}{2}}-\displaystyle{\mathbf i}\scriptstyle{\frac{\sqrt{2}}{8}\sqrt{5+\sqrt{5}}(5\sqrt{5}+1)} & \scriptstyle{\frac{\sqrt{5}}{2}}+\displaystyle{\mathbf i}\scriptstyle{\frac{\sqrt{2}}{8}\sqrt{5+\sqrt{5}}(5\sqrt{5}-1)} & & \\ \hline g(\varepsilon_k) & -1 & \scriptstyle{-\sqrt{5}-1}+\displaystyle{\mathbf i} \scriptstyle{\frac{\sqrt{2}}{4}\sqrt{5+\sqrt{5}}(3-\sqrt{5})} & \scriptstyle{\sqrt{5}-1}+\displaystyle{\mathbf i} \scriptstyle{\frac{\sqrt{2}}{4}\sqrt{5-\sqrt{5}}(3+\sqrt{5})} & & \\ \hline h(\varepsilon_k)=f(\varepsilon_k)g(\varepsilon_k) & -10 & \scriptstyle{7\frac{\sqrt{5}}{2}}+ \displaystyle{\mathbf i}\scriptstyle{\frac{\sqrt{2}}{8}\sqrt{5+\sqrt{5}}(31+3\sqrt{5})} & -\scriptstyle{7\frac{\sqrt{5}}{2}}+ \displaystyle{\mathbf i}\scriptstyle{\frac{\sqrt{2}}{8}\sqrt{5-\sqrt{5}}(31-3\sqrt{5})} & & \\ \end{array}

(здесь элементы 5_{}-го столбца будут комплексно-сопряженными к элементам 3_{}-го столбца, а элементы 6_{}-го столбца — комплексно-сопряженными к элементам 4_{}-го столбца). Собственно вычисление значений полиномов можно переписать в матричном виде:

\left(\begin{array}{l} f(\varepsilon_0) \\ f(\varepsilon_{-1})\\ f(\varepsilon_{-2}) \\ f(\varepsilon_{-3}) \\ f(\varepsilon_{-4}) \end{array} \right) = F \left(\begin{array}{l} a_0 \\ a_1\\ a_2 \\ 0 \\ 0 \end{array} \right),\quad \left(\begin{array}{l} g(\varepsilon_0) \\ g(\varepsilon_{-1})\\ g(\varepsilon_{-2}) \\ g(\varepsilon_{-3}) \\ g(\varepsilon_{-4}) \end{array} \right) = F \left(\begin{array}{l} b_0 \\ b_1\\ b_2 \\ 0 \\ 0 \end{array} \right) \quad npu \ F=\left( \begin{array}{lllll} 1 & 1 & 1 & 1 & 1 \\ 1 & \varepsilon_{-1} & \varepsilon_{-1}^2 & \varepsilon_{-1}^3 & \varepsilon_{-1}^{4} \\ 1 & \varepsilon_{-2} & \varepsilon_{-2}^2 & \varepsilon_{-2}^3 & \varepsilon_{-2}^{4} \\ 1 & \varepsilon_{-3} & \varepsilon_{-3}^2 & \varepsilon_{-3}^3 & \varepsilon_{-3}^{4} \\ 1 & \varepsilon_{-4} & \varepsilon_{-4}^{2} & \varepsilon_{-4}^{3} & \varepsilon_{-4}^{4} \end{array} \right) \ .

Матрица F_{} участвует также и в формуле для определения коэффициентов c_0,c_1,c_2,c_3,c_4 произведения h(x)= f(x)g(x):

\left(\begin{array}{l} f(\varepsilon_0)g(\varepsilon_0) \\ f(\varepsilon_{-1})g(\varepsilon_{-1})\\ f(\varepsilon_{-2})g(\varepsilon_{-2}) \\ f(\varepsilon_{-3})g(\varepsilon_{-3}) \\ f(\varepsilon_{-4})g(\varepsilon_{-4}) \end{array} \right) = F \left(\begin{array}{l} c_0 \\ c_1\\ c_2 \\ c_3 \\ c_4 \end{array} \right) \quad \Rightarrow \quad \left(\begin{array}{l} c_0 \\ c_1\\ c_2 \\ c_3 \\ c_4 \end{array} \right)=F^{-1} \left(\begin{array}{l} f(\varepsilon_0)g(\varepsilon_0) \\ f(\varepsilon_{-1})g(\varepsilon_{-1})\\ f(\varepsilon_{-2})g(\varepsilon_{-2}) \\ f(\varepsilon_{-3})g(\varepsilon_{-3}) \\ f(\varepsilon_{-4})g(\varepsilon_{-4}) \end{array} \right) \ .

Оказывается, что для обращения матрицы F_{} фактически достаточно переставить ее строки:

\left( \begin{array}{lllll} 1 & 1 & 1 & 1 & 1 \\ 1 & \varepsilon_{-1} & \varepsilon_{-1}^2 & \varepsilon_{-1}^3 & \varepsilon_{-1}^{4} \\ 1 & \varepsilon_{-2} & \varepsilon_{-2}^2 & \varepsilon_{-2}^3 & \varepsilon_{-2}^{4} \\ 1 & \varepsilon_{-3} & \varepsilon_{-3}^2 & \varepsilon_{-3}^3 & \varepsilon_{-3}^{4} \\ 1 & \varepsilon_{-4} & \varepsilon_{-4}^{2} & \varepsilon_{-4}^{3} & \varepsilon_{-4}^{4} \end{array} \right)^{-1}=\frac{1}{5} \left( \begin{array}{lllll} 1 & 1 & 1 & 1 & 1 \\ 1 & \varepsilon_{1} & \varepsilon_{1}^2 & \varepsilon_{1}^3 & \varepsilon_{1}^{4} \\ 1 & \varepsilon_{2} & \varepsilon_{2}^2 & \varepsilon_{2}^3 & \varepsilon_{2}^{4} \\ 1 & \varepsilon_{3} & \varepsilon_{3}^2 & \varepsilon_{3}^3 & \varepsilon_{3}^{4} \\ 1 & \varepsilon_{4} & \varepsilon_{4}^{2} & \varepsilon_{4}^{3} & \varepsilon_{4}^{4} \end{array} \right) \ .

В результате получаем ответ h(x)=-2-7\,x-7\,x^2-4\,x^3+10\,x^4, который соответствует реальности.

Обратим внимание на то, что в рассмотренном примере все этапы решения свелись к умножению столбцов на некоторую специфическую матрицу. Формализуем это посредством определения.


Пусть N_{} — произвольное натуральное число. Говорят, что cистема (последовательность) чисел \{ \mathfrak X_0,\dots, \mathfrak X_{N-1} \} получается из системы (последовательности) чисел \{ \mathfrak x_0,\dots, \mathfrak x_{N-1} \} в результате дискретного преобразования Фурье (ДПФ)1), если

\left( \begin{array}{l} \mathfrak X_0 \\ \mathfrak X_1 \\ \vdots \\ \mathfrak X_{N-1} \end{array} \right) = \left( \begin{array}{lllll} 1 & 1 & 1 & \dots & 1 \\ 1 & \varepsilon_{-1} & \varepsilon_{-2} & \dots & \varepsilon_{-N+1} \\ 1 & \varepsilon_{-1}^2 & \varepsilon_{-2}^2 & \dots & \varepsilon_{-N+1}^2 \\ & \dots & & & \dots \\ 1 & \varepsilon_{-1}^{N-1} & \varepsilon_{-2}^{N-1} & \dots & \varepsilon_{-N+1}^{N-1} \\ \end{array} \right)_{N\times N} \left( \begin{array}{l} \mathfrak x_0 \\ \mathfrak x_1 \\ \vdots \\ \mathfrak x_{N-1} \end{array} \right) \quad npu \quad \varepsilon_{-j}= \cos \frac{2\pi j}{N} - \mathbf i \sin \frac{2\pi j}{N} \ .

или, в эквивалентном виде,

\mathfrak X_k=\sum_{j=0}^{N-1} \mathfrak x_j \varepsilon_{-j}^k \quad npu \quad k\in\{0,\dots, N-1\} \ .

Будем обозначать ДПФ буквой \mathcal F:

\mathcal F(\mathfrak x_0,\dots, \mathfrak x_{N-1}) = (\mathfrak X_0,\dots, \mathfrak X_{N-1}) \ ;

оно, очевидно, является линейным преобразованием пространства \mathbb C^N. О последовательности \{\mathfrak x_0,\dots, \mathfrak x_{N-1} \} иногда2) говорят как о входной последовательности; а о последовательности \{ \mathfrak X_0,\dots, \mathfrak X_{N-1} \} — как о выходной последовательности. С точки зрения линейной алгебры, естественней говорить о векторах-столбцах пространства \mathbb C^{N}_{}:

\mathfrak x= \left( \begin{array}{l} \mathfrak x_0 \\ \vdots \\ \mathfrak x_{N-1} \end{array} \right) \quad u \quad \mathfrak X= \left( \begin{array}{l} \mathfrak X_0 \\ \vdots \\ \mathfrak X_{N-1} \end{array} \right) \ .

Еще раз подчеркнем смысл чисел \{\varepsilon_{-j}\}_{j=0}^{N-1}: это — система всех различных корней N_{}-й степени из 1_{} ; только в отличие от стандартного определения, они пронумерованы в обратном порядке (на единичной окружности комплексной плоскости движемся по часовой стрелке).

Дискретное преобразование Фурье \mathcal F полностью определяется матрицей F_{} преобразования Фурье. Последняя, равно как и ей обратная, строится на основе степеней корня N_{}-й степени из 1_{}, этот корень — либо \varepsilon_{-1} либо \varepsilon_{1} — являются первообразными корнями степени N_{} из 1_{}; существенность этого факта заключается в том, что степени \{\varepsilon^j\}_{j=0}^{N-1} такого элемента все различны и пробегают всё множество значений корней N_{}-й степени из 1_{}. Тот факт, что ДПФ последовательности из N_{} элементов строится на основе корня \varepsilon_{} будем записывать в виде

\mathcal F_{\varepsilon,N} \ ;

о таком преобразовании принято также говорить как о N-точечном преобразовании Фурье, построенном на числе \varepsilon_{}.


Теперь обсудим возможности развития предложенного подхода, а также возникающие при его применении сложности.

Метод позволяет решить более общую задачу, чем умножение полиномов. Рассмотрим полиномы f_{}(x) и g_{}(x) степеней меньших N_{} и поставим задачу вычисления остатка от деления f(x) g(x) на полином x^N-1, или, в модулярном формализме

h(x) = f(x) g(x) \pmod{x^N-1} \ .

Такая операция встречается в теории кодирования и известна там как циклическая свёртка полиномов f_{}(x) и g_{}(x). Тогда предложенная схема вычисления:

посредством дискретных преобразований Фурье, составленных на основе корня \varepsilon_{}=\cos \frac{2\pi}{N}+ \mathbf i \sin \frac{2\pi}{N}, позволяет вычислять подобную свертку. В самом деле, если

f(x) g(x) \equiv h(x) + Q(x) (x^N-1) \ ,

где Q(x) \in \mathbb Z[x] означает частное от деления f(x) g(x) на x^N-1, то на любом корне \varepsilon_{} полинома x^N-1 будет выполняться равенство h(\varepsilon)= f(\varepsilon) g(\varepsilon) и N_{} такими значениями полином h_{}(x) степени < N определяется однозначно.

§

В рассмотренном примере нам повезло произвести все вычисления безошибочно — потому что корни 5-й степени из 1_{} выражаются в квадратичных иррациональностях. Но это не выполняется для корней бóльших степеней, скажем, 7_{}-й. Поэтому в таблицах приходится иметь дело с приближенными значениями этих корней, что приводит к ошибкам округления. Так, в приведенном примере, «честный» ответ при сохранении десяти разрядов после точки будет

h(x)=-\scriptstyle{1.999999998}-(\scriptstyle{6.999999995}+\scriptstyle{0.2\cdot 10^{-8}}\displaystyle{\mathbf i})x-(\scriptstyle{6.999999991}+\scriptstyle{0.1\cdot10^{-8}}\displaystyle{\mathbf i})x^2-(\scriptstyle{3.999999992}+ \scriptstyle{0.3\cdot 10^{-8}}\displaystyle{\mathbf i})x^3 +(\scriptstyle{9.999999976}-\scriptstyle{0.33\cdot 10^{-8}}\displaystyle{\mathbf i}) x^4 \ ,

который — мало того, что не является целочисленным, — он не является вещественным! Как поступать в том случае, когда требуется абсолютная точность? На этот вопрос ответ дается в пункте ЧИСЛОВОЕ ПРЕОБРАЗОВАНИЕ ФЕРМА.

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

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

До сих пор не было сделано ограничений на выбор числа N_{} — преобразование Фурье \mathcal F_{\varepsilon,N} вполне определено как только вычислены корни степени N_{} из единицы. Попробуем упростить вычисления этого преобразования для случая выбора числа N_{} четным.

П

Пример. Для случая N_{}=6 сначала перегруппируем суммы, составляющие ДПФ, отделив слагаемые, содержащие \mathfrak x_j с четными индексами от слагаемых, содержащих \mathfrak x_j с нечетными индексами; потом воспользуемся свойствами корней 6-й степени из единицы:

\varepsilon_{-3}=-1, \varepsilon_{-4}=-\varepsilon_{-1}, \varepsilon_{-5}=-\varepsilon_{-2} \ .
\begin{array}{lccccc} \mathfrak X_0=& (\mathfrak x_0+\mathfrak x_2+\mathfrak x_4) & +(\mathfrak x_1+\mathfrak x_3+\mathfrak x_5) & & \\ \mathfrak X_1=& (\mathfrak x_0+\varepsilon_{-1}^2 \mathfrak x_2+\varepsilon_{-1}^4\mathfrak x_4) & +\varepsilon_{-1}(\mathfrak x_1+\varepsilon_{-1}^2\mathfrak x_3+\varepsilon_{-1}^4\mathfrak x_5) &= (\mathfrak x_0+\varepsilon_{-2} \mathfrak x_2+\varepsilon_{-2}^2\mathfrak x_4) & + \varepsilon_{-1}(\mathfrak x_1+\varepsilon_{-2} \mathfrak x_3+\varepsilon_{-2}^2\mathfrak x_5) \\ \mathfrak X_2=& (\mathfrak x_0+\varepsilon_{-2}^2 \mathfrak x_2+\varepsilon_{-2}^4\mathfrak x_4) & +\varepsilon_{-2}(\mathfrak x_1+\varepsilon_{-2}^2\mathfrak x_3+\varepsilon_{-2}^4\mathfrak x_5) &= (\mathfrak x_0+\varepsilon_{-4} \mathfrak x_2+\varepsilon_{-4}^2\mathfrak x_4) &+ \varepsilon_{-2}(\mathfrak x_1+\varepsilon_{-4} \mathfrak x_3+\varepsilon_{-4}^2\mathfrak x_5) \\ \mathfrak X_3=& (\mathfrak x_0+\mathfrak x_2+\mathfrak x_4) & - (\mathfrak x_1+\mathfrak x_3+\mathfrak x_5) & & \\ \mathfrak X_4=& (\mathfrak x_0+\varepsilon_{-4}^2 \mathfrak x_2+\varepsilon_{-4}^4\mathfrak x_4) & +\varepsilon_{-4}(\mathfrak x_1+\varepsilon_{-4}^2\mathfrak x_3+\varepsilon_{-4}^4\mathfrak x_5) &= (\mathfrak x_0+\varepsilon_{-2} \mathfrak x_2+\varepsilon_{-2}^2\mathfrak x_4) &- \varepsilon_{-1}(\mathfrak x_1+\varepsilon_{-2} \mathfrak x_3+\varepsilon_{-2}^2\mathfrak x_5) \\ \mathfrak X_5=& (\mathfrak x_0+\varepsilon_{-5}^2 \mathfrak x_2+\varepsilon_{-5}^4\mathfrak x_4) &+ \varepsilon_{-5}(\mathfrak x_1+\varepsilon_{-5}^2\mathfrak x_3+\varepsilon_{-5}^4\mathfrak x_5) &= (\mathfrak x_0+\varepsilon_{-4} \mathfrak x_2+\varepsilon_{-4}^2\mathfrak x_4) &- \varepsilon_{-2}(\mathfrak x_1+\varepsilon_{-4} \mathfrak x_3+\varepsilon_{-4}^2\mathfrak x_5) \end{array}

В получившихся формулах наблюдается следующая закономерность: вычисления \mathfrak X_0 и \mathfrak X_3 фактически дублируют друг друга (с различием лишь в знаке одного слагаемого); то же самое можно утверждать относительно пар \mathfrak X_1 и \mathfrak X_4, а также \mathfrak X_2 и \mathfrak X_5. Иными словами, грамотная организация процесса вычисления \mathfrak X_0, \mathfrak X_1, \mathfrak X_2 позволит фактически «даром» получить значения \mathfrak X_3, \mathfrak X_4, \mathfrak X_5.

Рассмотрим теперь повнимательней каждую из скобок ( \ ). Все скобки, содержащие \mathfrak x_0,\mathfrak x_2,\mathfrak x_4, имеют одинаковую структуру: они представляют суммы этих величин, домноженных на степени

\varepsilon_{-2} = \cos \frac{4\pi}{6} - \mathbf i \, \sin \frac{4\pi}{6} = \cos \frac{2\pi}{3} - \mathbf i \, \sin \frac{2\pi}{3} \ .

Последняя величина, с одной стороны является корнем 6_{}-й степени из единицы, а, с другой стороны, является корнем 3_{}-й степени из единицы. Если ввести обозначение этих корней:

\tilde \varepsilon_{0}=1,\ \tilde \varepsilon_{-1}=\cos \frac{2\pi}{3} - \mathbf i \, \sin \frac{2\pi}{3},\ \tilde \varepsilon_{-2}=\cos \frac{4\pi}{3} - \mathbf i \, \sin \frac{4\pi}{3} \ ,

то скобки,содержащие \mathfrak x_0,\mathfrak x_2,\mathfrak x_4, переписываются в виде

\begin{array}{rl} \mathfrak Y_0 = & \mathfrak x_0 +\mathfrak x_2\, \tilde \varepsilon_{0}+\mathfrak x_4\, \tilde \varepsilon_{0}^2, \\ \mathfrak Y_1 = & \mathfrak x_0 +\mathfrak x_2\, \tilde \varepsilon_{-1}+\mathfrak x_4\, \tilde \varepsilon_{-1}^2, \\ \mathfrak Y_2 = & \mathfrak x_0 +\mathfrak x_2\, \tilde \varepsilon_{-2}+\mathfrak x_4\, \tilde \varepsilon_{-2}^2, \end{array}

т.е. они представляют ДПФ набора \mathfrak x_0,\mathfrak x_2,\mathfrak x_4. Отличие от исходной задачи заключается в уменьшении размерности: исходная размерность была шестой, получившаяся — третьей. Совершенно аналогично расправляемся со скобками, содержащими \mathfrak x_1,\mathfrak x_3,\mathfrak x_5:

\begin{array}{rl} \mathfrak Z_0 = & \mathfrak x_1 +\mathfrak x_3\, \tilde \varepsilon_{0}+\mathfrak x_5\, \tilde \varepsilon_{0}^2, \\ \mathfrak Z_1 = & \mathfrak x_1 +\mathfrak x_3\, \tilde \varepsilon_{-1}+\mathfrak x_5\, \tilde \varepsilon_{-1}^2, \\ \mathfrak Z_2 = & \mathfrak x_1 +\mathfrak x_3\, \tilde \varepsilon_{-2}+\mathfrak x_5\, \tilde \varepsilon_{-2}^2; \end{array}

они также представляют ДПФ набора \mathfrak x_1,\mathfrak x_3,\mathfrak x_5. Окончательно:

\begin{array}{ccc} \mathfrak X_0=\mathfrak Y_0+\mathfrak Z_0, & \mathfrak X_1= \mathfrak Y_1+\varepsilon_{-1}\mathfrak Z_1, & \mathfrak X_2= \mathfrak Y_2+\varepsilon_{-2}\mathfrak Z_2, \\ \mathfrak X_3=\mathfrak Y_0-\mathfrak Z_0, & \mathfrak X_4= \mathfrak Y_1-\varepsilon_{-1}\mathfrak Z_1, & \mathfrak X_5= \mathfrak Y_2-\varepsilon_{-2}\mathfrak Z_2. \end{array}

Полученные формулы сводят вычисление шеститочечного ДПФ к двум трехточечным:

\mathcal F_{\varepsilon,6}(\mathfrak x_0,\mathfrak x_1,\mathfrak x_2,\mathfrak x_3,\mathfrak x_4,\mathfrak x_5) \Leftarrow \left\{ \begin{array}{c} \mathcal F_{\tilde \varepsilon,3}(\mathfrak x_0,\mathfrak x_2,\mathfrak x_4) \\ \mathcal F_{\tilde \varepsilon,3}(\mathfrak x_1,\mathfrak x_3,\mathfrak x_5) \end{array} \right.

Идея очевидно распространяется на случай произвольного четного числа N_{}:

\mathfrak X_k = \left\{ \begin{array}{lll} \mathfrak Y_k & +\varepsilon_{-k}\mathfrak Z_k & npu \ k\in \{0,1,\dots,N/2-1\}, \\ \mathfrak Y_{k-N/2}& -\varepsilon_{-(k-N/2)}\mathfrak Z_{k-N/2} & npu \ k\in \{N/2,N/2+1,\dots,N-1\}; \end{array} \right. \quad \varepsilon_{-k}=\cos \frac{2\pi k}{N}- \mathbf i \sin \frac{2\pi k}{N} \ .

Здесь в правых частях стоят ДПФ отдельно от четных разрядов и от нечетных разрядов вектора \mathfrak x_{}:

(\mathfrak Y_0,\mathfrak Y_1, \dots , \mathfrak Y_{N/2-1})= \mathcal F_{N/2} (\mathfrak x_0,\mathfrak x_2,\dots, \mathfrak x_{N-2}),\quad (\mathfrak Z_0,\mathfrak Z_1, \dots , \mathfrak Z_{N/2-1}) =\mathcal F_{N/2} (\mathfrak x_1,\mathfrak x_3,\dots, \mathfrak x_{N-1}) \ .

Число \varepsilon_{-k} называется дополнительным множителем3). Если число N/2 само является четным, то процедура повторяется. Если число N_{} представляет собой степень двойки: N=2^{\eta}, то рекурсивный алгоритм сведет вычисление N_{}-точечного ДПФ к вычислению двухточечного ДПФ. Этот алгоритм вычисления ДПФ называется бинарным прореживанием по времени4).

П

Пример. Для случая N_{}=8 схема вычислений будет следующей:

Подсчитаем теперь стоимость бинарного прореживания в элементарных операциях. Если считать, что (N/2)-точечные ДПФ уже вычислены, то для вычисления N_{}-точечных ДПФ требуется 2\, N операций ( N_{} умножений величин \mathfrak Z_{\ell} на дополнительные сомножители и N_{} сложений). Аналогично, для вычисления каждого (N/2)-точечного ДПФ потребуется 2(N/2) операций; и нам требуется вычисление двух таких ДПФ. Таким образом, в случае N=2^{\eta} общее количество операций будет равно

(2\,N)+2\left(2\frac{N}{2}\right)+4\left(2\frac{N}{4}\right)+\dots \ ,

и сумма завершается через \eta_{} шагов когда мы добираемся до двухточечных ДПФ, каждое из которых стоит две операции:

\mathcal F(\mathfrak x_j, \mathfrak x_k)= (\mathfrak x_j+ \mathfrak x_k,\mathfrak x_j- \mathfrak x_k) \ .

В результате общая стоимость вычислений 2\,N \, \eta = 2\, N \log_2 N, что существенно меньше стоимости вычислений всех чисел \{\mathfrak X_k\}_{k=0}^{N-1} «в лоб» по схеме предыдущего ПУНКТА — легко подсчитать, что та схема «стоит» N^2 умножений и N(N-1) сложений. Например, для N=2^{10} имеем:

N^2=1\, 048\, 576, \quad N \log_2 N = 10\, 240 \ .

Числовое преобразование Ферма

§

Для понимания материалов этого пункта крайне желательно просмотреть раздел ИНДЕКС (дискретный логарифм).

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

Если тождество h(x)\equiv f(x)g(x) справедливо для всех вещественных чисел, то оно остается справедливым если мы рассмотрим его по любому целочисленному модулю5):

h(x)\equiv f(x)g(x) \quad \ \Rightarrow \quad h(x)\equiv f(x)g(x) \pmod{M} \quad npu \ \forall M \in \mathbb N \ .

Обратное, вообще говоря, неверно. Но если мы обладаем дополнительной информацией о том, что коэффициенты полиномов f_{}(x) и g_{}(x) являются положительными, а коэффициенты их произведения h(x), вдобавок, строго меньше M_{}, то эти коэффициенты не изменятся, если взять их по модулю M_{}.

Аналогом уравнения в комплексных числах z^n -1 = 0 является сравнение x^n-1 \equiv 0 \pmod{M}. Возьмем в качестве модуля простое число: M=p. Тогда по этому модулю существуют числа A \in \{1,\dots,p-1\}, имеющие порядок \operatorname{ord}(A) равный произвольному делителю числа p-1:

A^{\operatorname{ord}(A)} \equiv 1 \pmod{p}, \ A^q \not\equiv 1 \pmod{p} \quad npu \ q\in\{1,2,\dots, \operatorname{ord}(A)-1 \}\ .
П

Пример. Пусть p=31. Делителем числа p-1=30 является 5_{}. Число \omega=2 имеет порядок 5_{} и, следовательно, по нему можно построить 5_{}-титочечное преобразование Фурье. Для полиномов f(x)=2+3\,x+5\,x^2 и g(x)=-1-2\,x+2\,x^2 примера из ☞ ПУНКТА получаем:

\left(\begin{array}{l} f(1) \\ f(2)\\ f(4) \\ f(8) \\ f(16) \end{array} \right) = F \left(\begin{array}{l} a_0 \\ a_1\\ a_2 \\ 0 \\ 0 \end{array} \right) \pmod{31},\quad \left(\begin{array}{l} g(1) \\ g(2)\\ g(4) \\ g(8) \\ g(16) \end{array} \right) = F \left(\begin{array}{l} b_0 \\ b_1\\ b_2 \\ 0 \\ 0 \end{array} \right) \pmod{31} \quad npu \ F=\left( \begin{array}{ccccc} 1 & 1 & 1 & 1 & 1 \\ 1 & 2 & 4 & 8 & 16 \\ 1 & 4 & 16 & 2 & 8 \\ 1 & 8 & 2 & 16 & 4 \\ 1 & 16 & 8 & 4 & 2 \end{array} \right) \ .

Значения интерполяционных таблиц:

\begin{array}{c|c|c|c|c|c} x & 1 & 2 & 4 & 8 & 16 \\ \hline f(x) \pmod{31} & 10 & 28 & 1 & 5 & 28 \\ \hline g(x) \pmod{31} & 30 & 3 & 23 & 18 & 14 \\ \hline f(x)g(x) \pmod{31} & 21 & 22 & 23 & 28 & 20 \end{array}

Теперь определяем коэффициенты полинома h_{}(x):

\left(\begin{array}{l} c_0 \\ c_1\\ c_2 \\ c_3 \\ c_4 \end{array} \right)=F^{-1} \left(\begin{array}{l} 21 \\ 22 \\ 23 \\ 28 \\ 20 \end{array} \right) \ .

Имеем \det F=5 \pmod{31} и (\det F)^{-1} =25 \pmod{31}. Так же, как и в случае пункта ☞ ДИСКРЕТНОЕ ПРЕОБРАЗОВАНИЕ ФУРЬЕ, обращение матрицы F_{} (по модулю 31) фактически сводится к перестановке ее строк:

F^{-1} = 25 \left( \begin{array}{ccccc} 1 & 1 & 1 & 1 & 1 \\ 1 & 16 & 8 & 4 & 2 \\ 1 & 8 & 2 & 16 & 4 \\ 1 & 4 & 16 & 2 & 8 \\ 1 & 2 & 4 & 8 & 16 \end{array} \right) \pmod{31} \quad \Rightarrow \quad \left(\begin{array}{l} c_0 \\ c_1\\ c_2 \\ c_3 \\ c_4 \end{array} \right) = \left(\begin{array}{l} 29 \\ 24\\ 24 \\ 27 \\ 10 \end{array} \right) \ .

Но этот результат не совпадает с ответом! В самом деле, при вычислениях по модулю положительного числа M_{} мы договорились результаты представлять в виде чисел множества \{0,1,\dots, M-1 \}, т.е. чисел положительных. Теперь часть этих положительных чисел надо сделать отрицательными, т.е. вычесть из них M_{}. Каким правилом следует при этом руководствоваться? — Можно взять модуль настолько большим, чтобы гарантировать попадание всех коэффициентов \{c_j\} в интервал \big] -(M-1)/2, (M-1)/2 \big]. При таком выборе полученное в ходе модулярных вычислений число \tilde c_j интерпретируется следующим образом:

c_j = \left\{ \begin{array}{ll} \tilde c_j & npu \ 0 \le \tilde c_j \le (M-1)/2, \\ \tilde c_j-M & npu \ (M+1)/2 \le \tilde c_j \le M-1. \end{array} \right.

В рассматриваемом примере следует действовать именно по этому правилу (если мы уже установили, что коэффициенты полинома находятся в интервале ]-15,15]).

§

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

Итак, весь алгоритм точного (безошибочного) вычисления произведения полиномов (циклической свертки) свелся к преобразованию Фурье, действующему также как и в пункте ☞ ДИСКРЕТНОЕ ПРЕОБРАЗОВАНИЕ ФУРЬЕ, но только теперь по модулю простого числа p_{}:

\mathcal F_{\omega, N} \pmod{p} \ .

Число N_{} в этом случае должно быть делителем числа p-1, а число \omega_{} \in \{2,\dots p-1\} — иметь порядок равный N_{}.

Теперь посмотрим как можно сэкономить на вычислениях.

П

Пример. Пусть p=257. Делителями числа p-1=256 являются степени 2_{}. Так любое число \omega\in \{4,\ 64,\ 193,\ 253 \} имеет порядок 8_{}:

\omega^8 \equiv 1 \pmod{257} \ .

Отсюда, в силу теоремы 1, приведенной ☞ ЗДЕСЬ,

\omega^4\equiv_{_{257}} -1,

ну и, дальше понятно, что

\omega^5\equiv_{_{257}} -\omega, \ \omega^6\equiv_{_{257}} -\omega^2, \ \omega^7\equiv_{_{257}} -\omega^3,\dots .

8_{}-точечное ДПФ \mathcal F_{\omega, 8} строится на числе \omega_{} по формулам:

\begin{array}{llllllll} \mathfrak X_0=\mathfrak x_0 &+ \mathfrak x_1&+ \mathfrak x_2 & + \mathfrak x_3 & + \mathfrak x_4 &+ \mathfrak x_5&+ \mathfrak x_6& + \mathfrak x_7, \\ \mathfrak X_1=\mathfrak x_0 &+ \mathfrak x_1\omega& + \mathfrak x_2\omega^2 & + \mathfrak x_3\omega^3 & - \mathfrak x_4 &- \mathfrak x_5\omega&- \mathfrak x_6\omega^2& - \mathfrak x_7\omega^3, \\ \mathfrak X_2=\mathfrak x_0 &+ \mathfrak x_1\omega^2& - \mathfrak x_2 & - \mathfrak x_3\omega^2 &+ \mathfrak x_4 &+ \mathfrak x_5\omega^2&-\mathfrak x_6& -\mathfrak x_7\omega^2, \\ \mathfrak X_3=\mathfrak x_0 &+ \mathfrak x_1\omega^3& - \mathfrak x_2\omega^2 & + \mathfrak x_3\omega &- \mathfrak x_4 &- \mathfrak x_5\omega^3&+\mathfrak x_6\omega^2& -\mathfrak x_7\omega, \\ \mathfrak X_4=\mathfrak x_0 &- \mathfrak x_1&+ \mathfrak x_2 & - \mathfrak x_3 & + \mathfrak x_4 &- \mathfrak x_5&+ \mathfrak x_6& - \mathfrak x_7, \\ \mathfrak X_5=\mathfrak x_0 &- \mathfrak x_1\omega& + \mathfrak x_2\omega^2 & - \mathfrak x_3\omega^3 &- \mathfrak x_4 &+ \mathfrak x_5\omega&-\mathfrak x_6\omega^2& -\mathfrak x_7\omega^3, \\ \mathfrak X_6=\mathfrak x_0 &- \mathfrak x_1\omega^2& - \mathfrak x_2 & + \mathfrak x_3\omega^2 &+ \mathfrak x_4 &- \mathfrak x_5\omega^2&-\mathfrak x_6& +\mathfrak x_7\omega^2, \\ \mathfrak X_7=\mathfrak x_0 &- \mathfrak x_1\omega^3& - \mathfrak x_2\omega^2 & - \mathfrak x_3\omega &- \mathfrak x_4 &+ \mathfrak x_5\omega^3&+\mathfrak x_6\omega^2& +\mathfrak x_7\omega. \end{array}

Обратим внимание на то, что вычисления чисел \mathfrak X_0 и \mathfrak X_4 фактически дублируют друг друга:

\begin{array}{lll} \mathfrak X_0 &= (\mathfrak x_0 + \mathfrak x_2+ \mathfrak x_4 + \mathfrak x_6) & + (\mathfrak x_1 + \mathfrak x_3+ \mathfrak x_5 + \mathfrak x_7), \\ \mathfrak X_4 &= (\mathfrak x_0 + \mathfrak x_2+ \mathfrak x_4 + \mathfrak x_6) & - (\mathfrak x_1 + \mathfrak x_3+ \mathfrak x_5 + \mathfrak x_7); \end{array}

аналогично \mathfrak X_1 и \mathfrak X_5:

\begin{array}{lll} \mathfrak X_1&=(\mathfrak x_0 + \mathfrak x_2\omega^2 - \mathfrak x_4- \mathfrak x_6\omega^2) &+ \omega(\mathfrak x_1+ \mathfrak x_3\omega^2 -\mathfrak x_5- \mathfrak x_7\omega^2),\\ \mathfrak X_5&=(\mathfrak x_0 + \mathfrak x_2\omega^2 - \mathfrak x_4- \mathfrak x_6\omega^2) &- \omega(\mathfrak x_1+ \mathfrak x_3\omega^2 -\mathfrak x_5- \mathfrak x_7\omega^2); \end{array}

\mathfrak X_2 и \mathfrak X_6:

\begin{array}{lll} \mathfrak X_2&=(\mathfrak x_0 - \mathfrak x_2+ \mathfrak x_4-\mathfrak x_6) & + \omega^2 (\mathfrak x_1-\mathfrak x_3+ \mathfrak x_5-\mathfrak x_7),\\ \mathfrak X_6&=(\mathfrak x_0 - \mathfrak x_2+ \mathfrak x_4-\mathfrak x_6) & - \omega^2 (\mathfrak x_1-\mathfrak x_3+ \mathfrak x_5-\mathfrak x_7); \end{array}

и \mathfrak X_3 и \mathfrak X_7:

\begin{array}{lll} \mathfrak X_3&=(\mathfrak x_0 - \mathfrak x_2\omega^2 - \mathfrak x_4+\mathfrak x_6\omega^2) &+\omega^3(\mathfrak x_1- \mathfrak x_3\omega^2 -\mathfrak x_5 + \mathfrak x_7\omega^2 ), \\ \mathfrak X_7&=(\mathfrak x_0 - \mathfrak x_2\omega^2 - \mathfrak x_4+\mathfrak x_6\omega^2) &-\omega^3(\mathfrak x_1- \mathfrak x_3\omega^2 -\mathfrak x_5 + \mathfrak x_7\omega^2 ). \end{array}

Введем обозначения для скобок ( \ \ ):

\begin{array}{rl} \mathfrak Y_0 = & \mathfrak x_0 + \mathfrak x_2+ \mathfrak x_4 + \mathfrak x_6, \\ \mathfrak Y_1 = & \mathfrak x_0 + \mathfrak x_2\omega^2 + \mathfrak x_4 \omega^4 + \mathfrak x_6\omega^6, \\ \mathfrak Y_2 = & \mathfrak x_0 + \mathfrak x_2\omega^4 + \mathfrak x_4 \omega^8 + \mathfrak x_6\omega^{12}, \\ \mathfrak Y_3 = & \mathfrak x_0 + \mathfrak x_2\omega^6 + \mathfrak x_4 \omega^{12} + \mathfrak x_6\omega^{18}, \end{array} \qquad \begin{array}{rl} \mathfrak Z_0 = & \mathfrak x_1 + \mathfrak x_3+ \mathfrak x_5 + \mathfrak x_7, \\ \mathfrak Z_1 = & \mathfrak x_1 + \mathfrak x_3\omega^2 + \mathfrak x_5 \omega^4 + \mathfrak x_7\omega^6, \\ \mathfrak Z_2 = & \mathfrak x_1 + \mathfrak x_3\omega^4 + \mathfrak x_5 \omega^8 + \mathfrak x_7\omega^{12}, \\ \mathfrak Z_3 = & \mathfrak x_1 + \mathfrak x_3\omega^6 + \mathfrak x_5 \omega^{12} + \mathfrak x_7\omega^{18}. \end{array}

В результате 8_{}-миточечное ДПФ можно представить формулами

\begin{array}{cccc} \mathfrak X_0=\mathfrak Y_0+\mathfrak Z_0, & \mathfrak X_1= \mathfrak Y_1+\omega\mathfrak Z_1, & \mathfrak X_2= \mathfrak Y_2+\omega^2\mathfrak Z_2, & \mathfrak X_3= \mathfrak Y_3+\omega^3\mathfrak Z_3, \\ \mathfrak X_4=\mathfrak Y_0-\mathfrak Z_0, & \mathfrak X_5= \mathfrak Y_1-\omega\mathfrak Z_1, & \mathfrak X_6= \mathfrak Y_2-\omega^2\mathfrak Z_2, & \mathfrak X_7= \mathfrak Y_3-\omega^3\mathfrak Z_3. \end{array}

Но каждая из последовательностей (\mathfrak Y_0,\mathfrak Y_1,\mathfrak Y_2,\mathfrak Y_3) и (\mathfrak Z_0,\mathfrak Z_1,\mathfrak Z_2,\mathfrak Z_3) также представляет собой ДПФ — только, в отличие от исходного, это преобразование — 4_{}-хточечное, составленное относительно числа

\omega^2 \in \{16 \equiv_{_{257}} 4^2 \equiv_{_{257}} 253^2,\ 241\equiv_{_{257}} 64^2 \equiv_{_{257}} 193^2 \} \ ,

имеющего порядок равный 4_{}:

(\mathfrak Y_0,\mathfrak Y_1,\mathfrak Y_2,\mathfrak Y_3)= \mathcal F_{\omega^2, 4}(\mathfrak x_0, \mathfrak x_2, \mathfrak x_4, \mathfrak x_6) \pmod{257},\quad (\mathfrak Z_0,\mathfrak Z_1,\mathfrak Z_2,\mathfrak Z_3)= \mathcal F_{\omega^2, 4}(\mathfrak x_1, \mathfrak x_3, \mathfrak x_5, \mathfrak x_7) \pmod{257} \ .










Статья не закончена!

Литература

Блейхут Р. Быстрые алгоритмы цифровой обработки сигналов. М. Мир. 1989

1) Discrete Fourier Transform (DFT).
2) Такая терминология принята в ТЕОРИИ СИГНАЛОВ.
3) Twiddle factor.
4) Radix-2 decimation-in-time (DIT), что дословно переводися как «бинарная децимация по времени»; не очень удачное, на мой взгляд, название, поскольку латинское слово decimatio означает казнь каждого десятого, а вовсе не второго
5) В первой формуле \equiv_{} означает тождественное равенство, а во второй — сравнение по модулю.

2016/08/21 09:05 редактировал au