КОМПЬЮТЕРНАЯ ОПТИКА И ОБРАБОТКА ИЗОБРАЖЕНИЙ
УДК 512.643:519.712.2:621.391
СИНТЕЗ ПАРАЛЛЕЛЬНЫХ АЛГОРИТМОВ ПРЕОБРАЗОВАНИЙ ФУРЬЕ-ГАЛУА В ПРЯМЫ1Х СУММАХ КОНЕЧНЫ1Х КОЛЕЦ
© 2000 В.М. Чернов
Институт систем обработки изображений РАН, г. Самара
Рассматривается метод синтеза быстрых алгоритмов дискретных преобразований, позволяющих, в частности, безошибочно и с минимальным числом умножений вычислять дискретные свертки целочисленных последовательностей. Предложенный подход базируется на неоднозначности разложения на простые множители в алгебраических кольцах. Выбор подходящего разложения определяется особенностями машинной реализации модулярных вычислений.
Введение
Модулярные аналоги дискретного преобразования Фурье (теоретико-числовые преобразования, ТЧП, преобразования Фурье-Галуа)
N-1
X(m) = ^ x(n)wmn (modp), wN ° l(modp); (1)
n =0
были впервые введены в работах Р.Г. Фарад-жиева и Я.З.Ципкина [1, 2], А.Штрассена и В. Шёнхаге [3] и переоткрыты в работе Ч.М. Рейдера [4]. Последняя работа получила наибольшую известность и явилась основой для различных модификаций и обобщений ТЧП. Наибольшая популярность этой тематики приходится на 70-е годы, что объясняется, в основном, вычислительными преимуществами целочисленной арифметики для использовавшихся цифровых устройств и возможностями эффективного решения массовых задач информатики - “безошибочного” вычисления циклической свертки (целочисленных) последовательностей [5] и быстрого умножения больших целых чисел [3]. Возрождение интереса к ТЧП, наметившееся в последние годы, связано с разработкой нового поколения СБИС-устройств, использующих модулярные вычисления для реализации арифметических операций [6 - 9].
Наиболее существенным недостатком ТЧП является то, что простые числа p с “удобной” машинной реализацией (простые числа Мерсенна, Ферма, Голомба и т.п) встречаются в натуральном ряду достаточно редко [10], что существенно ограничивает возможности ТЧП, например, в задачах цифровой обработки многомерных цифровых массивов
(в частности, изображений).
В настоящей работе рассматривается метод синтеза параллельных алгоритмов “безошибочного” вычисления сверток целочисленных последовательностей с использованием “ТЧП-подобных” преобразований, реализуемых в кольцах классов вычетов по модулю составных чисел. Метод базируется на двух независимых и хорошо апробированных идеях:
- вложение кольца классов вычетов по (составному) (mod m) в прямую сумму некоторых конечных колец;
- представление данных в «нетрадиционных» системах счисления.
Существенно новым в рассмотренном методе является согласованный выбор конечных колец, способа вложения кольца классов вычетов в прямую сумму этих колец, а также систем счисления «с иррациональным основанием», наследующих свойства двоичной системы счисления в кольцах классов вычетов по модулям простых чисел Ферма. Аналогичная задача для вычислений в кольцах классов вычетов по модулям составных чисел Мерсенна рассмотрена автором в работе [13].
Теоретико-числовое преобразование Ферма
Исторически первым ТЧП, нашедшим эффективные применения в задачах цифровой обработки сигналов, явилось преобразование (1) при
р = fs = 22 +1 = 2B +1, s = 0,1, 2, 3,4. (2)
Числа f называются числами Ферма и
являются простыми только при пяти указанных значениях s. Преобразование (1) при p = fs называется ТЧП Ферма или просто преобразованием Ферма.
Привлекательность преобразования Ферма обуславливается по меньшей мере двумя факторами:
- существование удобного для машинной реализации “битового” представления элементов поля (mod f) с просто реализуемыми операциями сложения и умножения элементов.
- наличие хорошей алгоритмической поддержки вычисления преобразования (1) в виде структурно простых модулярных аналогов классического быстрого алгоритма Кули-Тьюки дискретного преобразования Фурье.
Действительно, любое целое число х из диапазона 0 < х < fs - 1 =2B может быть представлено (В+1)-битовым разложением
x = xB 2B + xB-12B-1 +... + x121 + x020,
x=0,i. (3)
Ассоциированный с разложением (3) (В+1)-битовый вектор
< x >= (xB,xB-1,... xi,Xo) . (4)
будем называть кодом элемента х (modf).
Операции в кольце классов вычетов (mod f,) индуцируют формальные правила преобразования кодов, вполне определяемые соотношением
2 B =-1 (mod f ). (5)
Именно, сложение производится по почти обычным правилам двоичного сложения “с переносом в старший разряд”; при переполнении В-го разряда “лишняя” единица вычитается из В-разрядной части результата (или, в терминах кодов, “с инвертированным знаком переносится в самый младший разряд кода); умножение на 2 элемента х индуцирует преобразование кодов по правилу
< x >= (хв» xb-i»- xi» xc) a
a < 2x >= (x5-i,... xi,x^-x*); (5)
умножение элементов общего вида сводится к сложениям и умножением на степени 2.
Заметим, что код в правой части соот-
ношения (5) не обязательно является битовым вектором, но может быть легко преобразован к битовому виду с использованием тривиального соотношения
xj+12j+1 = xj+1^+ xj+12j (6)
Такие векторы, преобразованные по правилу (6) к битовому представлению, будем называть редуцированными кодами и обозна___ -*• *
чать < x > .
Далее, в отличие от поля комплексных чисел, в конечном поле (mod p) существуют корни не любой степени N единицы, а только удовлетворяющие условию делимости:
N || (p-1). (7)
Для чисел Ферма это стеснительное, в общем случае, ограничение (7) гарантирует существование структурно простых быстрых алгоритмов вычисления преобразования (1).
Наиболее просто реализуется преобразование (1) при w ° 2 (modf). В этом случае умножения на фазовые множители в модулярной версии алгоритма Кули-Тьюки (БПФ) реализуются без нетривиальных вещественных умножений [15].
К сожалению, в силу соотношения (5), элемент w ° 2 (modfS) является корнем степени N = 2В (modf), что ограничивает максимальную длину преобразования Ферма, реализуемого без умножений, числом N = 32.
Кроме того, для “безошибочного” вычисления свертки спектральным методом с использованием ТЧП, число p должно быть достаточно велико [16]. В частности, при решении типичной для цифровой обработки изображений задачи вычисления двумерной свертки двух целочисленных массивов размера (512x512) с динамическими диапазонами 0 255, для числаp должно выполняться
неравенство:
p > (512)2 (256)2 = 234 > f4 = 216 +1.
Это существенно ограничивает возможности применения ТЧП Ферма в задачах обработки многомерной цифровой информации. Использование в качестве модулей в преобразовании (1) составных чисел Ферма доставляет серьезные трудности, связанные с существованием в модулярных кольцах по составным модулям делителей нуля и, как
следствие, с необратимостью некоторых элементов соответствующих колец и/или с не-ортогональностью базисных функций преобразования (1). Кроме того, даже при условии частичной компенсации отмеченных проблем, например, при распараллеливании вычислений в системе остаточных классов [11, 12], характерные преимущества “битовой” реализации арифметических операций именно в полях по модулям чисел Мерсен-на и Ферма не наследуются для вычислений в полях по модулям целых делителей составных чисел Мерсенна или Ферма. Действительно,
f5 = 232+1 = 641 -6 700 417
и сомножители уже не являются числами Ферма.
Альтернативное разложение чисел Ферма
Пусть, для определенности, число Ферма f имеет вид
f = 2B +1, B = 2r = 3t +1 (8)
и обладает свойствами:
(a) при всех 1 < s < 3t +1 числа f и 2s -1 взаимно просты;
(b) элемент 2(31 +1) обратим в фактор-
кольце
(c) число f не делится на 3, то есть f Ф 0 (mod 3).
Рассмотрим кольцо S целых элементов поля F разложения полинома j(z) = z3 + 2 над Q. В кольце S з Z для числа f наряду с обычным представлением составного числа в виде произведения целых рациональных чисел возможно представление в форме
f = 2 ^2 +1)2 g 32 + 1)(2t g 32 +1). (9)
где g - примитивный корень третьей степени из единицы.
Лемма 3.1. Элементы fF f2,f е S:
/ =(2' \[2 +l) / =(2' gIII +1) / =(2' g 1/2 +1)
попарно взаимно просты в кольце S.
Доказательство. Допустим, что существуют a, b b2 е S, не являющиеся единицами кольца S, такие что f = a b f = ab2. Нор-
мальным полем для многочлена j(z) является поле F=Q(y, 32 ). Пусть NormYY (x) есть относительная норма элемента x є F в поле Q(Y). Tрехэлементная группа Галуа полинома j(z) над подполем Q(g) циклична и действует тождественно на Q. Относительная норма элемента
z = х + У 32 є S равна Normgg (z) =x3+2y3. Поэтому из очевидных равенств
= Yf + 1 - Y, f2 - /і = (/і -1)(Y - 1)
следует
Normgg (/І -1) Normgg (g - 1) =
= NormgY (a) Normgg (b1 - b2 ),
2B (YY - 1)3 =
= NormgY (a) Normgg (b1 - b2). (10)
Tак как NormgY (a) не может быть четным числом, что противоречило бы равенству
NormgY f) = NormgY (ab1) =
=f = 2B +1 = Normgg (a) Normgg (b1),
то равенство (10) может выполняться только при условии делимости
NormgY (a) || (g - 1)3.(11)
Tак как Normgg (a) є Q(g), то вычисляя норму элементов (11) над Q, получаем:
Norm (NormgY (a)) II Norm (g - 1)3 = 27,
что противоречит условиюf Ф 0 (mod 3). Аналогично доказывается взаимная простота остальных элементов в условии леммы.
Пусть элементы f1, f2, f3 определены в лемме 3.1. Введем для удобства новые обозначения: р = f1, Q=f2, R =f3 .
Лемма 3.1 гарантирует возможность
S /
представления фактор-кольца yf в виде прямой суммы
S/ ° S/ е S/ е S/
/f /p /q /г
фактор-колец
где f = (f),
p /q /г-
p =(P), q =(Q), r =(R) - главные идеалы, порожденные элементами f, P, Q, R, и возмож-
ность вложения подкольца
f в эту
прямую сумму. Поэтому следующая лемма является частным случаем китайской теоремы об остатках [14]. Её доказательство приводится только для явного описания такого вложения.
7/
Лемма 3.2. Для любого W е /у7 существуют эффективно определяемые элемен-
ты X, Y, Z є Z
(/)
и константы a, b, c є
(13)
Y1 +(-Y3 ) 2t+1 + Y2 22t+1 = (3b) -1 W,
Zj +(-Z3 ) 2t+1 + Z2 22M = (3c) -1 W. (14)
Из соотношений (14) эффективно определяются значения X, Y, Z. Действительно, пусть с есть наименьший неотрицательный вычет для (3а) -1 W (mod f). Пусть далее quot (u // v) и exc (u // v) - неполное частное и остаток от деления числа u на v, соответственно. Тогда, например, для X имеем:
хf Z такие, что:
(а) справедливо равенство
W = aXQR + bYPR + cZPQ, (12)
причем
X = W (mod (P)), Y = W (mod (Q)),
Z = W (mod (R));
(б) числа X, Y Z допускают представления в форме:
X1 = exc (c // 2t+1 ), X2 = quot (c - X1 // 22t+1 ), (-X3 ) = (c - Xi) 2
-t -1
X2 2t.
X = X1 + X2 32 + X3 ^4;
Y = Y1 + Y24l + Y3^4;
Z = Z1 + Z 2Ш + Z 3^4;
0 < X1, Y1, Z1 < 2*+';
0 < X2, Y2, Z2, X3, Y3, Z3 < 2t .
Доказательство. Соотношение (12) является следствием китайской теоремы об остатках. Для доказательства свойства (а) не-
Z
обходимо доказать, что a,b,c е /f z . Непосредственно проверяются равенства:
RQ = (P - 3)P + 3, PR = (Q - 3)Q + 3,
PQ = (R - 3) +3.
Поэтому для выполнения равенств
aQR = 1 (mod P), bPR ° 1 (mod Q),
cPQ = 1 (mod R) достаточно положить
a ° 3-1 (mod f), b ° 3-1 (mod f),
c ° 3-1 (mod f).
Опуская рутинные выкладки, связанные с применением метода неопределенных коэффициентов, приведем выражения для X, Y., Z. (j = 1, 2 ,3 ) в форме:
X, +(-X3 ) 2t+1 + X2 22t+1 = (3a) -1 W,
Нетрудно заметить, что элементы X, Х3 допускают ^-битовое представление, а элемент Х1 допускает (£+1)-битовое представление.
Рассмотрим первое из равенств (13). Пусть
X! = +... + X020, Х2 = Хг2_12г-1 +...
+... + X0220, Х3 = Хг3_12г_1 +... + X0320
есть битовые представления элементов Х1, X, X. Тогда соотношение (13) для X можно переписать в виде:
X = X,1 ((2 ) + X,2-, ((2 ) _‘ + X,3-, ((2 ) _2 +
X,-, ((2 Г +... + X0 ((2). (15)
+
Равенство (15) можно интерпретировать как представление элемента X “в системе
счисления с основанием (^2)”, равенства,
аналогичные (15), для Y и 2 - как представления элементов Y и 2 “в системах счисления с
основанием
уЩ) и у^2) , соответственно.
Замечание 3.1. Несмотря на не вполне привычную терминологию (“системы счисления с иррациональным основанием”), никаких “приближенных” вычислений в данной работе не производится. Элементы, обозначенные (32), у(32) и у(32), есть просто три различных корня уравнения = 1 в фактор
О /
кольце у^ . Как будет показано ниже, такая
(неформальная) интерпретация равенства (15) позволяет ввести простые правила преобразований ассоциированных кодов
< х >= (X1, X2-!, X3-!, х-,..., 4). (16)
Отметим также, что существует необходимый математический формализм для придания корректности понятию “системы счисления с иррациональным основанием”. Tакие системы счисления достаточно давно и успешно применяются в информатике [17, 18].
Арифметические действия над элемен-S
тами кольца Sf индуцируют правила преобразования кодов (1б): при сложении элементов коды преобразуются по правилу “перенос в старший разряд через две позиции”;
умножение элемента X на (^2 ) равносильно циклическому сдвигу кода с инвертированием знака младшего бита и т.д. Как и в случае обычной арифметики кольца вычетов (mod /), в результате таких преобразований получаются не обязательно битовые векторы, которые, тем не менее, могут быть легко преобразованы к битовому виду с использованием соотношения (б). Tакие векторы, преобразованные по правилу (б) к битовому представлению, будем также называть редуцированными кодами компонентов элемента
< X(m) >
< X(n) >
кольца Sf и обозначать <X>*
<Y>* <Z>*.
(18)
N_1
:ІЗ'
n=0
где m = 0,1,..., N-1.
Замечание 4.1. При интерпретации редуцированных кодов как битовых представлений элементов поля классов вычетов по модулю f простого числа Ферма введенное шифт-преобразование совпадает с TЧП Ферма при w ° 2 (modf). В этом случае при m Ф к (mod 2B) равенство
T_1
0 =ХЗ-
n=0
r_mk < 3mn < е >*>*
(19)
Шифт-преобразования Ферма
Определение 4.1. Оператор 3, определенный на множестве (к+1)-мерных редуцированных кодов, такой что
3: <■¥>* = X Х-.....Х,, X) ®
<3 Х>- = (Х-,,...Х,,Х„ -X,)* (17)
будем называть оператором левого сдвига Ферма. Аналогично определяется оператор правого сдвига Ферма, являющийся обратным к 3.
Следующие два утверждения, сформулированные как леммы, приводятся без доказательств.
Лемма 4.1. Для любого (к+1)-мерного битового вектора X период последовательности <3П Х>* является делителем числа 2(к+1).
Определение 4.2. Пусть <Х(п)>* есть Ы-периодическая последовательность (к+1)-мерных редуцированных кодов. Определим шифт-преобразование Ферма соотношением
равносильно очевидному равенству 2B-1
0 ° £2”(”-‘> (modf). (20)
n=0
При m ° к (mod 2В) значение сумм (19) и (20) равны 2В (mod f). В общем случае шифт-преобразований значение суммы (19) определяется конкретной интерпретацией действий над кодами, связанной с арифметическими операциями в ассоциированном конечном кольце.
Лемма 4.2. Пусть <Е >* = (0, 0,..., 0, 1), T= 2В есть период последовательности
<3n Е >*. Тогда для колец K = ^ ^ , ^/r
при m Ф к (mod 2В) справедливо соотношение (19). При m ° к (mod 2В) значение суммы (19) равно 2В (mod p), (mod q), (mod r), соответственно.
Доказательство. Если рассматривать соотношение (17) как преобразование, индуцированное умножениями на элементы (^2 ),
у\у 2) и у^2; соответственно, то единственной причиной нарушения равенства (19) мос/
жет служить, например, для кольца необратимость элементов
г к _ 1 = (2) _ 1 при суммировании членов геометрической прогрессии. Этого не
происходит, если NormgY число, взаимно простое с f.
есть
Но при t ф 0 (mod 3) имеем
Norm-g х
х
(2 )-1 ) = [ ((2 )-1 I -((2 )-1 I-((2 )-1 | = 2Т -1.
При t ° 0 (mod 3) имеем
Norm YY
2 -1 =
У3 -1
чЗ
и существование нетривиального общего де-
лителя чисел Norm-
V 32 J-
невоз-
можно при и =0, І, 2 и t > 3. При t = І, 2, З
элементы
V 312 J-
являются единица-
Z/
кольце f) ;
(f) число f не делится на 3, то есть f ф О (mod 3),
то достаточно рассмотреть разложение числа f на множители в кольце S целых элементов поля разложения полинома y(z) = 2z3 + І над Q.
В S для числaf наряду с представлением в виде произведения целых рациональных чисел возможно представление в форме, аналогичной (9):
f =[ 2' Ш + ф'-
^+1][2'- з/Х
+1
ми кольца 8 и, следовательно, обратимы в соответствующих фактор-кольцах.
Лемма 4.2 позволяет рассматривать “правое” шифт-преобразование как обратное по отношению к “левому” шифт преобразованию (18) и наоборот.
Вычисление свертки
Рассмотренные выше шифт-преобразо-вания позволяют вычислять циклическую свертку целочисленных 2В-периодических последовательностей с помощью (В+1) -битовых процессоров по обычной спектральной схеме (см., например, [16]). Опуская детали описания такой схемы, сформулируем окончательный результат.
Теорема. Если для числа Ферма (8) выполняются условия (а)-(с),то для вычисления циклической свертки длины N = 2(31 + 1) достаточно выполнения:
1. девяти (шести левых и трех обратных) шифт-преобразований;
2. вычисления произведений компонентов спектров шифт-преобразований;
3. реконструкции значений свертки по китайской теореме об остатках в форме (12).
Если число Ферма имеет вид
/ = 2В +1, В = 2Г = Ъг -1 и обладает свойствами:
(ё) при всех 1 < ^ < 3t -1 числа/и 2* -1 взаимно просты;
(е) элемент 2(31 -1) обратим в фактор-
где у - примитивный корень третьей степени из единицы. Как и в лемме 3.1, доказывается взаимная простота элементов
2'-ш+Л (2'-з/Х+1
Доказательства аналогов соответствующих лемм и теоремы также существенно не отличаются от приведенных выше.
Непосредственная численная проверка показывает, что числа Ферма f , f , f7 удовлетворяют условиям (a)-(c) или (d)-(f), что позволяет вычислять свертки длин 64, 128, 256 с помощью описанного метода.
Заключение
Возможности рассмотренного метода, по мнению автора, не ограничиваются задачей «безошибочного» вычисления свертки. Представляется перспективным его использование, например, при быстром параллельном возведении в степень в конечных полях (массовая криптографическая задача), при сжатии информации и т.п.
В этой связи основную техническую трудность при обобщении на случай иррациональностей высших порядков представляет отыскание явного вида вложения кольца вычетов (modf) в прямую сумму колец (лемма 3.2).
Экстраполяция результатов на случай модулейp более общего видаp = bk ± 1 также не вызывает принципиальных затруднений. В [19] приведены разложения чисел p на простые множители. Практическая целесообразность такого обобщения ограничивается ис-
І33
ключительно возможностями вычислительных средств, ориентированных на недвоичное представление информации [10].
СПИСОК ЛИТЕРАТУРЫ
1. Фараджиев Р.Г. Аналитические способы вычисления процессов в линейных последовательных машинах // Известия АН СССР. Техническая кибернетика. 1965. N 5.
2. Фараджиев Р.Г., Ципкин Я.З. Преобразование Лапласа-Галуа в теории последовательных машин // Доклады Академии наук СССР. 1966. Т.166. N 36.
3. Schoenhage A., Strassen V. Schnelle Multiplication grosser Zahlen // Computing. B.7. 1966. N.3/4.
4. Rader C.M. Discrete convolution via Mersenne transorm // IEEE Trans. Comp. C-21. 1972.
5. Rader C.M. On the application of the number theoretic methods of high-speed convolution to two-dimensional filtering // IEEE Trans. on Circuits and Systems. 1975. V.22.
6. AlfredsonL.-I. A fast Fermat number transform for long sequences // Proc. EUSIPCO-94, Edinburg, Scotland. 1994. V.111.
7. Alfredson L.-I. VLSI architectures and arithmetic operations with application to the Fermat number transform. Link»oping Studies in Sci. and Technology, Dissertation. 1996. N.425.
8. Boussakta S., Holt A.G.J. Calculation of the discrete Hartley transform via F ermat number
transform using VLSI chip // IEE Proc. 1988. V.135. N 3.
9. Towers P.J., Pajayakrit A., Holt A.G.J. Cascadable NMOS VLSI circuit for implementing a fast convolver using the Fermat number transform // IEE Proc. 1987. V.135. N 2.
10. Вариченко Л. В., Лабунец В. Г., Раков М. А. Абстрактные алгебраические системы и цифровая обработка сигналов. Киев: Наукова думка, 1986.
11. Дэвенпорт Дж., Сирэ И., Турнье Э. Компьютерная алгебра. М.: Мир, 1991.
12. Торгашев В.А. Система остаточных классов и надежность ЦВМ. М: Сов. радио, 1973.
13. Chernov V.M.,Pershina M.V «Error-free» calculation of the convolution using generalized Mersenne and Fermat transforms over algebraic fiels // Proc. CAIP’97. Springer. LNCS. 1296.
14.Ленг C. Алгебра. М.: Мир, 1965.
15. Блейхут Р. Быстрые алгоритмы цифровой обработки сигналов. М.: Мир, 1987
16. Нуссбаумер Г. Быстрое преобразование Фурье и алгоритмы вычисления сверток. М.: Радио и связь, 1985.
17. Кнут Д. Искусство программирования для ЭВМ. Т. 2. М.: Мир, 1977
18. Стахов А.П. Коды золотой пропорции. М.: Радио и связь, 1984.
19. Brillart J., Lehmer D.H., Selfridge J.L., Tuckerman B., Wagstaff S.S. Factorization of bk +1, b=2,3,5,6,7,10,11, 12 up high powers // Contemp.Math. AMS. 1988. V.22.
THE PARALLEL ALGORITHMS OF FOURIER-GALOIS TRANSFORMS SYNTHESIS IN DIRECT SUMS OF FINITE RINGS
© 2000 V.M. Chernov
Image Processing System Institute of Russian Academy of Sciences, Samara
The method of fast algorithms of discrete transforms synthesis is considered. The algorithms synthesized with this method are the algorithms of «error-free» calculations of integer-valued discrete convolutions. The approach proposed in the paper is based on the following fact. The factorization of the elements in algebraic rings is can be performed in a number of ways. The choice of factorization is determined by the way in which the modular calculations are organized.