Об умножении скалярных и матричных многочленов с целыми коэффициентами Текст научной статьи по специальности «Математика»

УДК 519.61

А.С. Бусленко, Х.Д. Икрамов

ОБ УМНОЖЕНИИ СКАЛЯРНЫХ И МАТРИЧНЫХ МНОГОЧЛЕНОВ С ЦЕЛЫМИ КОЭФФИЦИЕНТАМИ

(кафедра общей математики факультета ВМиК, e-mail: [email protected])

Показано, что в системе Mathematica 5 точное произведение степенных рядов, коэффициенты которых являются целыми числами или целочисленными матрицами, можно вычислить значительно быстрее, если вместо библиотечной функции для этого произведения использовать методы теоретико-числовых преобразований, представляющие собой аналоги быстрого преобразования Фурье в кольце целых чисел Z.

1. Введение. Пусть заданы многочлены степени п - I

п-1

Р(х) =

г=0

И

и — 1

Q(x) = ^ ft®'

г=0

и требуется вычислить их произведение

2п-2 г

V(x) ■ Q(x) = Щх) =

г=0 j=О

Коэффициенты многочлена-произведения TZ{x) суть компоненты свертки векторов

РТ = (Ро Pi ••• Pn-i) , QT = (go <?i ••• Qn-i) ■

Прямой метод решения этой задачи требует порядка п2 умножений. Хорошо известно, что использование техники быстрого преобразования Фурье (БПФ) позволяет уменьшить число операций до 0(п logn). Применение этого алгоритма, реализованного над полем комплексных чисел С, к перемножению многочленов рассмотрено в [1]. Однако стандартное БПФ, выполняемое в комплексной арифметике, не позволяет получить точное произведение многочленов с целыми коэффициентами. Такое произведение может быть найдено обычным алгоритмом умножения за 0(п2) операций с целыми числами. В то же время существуют аналоги дискретного преобразования Фурье, называемые методами теоретико-числовых преобразований, которые работают над Ъм и позволяют обеспечить такую же экономию вычислений, как БПФ. В разделах 2 и 3 напоминаются теоретико-числовые факты, лежащие в основе этих методов. Опыт их применения к перемножению многочленов и полиномиальных матриц с целыми коэффициентами обсуждается в разделах 4 и 5. Приведено сравнение эффективности таких методов с эффективностью библиотечных процедур системы Mathematica 5 для перемножения многочленов и полиномиальных матриц.

2. Теоретике-числовые преобразования. Пусть а и т — взаимно простые натуральные числа, причем т ^ 1. Рассмотрим натуральные степени числа а:

2 f СЬу Qi ► » » » ► Qi »

Наименьшее /, такое, что

= 1 mod m,

называется экспонентой а по модулю т (или порядком а по отношению к т):

f = ехрто а.

Если m ) 1, то число ф(т) натуральных чисел, не превосходящих т и взаимно простых с т, называется функцией Эйлера. Известно, что

к

ф(т) = т Д(1 - l/pi), i= 1

к

где т = П рГ есть разложение т на простые множители. Если р простое иг) 1, то ф(рг) = рг —рг~1.

г= 1

Если ехртош = ф(т), то ш называется первообразным корнем из единицы по модулю т. Если шп = 1 mod m, то п должно делить ф(т).

Дискретное преобразование Фурье, в котором вместо комплексных экспонент используются целые числа, являющиеся первообразными корнями из единицы, а операции выполняются по модулю целого числа М (простого или составного), известно как теоретико-числовое преобразование, или NTT (Number Theoretic Transform). Приводимые ниже утверждения описывают свойства таких преобразований. Формулировки этих утверждений цитируются по книге [2], а доказательства могут быть найдены в [3].

Результат применения теоретико-числового преобразования к целочисленному вектору w размерности п будем обозначать NTT(*w) и называть число п длиной преобразования. Говорят, что теоретико-числовое преобразование обладает свойством циклической свертки, если

NTT (и *v) = NTT (и) • NTT(v)

для любых целочисленных векторов и и v. Символ * обозначает операцию свертки, а символ • — покомпонентное умножение двух векторов.

Утверждение 1. Теоретико-числовое преобразование длины п с первообразным корнем ш, определенное по модулю простого р, обладает свойством циклической свертки тогда и только тогда, когда ш — корень из единицы по модулю р порядка п. При этом п делит (р — 1).

Следующие два утверждения характеризуют длину п теоретико-числового преобразования для составного модуля М. При заданных М шп теоретико-числовое преобразование полностью определено, если выбран корень порядка п.

Утверждение 2. Теоретико-числовое преобразование длины п с корнем ш, определенное по модулю составного Л /. обладает свойством циклической свертки тогда и только тогда, когда:

1) шп = 1 mod М;

2) существует целое число п-1, обратное к п по модулю М:

п ■ п~1 = 1 mod Л/;

3) gcd(wd — 1 ,М) = 1 для всякого d, такого, что n/d, простое.

Замечание. В случае М = рг, гдер — простое число, эта теорема означает, что п делит р—1. Действительно, п-1 существует тогда и только тогда, когда п и М взаимно простые; п должно делить ф(М) = рг~1(р — 1), следовательно, п делит (р — 1).

к

Утверждение 3. Для М = П Р? теоретико-числовое преобразование длины п, определенное

г= 1

по модулю Л /. обладает свойством циклической свертки тогда и только тогда, когда

п\ gcd{(pi - 1), (р2 - 1),-.., (рл - !)}• Таким образом, в этом случае максимальная длина преобразования есть gcd{(pi — 1),..., (рн — !)}•

3. Преобразования Мерсенна, Ферма и Радера. Чтобы теоретико-числовые преобразования имели преимущество перед прямым вычислением свертки, их нужно уметь выполнять эффективно. Если п есть произведение большого числа простых сомножителей (наиболее выгоден случай, когда п — это степень двойки), то возможно применение быстрых алгоритмов, аналогичных алгоритму БПФ. Для практических целей длина п должна быть достаточно большой, а арифметика по модулю М — достаточно простой. В частности, несложным должно быть умножение на степени числа ш. Например, если ш есть 2 или степень двойки, то такие умножения сводятся к сдвигам разрядов второго сомножителя.

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

Если М четное, то оно имеет делитель 2; поэтому, согласно утверждению 3, птах = 1. Таким образом, для практического использования М должно быть нечетным числом. Если М простое, то максимальная возможная длина преобразования птах = М — 1. Числа Л/ 2>' - I. где р нечетное простое, называют числами Мерсенна, а теоретико-числовые преобразования по модулю таких М — преобразованиями Мерсенна. Для таких М число ш = 2 является корнем порядка р, что дает длину преобразования п = р. Выбор ш = —2 дает преобразование длины п = 2р. Достоинство преобразований Мерсенна состоит в том, что умножение на степени числа ш заменяется операциями сдвига. С другой стороны, имеются и серьезные недостатки — отсутствие быстрой модификации и жесткая связь между длиной слова М и длиной преобразования п.

Пусть Л/ /<', 21' + 1 = 2Ь + 1, где Ь = 2г ш I — неотрицательное целое число. Такие М называют числами Ферма. Первые пять из этих чисел (соответствующие значениям £ = 0,1,2,3,4) просты; большие числа Ферма могут быть и составными. Теоретико-числовые преобразования по модулю чисел Ферма называют преобразованиями Ферма. Поскольку числа Ферма М = 2Ь + 1 вплоть до ¿<4 простые, то максимальная длина преобразования птах = 2Ь; возможно также преобразование любой длины п = 2™, т ^ Ь. Для чисел Ферма максимальная длина реализуется выбором ш = 3. При выборе ш = 2 порядок равен п = 2Ь = 2*+1, т.е. длина преобразования в этом случае меньше. Однако, поскольку п является степенью двойки, может быть применен быстрый алгоритм вычисления теоретико-числового преобразования, аналогичный алгоритму БПФ. Алгоритмы этого типа называют преобразованиями Радера.

4. Умножение скалярных многочленов с целыми коэффициентами. Поставленная в пункте 1 задача точного и эффективного умножения двух многочленов с целыми коэффициентами от одной переменной решается следующим образом: каждый из векторов длины п, составленных из коэффициентов исходных многочленов,

РТ = (ра Р1 ■■■ рп-1) , ЯТ = (до 91 ••• Чп-1) дополняется справа п нулями, вычисляются их теоретико-числовые преобразования длины 2п:

Р = ШТ((ро Р1 ... рп_г 0 ... 0)), д = ЭТТ((д0 91 ... Яп-1 0 ... 0)),

компоненты полученных векторов попарно перемножаются:

Т1к = Рк- (2к, А; = 0, 2п - 1, после чего к вектору Л применяется обратное преобразование длины 2п:

г = ЫуегзеЫТТ (11). Вектор г составлен из коэффициентов многочлена-произведения

2п-2

Щх) = ^ Пх\

г=0

Оказалось, что алгоритм теоретико-числового преобразования не реализован в компьютерно-алгебраической системе МаШетаМса 5. Поэтому мы составили собственную программу VI п11 \\'п Ы'ХТТ для описанного выше метода перемножения многочленов. В табл. 1 показаны результаты сравнения этой программы с библиотечной функцией умножения многочленов (оператор Ехрапс1[р1 * р2}) в МаШетаМса 5. Рассматривались скалярные многочлены с целыми коэффициентами степени п — 1 для п = 1024, 2048,4096. Прямые и обратные теоретико-числовые преобразования длины 2п проводились по модулю чисел Ферма Л/ 2" + I с корнем ш = 2. Длина исходных векторов п выбиралась равной степени двойки, поэтому применялся быстрый алгоритм, аналогичный алгоритму БПФ. Чтобы результаты не выходили за пределы кольца вычетов по модулю Л/. коэффициенты многочленов Т(х)

и Q(x) должны соотноситься с М определенным образом, а именно: коэффициенты г^ произведения TZ должны удовлетворять условию

-М/2 <п< М/2. Если \pi\ ^ а и \qi\ ^ /3, то |rj|max ^ п \aß\ < М/2, или

М > 2n\aß\ .

Вместо и п2/2 скалярных умножений в случае прямого перемножения двух многочленов степени n — 1 алгоритм, реализованный в MultWithFNTT, требует порядка п log2 п + п умножений (три быстрых теоретико-числовых преобразования и покомпонентное умножение двух векторов длины 2п). Из табл. 1 видно, что уже при п = 4096 MultWithFNTT вычисляет произведение за существенно меньшее время, чем библиотечная функция. (Все значения времени в обеих таблицах указаны в секундах. Вычисления проводились на компьютере со следующими характеристиками: процессор Intel Pentium 4, частота 1,80 ГГц, 512 Мб оперативной памяти, кэш данных первого уровня 8 Кб, кэш второго уровня 256 Кб.)

Таблица 1

Скалярные многочлены

п Процедура MultWithFNTT Expand [pi * р2]

1024 11,141 11,625

2048 72,765 85,547

4096 440,468 767,985

5. Умножение матричных многочленов с целыми коэффициентами. Предположим, что вместо скалярных многочленов имеем два матричных многочлена от одной переменной:

п-1

А(х) = xiAi i=0

и

п-1

в(х) = J2xiBi-> i=0

где

AuBi(EZdxd, г = 0, п — 1. Требуется вычислить матричный многочлен

2п-2

С(х) = А(х)В(х) = Y^ c'i е zdXd' * = 0,2п — 2.

г=0

Эту задачу можно свести к задаче перемножения скалярных многочленов. Рассмотрим А(х) и В(х) как dxd-матрицы, элементами которых являются многочлены (п—1)-й степени с целыми коэффициентами. Тогда для элементов матрицы-произведения С(х) имеем

d d crs(x) = (C(x))rs = ^2(A(x))rt(B(x))ts = ^ art(x)bts(x), r, s = 1, d. t=l t=l

Для вычисления каждого элемента crs нужно найти d произведений элементов исходных матриц; для каждого из этих произведений применяется описанная выше техника перемножения скалярных многочленов.

Описанный метод можно дополнительно ускорить. Каждый из элементов исходных матриц участвует в вычислении d элементов матрицы-результата. Поэтому для экономии вычислений выгодно, единственный раз применив теоретико-числовое преобразование к вектору коэффициентов каждого многочлена art(x) или bts(x), сохранить полученные векторы в памяти компьютера и использовать их в дальнейшем при вычислении сверток. Аналогичный подход был использован в [1] при перемножении степенных рядов, коэффициентами которых являются вещественные или комплексные матрицы.

Численные эксперименты проводились для матричных многочленов степени п — 1 с коэффициентами из Zdxd, где d = 2, 3, а п = 1024, 2048. Для перемножения этих многочленов, с одной стороны, использовалась составленная нами процедура FastMultMatrWithFNTT, а с другой — библиотечная функция перемножения матриц Expand\pl.p2]. Из табл. 2 видно, что, как и в случае скалярных многочленов, использование методов, основанных на применении теоретико-числовых преобразований, позволяет значительно сократить время вычислений при больших п.

Таблица 2

Матричные многочлены

п Порядок d Процедура FastMultMatrWithFNTT Expand \pl.p2]

1024 2 93,547 111,703

1024 3 306,437 353,766

2048 2 540,391 791,125

СПИСОК ЛИТЕРАТУРЫ

1. Буеленко А. С., Икрамов X. Д. Об умножении числовых и матричных степенных рядов // ЖВМиМФ. 2005. 45. № 1. С. 3-7.

2. Krishnamurthy Е. V. Error-free polynomial matrix computations. N. Y.: Springer-Verlag, 1985.

3. Nussbaumer H.J. Fast Fourier transform and convolution algorithm. N.Y.: Springer-Verlag, 1980.

4. Pan V. Y. Structured matrices and polynomials. Boston: Birkhauser, 2001.

5. Каханер Д., Моулер К., Нэш С. Численные методы и программное обеспечение. М.: Мир, 1998.

Поступила в редакцию 10.04.07

УДК 519.233.32, 519.257 Р.Е. Игла

ПАРАМЕТРИЧЕСКИЙ АНАЛИЗ ДАННЫХ ПРОТОЧНОЙ ЦИТОМЕТРИИ

(кафедра математической статистики факультета ВМиК, e-mail: [email protected])

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

1. Введение. Проточная цитометрия представляет собой технику для быстрого оптического анализа отдельно взятых клеток.

Принципы проточной цитометрии весьма просты. Клетки или ядра поодиночке пересекают сфокусированный световой пучок. Свет определенной длины возбуждает молекулы флуоресцирующих красителей, связанных с различными клеточными компонентами. Свет, испускаемый красителями, собирают с помощью системы линз и зеркал и разлагают на компоненты. Световые сигналы детектируют, преобразуют в электрические импульсы и далее в форму, удобную для компьютерной обработки и хранения информации. В результате процедуры, описанной выше, для каждой частицы получается несколько параметров. Анализируя эти параметры, можно получить некоторую общую информацию о совокупности изучаемых частиц [1].

Одной из важнейших задач проточной цитометрии является определение доли клеток, обладающих определенными свойствами, например содержащих некоторые интересующие нас гены. В данной статье будет рассмотрена одна из таких задач: найти процент CD3-положительных (CD3+), CD8-отрицательных (CD8-) лейкоцитов, содержащих INF7 (интерферон-гамма) или IL-4 (интерлейкин-4) среди всех ( 'l).'l' -. СБ8~-лейкоцитов.