и V {х, у)1к t (х, у) представляют собой многочлены от двух переменных, задаваемые входными и выходными данными соответственно равенствами

(Х !/)= Ti S хУ

Нх. !/)= S S t rxy-.

Чтобы задать многочлен v (х, у) по компонентам vt+i, нужно выполнить только перестановки. Аналогично, для вычисления компонент 2+1 по многочлену t (х, у) нужны только перестановки. Структуру легче увидеть, если переписать многочлены в виде к (X, у) = V, (X) + № (х), i (х, у) = t, (X) + yi, (X) -

g{x.y)= S

где g{x) = £ Тогда

xy = g{x) + yg-(x).

Если входные данные вещественны, то t, (х) = tq (х) и вычислять надо только

/о W = g W V, (х) + g* (х) И1 (х) (mod X * 1). На следующем щаге воспользуемся разложением х / - 1 = (х /- 1)(х + 1) и китайской теоремой об остатках. Пусть

g )(x) = g(x) (modx- - I), g> (х) = g (х) (mod x + 1), и остальные необходимые многочлены определены аналогичным образом. Так как согласно теореме 4.6.3 g (х) = О, то вычислять надо только члены, содержащие g** (х):

Далее воспользуемся алгоритмом 2-точечной циклической свертки

\i[b (x) +

f (JT) - 1

рис. 10.5- Рекурсивная форма БПФ-алгоритма Винограда по основанию два.

.Многочлен 4 *) + ()] содержит только вещественные коэффициенты, а многочлен -[ W -g**(x)] содержит только чисто мнимые коэффициенты. Следовательно, задача свелась к вычислению четырех произведений по модулю х + 1 многочленов с вещественными коэффициентами. Теоретически, каждая из этих задач требует п/4 - 1 вещественных умножений; практически из-нестные алгоритмы требуют несколько большего числа умножений.

На рис. 10.5 приведена блок-схема описанного рекурсивного вычисления; на рис. 10.6 затабулированы получаемые прн этом характеристики.

Рис, 10.6. Характеристики некоторых БПФ-алгоритмов Винограда по основанию два.

Вход в 1 -точечиое БПФ

Устаиобить вычисления

Вызоб 2 -тоиечиого БПФ

Вызов 2 --тоиециого БПФ

Мкратный вызоВ полиномиального произведений по модулю x *f

Сочетание частичных результатов

Выход

Вещественные входные да1 Число вещественных

356 Гл. Ю. Быстрые алгоритмы, основанные на стратегии дублирования

10.5. Быстрая транспозиция

Во время обработки больших двумерных таблиц, подобных оцифрованным изображениям, процессор может воспринять только малую часть этой таблицы. Например, для запоминания (1024X X 1024)-таблицы требуется более, чем миллион ячеек памяти и более, чем два миллиона, если таблица комплексная. В большинстве случаев таблица запоминается во внешней памяти системы и считывается по частям в оперативную память процессора. Как правило, (л X п)-таблица запоминается по столбцам (или по строкам) во внешней памяти, и в каждый момент времени между оперативной памятью и внешней памятью происходит обмен одним столбцом. Обмен двух элементов из двух различных столбцов столь же трудоемок, сколь обмен двух целых столбцов. Таким образом, для формирования одной строки матрицы приходится считывать п столбцов, а для формирования всех строк приходится считывать столбцов, так что прямое транспонирование матрицы приводит к считыванию п столбцов,

В тех приложениях, в которых оперативная память мала, для перестановки строк и столбцов во внешней памяти можно пользоваться алгоритмом быстрой транспозиции. Мы рассмотрим случай, когда оперативная память допускает запоминание двух столбцов из квадратной таблицы. Это наиболее интересный и позволяющий продемонстрировать все основные идеи случай. Для транспонирования (я X я)-таблицы во внешней памяти алгоритм делает п logj я считываний столбцов в оперативную память. Если в оперативную память нельзя записать двух столбцов, то алгоритм приходится дублировать, что снижает его эффективность; если оперативная память допускает запись более двух столбцов, то алгоритм можно улучшить на некоторую константу.

Транспонирование матрицы, записанной в оперативной памяти прямого доступа, является тривиальной задачей, так как реализуется простым изменением адресов при вызове данных. Мы будем полагать, что транспонирование малых матриц сводится к считыванию таблицы в оперативную память и последующему считыванию ее во внешнюю память. Пусть (2 Х2 )-матрица А разбита на блоки размерности 2 - в виде

Ala А21 Аи,

Предположим, что у нас имеется способ вычисления матрицы Аи Ап]

А21 Л2:

10.6. Умножение матриц

Тогда для вычисления матрицы ГА(

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

Теперь применим эту идею рекурсивно к вычислению Aui A[i. А?2 и аГз, разбивая каждую т них на подблоки. Так как An !i Ал расположены в одних и тех же столбцах таблицы, то входящие в них транспозиции выполняются перестановкой одних и тех же столбцов. Такнм образом, каждый уровень рекурсии затрагивает п столбцов, а всего уровней имеется loga п. Следовательно, в противоположность прямому методу транспонирования матрицы, содержащему перестановок столбцов, построенный быстрый алгоритм транспозиции содержит п logg п перестановок столбцов.

10.6. Умножение матриц

Произведение матриц размерности два

можно вычислить по правилам

и = 11*11 + <l2*21.

См = fiaibii + 22*21,

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

семь умножений,

В алгоритме Штрассена сначала выполняются следующие

вычисления:

mi = (аг - Щг) (b-ti + bn). mj = (а -f aja) (йц + Ьгг), ins = (a,i - Oji) (Ьц + bia), 4 = (ill + Ois) *2!.

т, = (bj, - 6 ).

m, = (da + ajj) йц-Затем элементы произведения матриц вычисляются по формулам

Сц = m-i + Шг - nil + ni

Сп = mi + mi,

Си = me + m

= ma - тз + - m. В матричном виде алгоритм Штрассена записывается равенством !j

1

[о 1-1

о -1 о 1

0 1 о]с

1 о о:

о 1 о о

1 о о 16,

0 о 1

1 0-1

о 1 о

в котором элементы стоящей в центре диагональной матрицы вычисляются по правилу

Алгоритм Штрассена содержит семь умножений и 18 сложений. Если одна из двух перемножаемых матриц является постоянной и используется много раз, то некоторые сложения можно вычислить заранее вне процесса умножения и тогда число сложений становится равным 13.

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

Теперь рассмотрим задачу умножения двух матриц размерности п. Прямой метод-такого вычисления требует п умножений и (п - 1) п сложений. Предположим, что п равно степени двух: п = 2 для некоторого целого т; в противном случае дополним матрицу справа таким количеством нулевых столбцов и снизу таким количеством нулевых строк, чтобы п стало равным степени двух.

Произведение матриц С = АВ можно разбить на блоки

где каждый блок представляет собой матрицу размерности 2 -. Вычисление блоков слева прямым умножением блоков справа содержит восемь умножений матриц размерности п/2 и четыре сложения матриц размерности л/2. Если все эти вычисления проводить стандартными способами, то полное число умножений будет равно

M{n) = 8(Jy = п\

А{п)8(-1)(У + 4(У=(п-1)п

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

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

Ml = (Ai2 - Am) (Вг1 + Вм),

Мз = (А1, - AjiKBu-f Bi,), Ml = (All + Ак) Вгг, Mb = Aii(Bi2-B,J, Mj = Az2 (Bji - Bii), M, = (A2i + A22)Bu.

Б.чоки матрицы С вычисляются по правилам

Си = Ml + - + Ме,

Ci3 = Ml + Мб,

Cjl = Mg -t- M7,

Саз = Ma - M3 M5 - M7.

Если матрица A является матрицей констант, то вместо имевшихся ранее восьми умножений и четырех сложений матриц размерности п/2, мы получаем семь умножений и тринадцать сложений матриц размерности п/2. Если каждое из этих вычислений проводить стандартным способом, то для полного числа умножений получаем