Алгоритмы умножения карабуцы и Штрассена для разреженных структур Текст научной статьи по специальности «Математика»

ЛИТЕРАТУРА

1. Кнут Д.Э. Искусство программирования. Получисленные алгоритмы. М.: Издат. дом «Вильямс», 2001. Т. 2. С. 553.

2. Mulders Т., Storjohann A. On Short Multiplications and Divisions // Applicable Algebra in Engineering, Communication and Computiong. 2000. V. 1.

3. Jebelean T, Dragan M, Tepenue D, Negru V. Parallel algorithms for practical multiprecision arithmetic using the Karatsuba method //

Proceedings of Algoritmy. 2000. P. 1-10.

4. Jebelean T. Practical integer division with Karatsuba complexity // W. Kuechlin, ed. ISSAC’97. ACM Press, 1997. P. 339-341.

БЛАГОДАРНОСТИ: Работа частично финансировалась грантами Министерства образования РФ Е02-2.0-98 и «Научный потенциал» 23-03-24.

АЛГОРИТМЫ УМНОЖЕНИЯ КАРАЦУБЫ И ШТРАССЕНА ДЛЯ РАЗРЕЖЕННЫХ СТРУКТУР

© Г.И. Малашонок, Е.С. Сатина

Обсуждается применение быстрых алгоритмов умножения к разреженным структурам. Даются оценки для допустимой степени разреженности матриц и полиномов, при которой быстрые алгоритмы умножения имеют преимущество по мультипликативной сложности по сравнению со стандартными алгоритмами. Исследование ведется для алгоритмов Карацубы [1] и Штрассена [2] и предполагается равномерное распределение нулевых элементов.

Исследуемым параметром является коэффициент плотности р - отношение числа ненулевых коэффициентов к числу всех коэффициентов матрицы или полинома. При р = 1 все коэффициенты ненулевые и быстрые алгоритмы эффективнее стандартных. При уменьшении р эффективность быстрого алгоритма по отношению к стандартному алгоритму уменьшается. Задача состоит в указании граничного значения р, при котором сложность быстрого алгоритма становится равной сложности стандартного алгоритма

АЛГОРИТМ КАРАЦУБЫ

Пусть R - поле. Рассматривается задача умножения двух полиномов Дх) и g(x) степени п (п = 2n - 1) из кольца Л[х], у которых ненулевые элементы распределены равномерно и их плотность равна а. Применяется рекурсивный алгоритм:

fg = (А + Вх)(С + Dx) =

= АС + (АС + BD - (А - В)(С - D))x + BDx,

где А, В, С, D - многочлены степени 2ЛМ - 1. Чтобы учесть разреженность полиномов, будем полагать, что глубина рекурсии равна к. После чего вычисляется произведение полиномов степени 2'v к - 1 стандартным способом и при этом требуется (2N~ka)2 операций умножения.

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

сі

5=0

р-1

2J-

Для поля характеристики 0 достаточно перейти к пределу при р —> о0, и в выражении для ф'(а) исчезнет дробь.

Сравним полученную оценку сложности со сложностью стандартного умножения - (а2л)2 для разной глубины рекурсии (1, 2, ..., 15) и для разной плотности исходных многочленов. Полученные результаты представлены в таблице плотностей (табл. 1). В левом столбце указана глубина рекурсии, в верхней строке -отношение сложности алгоритма Карацубы к сложности стандартного умножения. На пересечении соответствующих столбца и строки - плотность полиномов-сомножителей.

Например, при глубине рекурсии 15 алгоритм Карацубы имеет вдвое меньшую мультипликативную сложность, чем стандартный алгоритм, когда плотности полиномов сомножителей равны 0,1503.

В практических целях можно пользоваться графиком, приведенным на рис. 1. На этом графике изображены линии уровня (1, 1/2, 1/4, 1/8, 1/16) с постоянным отношением сложности алгоритма Карацубы к сложности стандартного алгоритма. По горизонтальной оси откладываются логарифмы степеней полиномов сомножителей N = \0g2n, по вертикальной оси - плотности.

АЛГОРИТМ ШТРАССЕНА

Формула Штрассена [2] для умножения матриц дает алгоритм для вычисления произведения матриц порядка п = 2^, в котором ~ и2 8 операций умножения.

В левом столбце указана глубина рекурсии, в верхней строке - отношение сложности алгоритма Штрассена к сложности стандартного умножения. Число, стоящее на пересечении соответствующих столбца и строки - это плотность матриц сомножителей.

Таблица

Таблица 2

Плотности полиномов при заданной глубине рекурсии (к) и заданном отношении сложности алгоритма Карацубы и стандартного алгоритма

к 1 1/2 1/4

1 0,5858

2 0,5236

3 0,4668 0,8819

4 0,4150 0,7344

5 0,3680 0,6237 0,9702

6 0,3256 0,5353 0,8247

7 0,2875 0,4621 0,7054

8 0,2534 0,4002 0,6061

9 0,2229 0,3473 0,5224

10 0,1959 0,3018 0,4512

11 0,1719 0,2624 0,3903

12 0,1507 0,2283 0,3379

13 0,1320 0,1985 0,2928

14 0,1155 0,1727 0,2539

15 0,1009 0,1503 0,2202

0.8

0.6

0.4

0.2

Рис. 1.

к 1/8 1/16 1/32

7 1,0

8 0,8895

9 0,7654 1,0

10 0,6597 0,9482

11 0,5694 0,8186

12 0,4920 0,7072 1,0

13 0,4254 0,6112 0,8709

14 0,3681 0,5285 0,7532

15 0,3187 0,4572 0,6516

ч \ -

ч \ \

ч \ \ ■ V \

X \

ч

10

Если плотность матриц более 2/3, то алгоритм Штрассена при любой глубине рекурсии предпочтительнее, чем стандартный алгоритм. Если плотности сомножителей менее 0,39, а порядок матриц менее 1000, то алгоритм Штрассена проигрывает по сравнению со стандартным алгоритмом. В практических целях можно пользоваться графиком (рис. 2).

Плотности матриц при заданной глубине рекурсии (к) и заданном отношении сложности алгоритма Штрассена и стандартного алгоритма

к 1 1/2 1/4 1/8

1 2/3

2 0,6637 1,0

3 0,6345 0,9272

4 0,5956 0,8555

5 0,5564 0,7928

6 0,5187 0,7363

7 0,4833 0,6847 0,9690

8 0,4503 0,6374 0,9016

9 0,4196 0,5938 0,8398

10 0,3912 0,5534 0,7827

11 0,3649 0,5160 0,7298

12 0,3404 0,4813 0,6808 0,9627

13 0,3176 0,4491 0,6352 0,8983

14 0,2964 0,4192 0,5928 0,8384

15 0,2767 0,3913 0,5535 0,7826

0.8

0.6

0.4

0.2

Рис. 2.

10

Чтобы учесть разреженность матриц будем полагать, что вычисления производятся по рекурсивному алгоритму Штрассена и глубина рекурсии равна к. После чего вычисляются произведения матриц порядка 2м~к стандартным способом, при этом для произведения матриц порядка п = 2ы~к, имеющих плотности а, требуется г?а1 операций умножения.

Тогда общее число операций умножения составит:

/=0 j=0\J

]

І +І-

:У+]~к(2^к)Ц1(а)^(а),

где функции ф'(а) были определены выше.

Сравним полученную оценку сложности со сложностью стандартного умножения матриц (а2п3) для разной глубины рекурсии (1, 2, ..., 15) и для разных плотностей исходных матриц. Полученные результаты приведены в таблице плотностей.

На этом графике изображены линии уровня (1, 1/2, 1/4) с постоянным отношением сложности алгоритма Штрассена к сложности стандартного алгоритма. По горизонтальной оси откладывается логарифм

порядка матрицы N = ^2«, по вертикальной оси -плотность.

Отметим, что алгоритм Винограда [3], в котором столько же умножений, но меньше сложений проигрывает алгоритму Штрассена для разреженных матриц. Каждая рекурсия в этом алгоритме требует суммирования четырех подматриц, а не двух, как в алгоритме Штрассена, поэтому плотность операндов растет быстрее.

ЛИТЕРАТУРА

1. Strassen V. Gaussian Elimination is not optimal // Numerische Mathematik. 1969. V. 13. P. 354-356.

2. Карацуба A. II ДАН СССР. 1962. Т. 145. С. 293-294.

3. Кнут Д.Э. Искусство программирования. Получисленные алгоритмы. М.: Издат. дом «Вильямс», 2001. Т. 2. С. 553.

БЛАГОДАРНОСТИ: Работа частично финансировалась грантами Министерства образования РФ Е02-2.0-98 и «Научный потенциал» 23-03-24.

ОБ ОЦЕНКЕ ДЕТЕРМИНАНТА ПОЛИНОМИАЛЬНОЙ МАТРИЦЫ © Г.И. Малашонок, П.А. Азарова

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

Пусть А = (а,и) - матрица порядка п, с коэффициентами а,, e 2 [х], max deg a,j = к и модуль максимального коэффициента полиномов а,у не превосходит а (V0 <

i,j < п). Тогда det(/i) - это полином степени не выше пк с целыми коэффициентами. Задача стоит в получении верхней оценки для коэффициентов этого полинома.

Для оценки полинома det А воспользуемся неравенством Адамара:

(det А)2 < П,„0" W а,/. (1)

Вспомогательной задачей является задача оценки коэффициентов bj полинома

Из равенства х = (У’+1 - \)/(х - 1), возведением в степень s, получим

aw ь, х1 )(£;=</ с‘ (-1 г )=w с; (-1 г ^+1) •

Приравняем коэффициенты при одинаковых степенях неизвестного и получим систему из s(p + 1) + 1 линейных уравнений с ps + 1 неизвестным bj (i = 0, .., ps)

bo — 1

ЬоСД-i)1 + г>,сД-1)° = о

U^(-\rcrb,=g,j<s.

Zi=jJ(-1 У~'Сi~' bi=gl,j>s.

bps-1

Здесь gj = СД-l/ при/ = к(р + 1), 0 < к < s, и gi = 0 для всех остальных j. При этом полагаем, что С/ = 0 при v < 0. Поэтому в каждом уравнении не более чем 5 + 1 коэффициент отличен от нуля, и матрица коэффициентов имеет нижний треугольный вид с единичной главной диагональю. В силу симметрии Ьк = bps_k Поэтому достаточно найти первые г = [ps/2 + 1 ] неизвестных, то есть ограничиться системой из первых г уравнений.

Для решения такой системы требуется вычислять определители dk порядка к

dk = det (dij), dij = 0 приj> i+ 1 ,dij= 1 приj = i+ 1, dij = (-iy~'+1 (y~'+\ при j < i. Например,

| с, 1 0 0| d4 - | c2 Ci 1 0 |, где Cj = (-1)1 | C3 C2 C] 1 |

I C4C3C2CI I

Получим dn = C\d„A - C2d„ 2 +■■ (-l)'c,</„_, + .. + + (-1 )"c„, разложением по первому столбцу.

Следовательно, dt = C\, d2 = c2 - c2, d3 = C|3 - 2 Cj c2 +

+ С)ИТ. д.

Прямые вычисления дают:

dx=-s,d2 = (s2 + s)/2!,

d3 = -(s3 + 3s2 + 2!s)/3!,

rf4 = (s4 + 6s3+ lls2 + 3!s4)/4!,

ds = -(s5 + 10s4 + 35s3 + 50s2 + 4!s)/5!,

d6 = (s6 + 15s5 + 85s4 + 225s3 + 274s2 + 5!s)/6!,

d1 = -(s7 + 21s6 + 175s5 + 735s4 + 1624s3 + 1764s2 +

+ 6!s)/7!.

Так как диагональные элементы в матрице коэффициентов единичные, то решениями bj системы будут определители d’j, которые получаются дополнением определителя dj сверху единичной строкой (первый элемент - единица, остальные - нули) и справа -столбцом g’„ который образован первыми i элементами столбца g = (gj).

Отсюда следует, что b(i = 1; 6,- = (-1)' d, для 1 < i <р, b, = (-1)'' d, + 0т (-1)'+^+1> С/ dt_k(p+x) для р < I <

^ Vps!2\, т = \_И(р + 1)1