будет содержаться текущий номер уровня в дереве рекурсии (растущий от s = 1 для листьев до lgra для корня, когда мы получаем га-элементный результат). Вот схема программы:

1for $s \leftarrow 1$ to $\lg n$

2\quad do for $k \leftarrow 0$ to $n-l$ by $2~s$

3\qquad do преобразовать два $2~{з-1}$-элементных куска \quad\qquad $A[k..k+2~{s-l}-l]$ и $A[k+2~{s-l>..k+2~s-l]$ \quad\qquad в $2~з$-элементный кусок $A[k..k+2~s-l]$

Опишем тело цикла (строка 3) более подробно. Мы воспроизводим цикл for из процедуры Recursive-FFT, используя вместо у И кусок A[k..k + 2S~1 - 1], а вместо у№ - кусок А[к + 25"1 ..к + 2S - 1]. Значение и (используемое в преобразовании бабочки) зависит от s: оно есть ujm, где т = 2s. Временная переменную и используется в преобразовании бабочки. Таким образом, развёртывая строку 3 в предыдущей процедуре, получаем такую программу:

\textsc{FFT-Base}

1$п \leftarrow length[а]$ \qquad $п$ есть степень 2

2for $s \leftarrow 1$ to $\lg n$

3\quad do $m \leftarrow 2~s$

4\qquad $\omega m \leftarrow e~{2\pi i/m}$

5\qquad for $k \leftarrow 0$ to $n-l$ by m

6\qquad\quad do $\omega \leftarrow 1$

7\qquad\qquad for $j \leftarrow 0$ to $m/2-l$

8\qquad\qquad\quad do $t \leftarrow \omega A[k+j+m/2]$

9\qquad\qquad\qquad $u \leftarrow A[k+j]$

10\qquad\qquad\qquad $A[k+j] \leftarrow u + t$

11\qquad\qquad\qquad $A[k+j+m/2] \leftarrow u - t$

12\qquad\qquad\qquad $\omega \leftarrow \omega \omega m$

Наконец, представим окончательную версию итеративной процедуры быстрого преобразования Фурье, в которой два внутренних цикла переставлены (экономия при вычислении uj) и использована дополнительная процедура Bit-Reverse-Copy(<2, А), копирующая элементы вектора а в массив А в нужном нам порядке.

\textsc{Iterative-FFT>$(a)$

1\textsc{Bit-Reverse-Copy>$(a, А)$

2$п \leftarrow length[а]$ \quad $п$ --- степень $2$

3for $s \leftarrow 1$ to $\lg n$

4\quad do $m \leftarrow 2~s$

5\qquad $\omega m \leftarrow e~{2\pi i/m}$

6\qquad $\omega \leftarrow 1$

7\qquad for $j \leftarrow 0$ to $m/2-l$

8\qquad\quad do for $k \leftarrow j$ to $n-l$ by $m$

9\qquad\qquad\quad do $t \leftarrow \omega A [k+m/2]$

10\qquad\qquad\qquad $u \leftarrow A[k]$

11\qquad\qquad\qquad $A[k] \leftarrow u+t$

12\qquad\qquad\qquad $A[k+m/2] \leftarrow u-t$

13\qquad\qquad $\omega \leftarrow \omega \omega m$

Каким образом процедура Bit-Reverse-Copy переставляет элементы входного вектора а, помещая их в массив А? Посмотрев на листья дерева на рис. 32.4, можно заметить, что они расположены «перевёрнутом двоичном» порядке ("bit-reverse binary"). Объясним, что имеется в виду. Через rev (к) при к = 0,1, 2,... , 2п - 1 обозначим (lg га)-битовое целое число, двоичная запись которого получается, если написать двоичную запись числа к (всего lgra битов, включая начальные нули, если они есть) «задом наперёд». Так вот, процедура Bit-Reverse-Copy помещает элемент ак в А[геи(&)]. На рисунке 32.4, например, листья ао,а\,а2, ,а? появляются в позициях 0, 4, 2, 6, 1, 5, 3, 7; в двоичной системе номера позиций записываются как ООО, 100, 010, 110, 001, 101, 011, 111. (Отметим, что rev (rev (к)) = к.)

Чтобы понять, почему листья идут в перевёрнутом двоичном порядке, заметим, что в корне дерева налево идут индексы с нулевым младшим битом, а направо - с единичным, что после обращения битов соответствует обычному порядку (сначала числа с нулевым старшим битом, потом с единичным). Отбрасывая на каждом уровне младший бит, мы продолжаем этот процесс вниз по дереву, и в конце концов получаем перевёрнутый двоичный порядок в листьях.

Поскольку функцию rev(k) легко вычислять, можно записать процедуру Bit-Reverse-Copy (копирование в перевёрнутом двоичном порядке) следующим образом:

\textsc{Bit-Reverse-Copy}$(a,А)$

1$n \leftarrow length[а]$

2for $k \leftarrow 0$ to $n-l$

3\quad do $A [{\mathrm rev}(k)] \leftarrow a k$

Построенная итеративная реализация быстрого преобразования Фурье требует времени О (га lgra). В самом деле, вызов Bit-Reverse-Copy заведомо укладывается в О (га lgra) операций, поскольку цикл исполняется га раз, а целое число от 0 до га - 1 занимает lgra битов, которые могут быть обращены за время О (lgra). (На практике мы обычно заранее знаем значение га, поэтому можем изготовить таблицу значений функции rev (к) заранее и пользоваться этой таблицей. Можно также использовать приём из задачи 18-1, позволяющий обратить все га чисел от 0 до га - 1 за О (га) операций.)

Мы должны ещё проверить, что остальная часть процедуры Iterative-FFT требует времени ©(гаlgra), Это ясно: на каждом из lgra уровней дерева рекурсии (рис. 32.4) имеется га чисел, на

Рис. 31.6 32.5 Схема parralel-FFT, вычисляющая преобразование Фурье для п = 8 входов. Три уровня соответствуют значениям s = 1,2,3. Для п входов подобная схема имеет глубину 0(lgn) и размер 0(nlgn).

получение каждого из них расходуется 0(1) операций.

Параллельная схема быстрого преобразования Фурье

Изложенные идеи позволяют построить схему из функциональных элементов, эффективно выполняющую преобразование Фурье. (Схемы из функциональных элементов описаны в главе 29; в нашем случае элементы будут выполнять арифметические операции.)

Схема, выполняющая преобразование Фурье для га входов, которую мы назовём Parallel-FFT, изображена на рис. 32.5 (для п = 8). Значения на вход поступают в перевёрнутом двоичном порядке; схема содержит lg га каскадов, в каждом из которых параллельно выполняются га/2 преобразований бабочки Таким образом, глубина схемы есть в (lgra).

Такое распараллеливание возможно, поскольку на каждом уровне рекурсии (рис. 32.4) имеется га/2 независимых преобразований бабочки, которые можно выполнять параллельно.

Упражнения

32.3-1

Покажите работу процедуры Iterative-FFT на примере входного вектора (0, 2, 3, -1,4, 5, 7, 9). 32.3-2

Как реализовать алгоритм быстрого преобразования Фурье, выполняя перестановку в перевёрнутом двоичном порядке в конце, а не в начале? (Указание. Рассмотрите обратное преобразование Фурье.)

32.3-3

Сколько проводов, элементов сложения, вычитания и умножения имеется в схеме Parallel-FFT описанной в этом разделе? (Каждый провод соединяет один выход с одним или несколькими входами.)

32.3-4*

Предположим, что сумматоры в схеме Parallel-FFT иногда ломаются и выдают на выходе ноль независимо от входных значений. Предположим, что сломался ровно один сумматор, и мы не знаем, какой. Как можно быстро его найти, подавая на вход схемы разные наборы значений и наблюдая за выходными значениями?

Задачи

32-1 Умножение методом разделяй и властвуй

a.Покажите, как умножить два линейных многочлена ах + Ь и cx-\-d, сделав только три умножения. (Указание: одним из умножений будет (а + Ъ) (с + d).)

b.Предложите два алгоритма для вычисления произведения двух