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

Информатика

УДК 519.677

ВЫЧИСЛЕНИЕ ГИПЕРГЕОМЕТРИЧЕСКИХ РЯДОВ С КВАЗИЛИНЕЙНОЙ ВРЕМЕННОЙ И ЛИНЕЙНОЙ ЁМКОСТНОЙ СЛОЖНОСТЬЮ

С. В. Яхонтов

Санкт-Петербургский государственный университет, математико-механический факультет,

198504, Санкт-Петербург, Старый Петергоф, Университетский пр-т, 28.

E-mails: sergey_home_mail@inbox. ru

Проводится построение простого для практической реализации алгоритма со сложностью 0(М(те) log (те)2) по времени и О (те) по памяти для вычисления ги-пергеометрических рядов с рациональными коэффициентами на машине Шёнха-ге, где М(те) — сложность умножения целых чисел. Показывается, что данный алгоритм пригоден в практической информатике для. построения конструктивных аналогов часто используемых констант математического анализа.

Ключевые слова: конструктивные вещественные числа, гипергеометрические ряды, квазилинейная временная сложность, линейная ёмкостная сложность.

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

С помощью Sch(FQLINTIME//LINSPACE) будем обозначать класс алгоритмов, вычислимых на машине Шёнхаге, квазилинейных по времени и линейных по памяти. Основной особенностью машины Шёнхаге является возможность рекурсивных вызовов процедур. Под квазилинейностью понимается ограниченность сверху функцией вида 0(nlog(n)fc) при некотором к.

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

Известно [2], что линейно сходящиеся гипергеометрические ряды с рациональными коэффициентами можно вычислять с помощью метода двоичного деления с временной сложностью 0(M(n)(log(n))2) и ёмкостной слож-

Яхонтов Сергей Викторович (к.ф.-м.н.), доцент, каф. информатики.

ностью 0(п1(щ(п)) (М(п)—сложность умножения п-битовых целых чисел). В последнее время появились работы, например [3], в которых описываются алгоритмы (на основе модифицированного метода двоичного деления) вычисления значений линейно сходящихся гипергеометрических рядов с временной сложностью 0(М(п)(к^(п))2) и ёмкостной сложностью О(п). В [4] предложены простые для практической реализации алгоритмы для вычисления степенных рядов с полиномиальной временной и линейной ёмкостной сложностью.

Недостатком метода из [3] является относительная сложность практической реализации. В настоящей работе предлагается простой для практической реализации алгоритм, который, так же как и алгоритм из [3], является квазилинейным по времени и линейным по памяти. Идея работать с некоторой заданной точностью для того, чтобы вычислить значение гипергеомет-рического ряда с линейной ёмкостной сложностью, предложена в [5].

Вычислительная сложность конструктивных вещественных чисел и функций подробно освещается в монографии [6]. Множество конструктивных вещественных чисел с квазилинейной по времени и линейной по памяти сложностью вычисления двоично-рациональных приближений будем обозначать Бс}г(РС2ЫМТ1МЕ//ЫИБРАСЕ)ср (данное множество можно определить по аналогии с множеством РЫКБРАСЕ вычислимых действительных чисел [4] с учётом другой вычислительной модели).

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

1. Метод двоичного деления. Данный метод предназначен для вычисления значений линейно сходящихся рядов с рациональными коэффициентами, в частности, для вычисления гипергеометрических рядов вида

где а, Ь, р, д — полиномы с целыми коэффициентами; данный ряд линейно сходится, если его частичная сумма

где !л{к) — линейная функция от к, отличается от точного значения не более

В классическом варианте метод двоичного деления состоит в следующем. Обозначим к\ = /л(к). Рассмотрим частичную сумму (2) для некоторых чисел Ч и г2, 0 ^ п ^ к\, 0^г2 ^ к\, 1\ ^ г2:

Будем ВЫЧИСЛЯТЬ ВеЛИЧИНЫ Р(%1,12) = Р{Ч) ■ ■ ■ Р(,Ь), Qi.il, ^2) = Я(,Ч) ■ ■ •

В(г 1, гг) = Ь(п)... Ь(г2) и Т(гь г2) = Б(гь г2)<3(*ъ г2)5'(*1, г2). Если гх = г2, то

(1)

(2)

чем на 2 к: |5 — 3(/л(к))\ ^ 2 к.

(3)

величины вычисляются напрямую. Иначе ряд делится на две части, левую и правую, И величины вычисляются для каждой из

частей рекурсивно. Затем полученные величины комбинируются:

P{iii ^2) — PiPr, Q(i l) ^2) — QiQr, ^2) — BiBr,

T(i\,i2) = BrQrTi + BiPxTr.

Алгоритм начинает свою работу с i\ = 0, 12 = k\. После вычисления величин Т(0, fci), В{0, fci), Q(0, k\) осуществляется деление Т(0, fci) на Б(0, fci)Q(0, fci), чтобы получить результат с заданной точностью. Длины чисел Т(0,к\) и B(0,ki)Q(0,k\) пропорциональны k\og(k), то есть алгоритм двоичного деления является квазилинейным по памяти; временная сложность данного алгоритма — 0(М(к) log(fc)2) [2].

2. Алгоритм вычисления гипергеометрических рядов из [3]. В алгоритме из [3] вычисляются редуцированные числитель и знаменатель дроби Т/Q, где Т = TBQ, Q = BQ, имеющие длину О(п), тогда как в классическом варианте вычисляются числитель и знаменатель длины O(nlog(n)). Затем производится деление редуцированного числителя на редуцированный знаменатель с тем, чтобы получить значение суммы ряда с заданной точностью. Параметром усовершенствованного алгоритма из [3] является модуль m — верхняя граница длины редуцированных Т и Q; модуль m вычисляется как произведение достаточного числа простых чисел. Ряд разбивается на N/G групп, где С— величина, зависящая от т, и для каждой группы применяется классический алгоритм двоичного деления. Далее вычисленные значения комбинируются по модулю m с использованием реконструкции рациональных чисел с тем, чтобы получить Т и Q.

Недостатком данного алгоритма является необходимость оценивать величину m для каждого вида ряда отдельно; в [3] для этого предлагается достаточно трудоемкий метод, на примере которого оценивается m для ряда £(3)-Недостатком алгоритма из [3] является также необходимость реализовывать достаточно сложный алгоритм реконструкции рациональных чисел (rational number reconstruction), имеющий временную сложность О(M(n) log(n)) и ёмкостную сложность О (п).

3. Основной алгоритм класса Sch(FQLINTIME//LINSPАСЕ). Модифицируем метод двоичного деления для гипергеометрических рядов так, чтобы алгоритм получился простой и при этом вычисления находились в пределах класса Sch(FQLINTIME//LIN SPACE).

Пусть требуется вычислить значения гипергеометрического ряда (1) с точностью 2~п. Для этого достаточно вычислить частичную сумму (2) с точностью 2_(га+1\ так как

\S — + 1))*| ^ \S — + 1))| + \S(/j,(n + 1)) — + 1))*| ^

< 2~(n+l) _|_ 2-(га+1) = 2~п

Здесь ЗДп + 1))* — обозначение для приближенного значения + 1)).

Обозначим г = ц(п + 1). Возьмем минимальную величину к\ такую, что 2kl ^ г; пусть г\ = \г/к\\. Запишем частичную сумму (2) в виде

P{r) = ai+T2 [а2 + т3[(73 + ... + Tkl-i[(Tkl-i +Tklakl]]] , (5)

где <

£з + о- з,

Величины а[ будем вычислять классическим методом двоичного деления для суммы (3), где г\ = (£ — 1)г\ +1, *2 = t•r\ — 1; для а\ возьмем г\ = 0, г2 = Г\ — 1.

Введём следующие обозначения: ш = 1{Ш) + 1, IV = тах(А,В), где А = = тах(|а(г)|, |Ь(г)|), В = тах(|р(^')|,\д^)\), 1(и) — длина битового представ-

г=0 ..г j=0..r

ления и. Оформим оценку вычислительной сложности расчёта а* и г* в виде двух лемм.

Лемма 1. Временная сложность алгоритма двоичного деления для вычисления <т* ограничена сверху 0(М(г) к^(г)); ёмкостная сложность ограничена сверху О (г).

Доказательство. Рассмотрим произвольную максимальную цепочку рекурсивных вызовов, полученных в результате вычислений в соответствии с формулами (4). Условимся, что нумерация в цепочке начинается с самого глубокого элемента цепочки: г = 1,... где — длина цепочки, £ ^ [~к^(2г1)~|.

Пользуясь методом математической индукции по г, покажем, что длина представления числа Т на глубине г при вычислении а[ удовлетворяет соотношению 1(Тг) < 2г+1ш + 2\ При этом заметим, что 1(Рг) ^ 2гш, ^ 2гш, 1(Вг) ^ 2гСс>, так как при увеличении г происходит удвоение длины числа. База индукции: г = 1, 1{Т\) ^ 2ш < 22ш + 21; индукционный переход:

здесь учитывается тот факт, что все коэффициенты а(г), Ь(г), р(.]), д{у) являются полиномами.

Оценим временную сложность вычисления а[. При этом будем учитывать свойство функции сложности умножения: 2М(2~1т) ^ М(т) (квазилинейные и полиномиальные функции удовлетворяют этому свойству). Так как в узле дерева вызовов на уровне г производится Сз умножений 2?-г чисел, длины которых не превосходят 2г+1ш + 2г, получаем следующую оценку количества операций, требуемых для вычисления а'г:

1(ТШ) < 2'1ш + 2*и + 2+ 2* + 1< 2(т)+1ш + 2Ш. Отметим, что исходя из этой оценки Тя имеет длину О (г), так как

ВД < С&и < с2-^— \оё(г) = С2г; \ogir)

Т1те((т0 < С3 ^ 2?_*М(2*+1Сс) + 2*) < С4 ^ 2^М(2*+2и)

г=1

г=1

я я

= С4 ^ ш) 2?-*2-?+(*+2)М(2?(х;) <

г=1 г=1

^ СеяМ(г) ^ С7 к^(г)М(г)

(здесь используется неравенство для 2?Сс> из оценки для 1(ТЯ)). Заключительное деление дает О(М(г)) операций.

Теперь оценим ёмкостную сложность вычисления а\. Учитывая, что в узле дерева вызовов вложенности г количество памяти, расходуемой на временные переменные, — это С8(2%+1ш + 2г), получаем, что количество памяти во всех одновременно существующих рекурсивных вызовах оценивается следующим образом:

8расе(а[) ^ ^ С8(2г+1ш + 2г) ^ Сд(2яш + 2?) ^ Сюг = О (г). □

г=1

Лемма 2. Временная сложность алгоритма двоичного деления для вычисления числителя и знаменателя т* ограничена сверху О(М(г) к^(г)); ёмкостная сложность ограничена сверху О (г).

Доказательство. Оценки вычислительной сложности т* совпадают с таковыми для <т* в силу того, что неравенства из доказательства леммы 1 также актуальны и для г*, если вычислять числитель и знаменатель <т* с помощью алгоритма двоичного деления для произведений. □

Рассчитывать приближённые значения Р(г)* с точностью 2_(га+1) по формуле (5) будем в соответствии со следующим итеративным процессом:

Ъ>\{т) = я*к1,

Ь*(т) = о*к1-г+\ +г/с1-г+Л-1, 1 = 1,...,кг, (6)

Ы(т) = Ы(т) +£*;

при г = к\ полагаем Р(г)* = Ьк1{т). Здесь т ^ г (величину т выберем

позже); а*, т* — приближения сг^, т% с точностью 2~т. Величины /гДт) полу-

чаются отбрасыванием битов <?т+1<?т+2 • • • Чт+з чисел Ы(т) после двоичной точки, начиная с (т + 1)-го, т. е.

М = IЫ(т) -Ы(т)I = 0.0...0qm+1qm+2...qm+j, (7)

а знак ei совпадает со знаком Ы(т) (ясно, что \вг\ < 2~т).

Предположим, что выполняются следующие условия:

\Ь(г)\ ^ 2 для всех г, \рШ/\яШ ^ 1 для всех 3. (8)

Покажем, что тогда верны следующие две леммы.

Лемма 3. Для любого г = 1,2,... ,к\ справедлива оценка

\Ы(т)\ < (г + 1)г{\¥. (9)

Доказательство. Применим математическую индукцию по j для hj(m). База индукции при j, равном 1: \h\(m)\ ^ r{W + 2~т < 2r{W. Индукционный переход для (j + 1) ^ 2:

I^J + l(m)l = Wk-i-ij+Vj + l + Tfci-(j + l)+2^i +£J + ll ^

^ ^riW + (j + l)r\W + 2~m < {{j + 1) + l)r\W. □

JIemma 4. Погрешность вычисления hkl(m) no схеме (6) оценивается как

A(k\, m) < 2~mrnk\W.

Доказательство. Обозначим Н\ = akl, Hi = akl-i+\ + r)(i,m) = \hi(m) — Hi\. Воспользуемся методом математической индукции для r](j,m) по j. База индукции при j, равном 1:

г](1,т) = \hi(m) - Hi\ = \a*kl - akl \ < 2~m.

Индукционный переход для (j + 1) ^ 2:

r](j + 1 ,т) = \a*ki_{j+)+1 + т*к1_и+1)+2^{т) + ej+l-

~ <Jfci-(j+l)+l — rfci-(j+l)+2-^il <

< |T*hj(m) — Tvhj(m) + Tvhj(m) — rvHj\ + 2 • 2~m ^

^ 2~mhj{m) + r](j, m) + 2 • 2~m.

Так как из (9) \hj(m)\ < (j + l)r\W и, по индукционному предположению, r](j,m) < 2~mmj2W, можно записать

r](j + 1, m) < 2~m(j + 1 )nW + 2~mmj2W + 2 • 2~m <

< 2~m(j + 1 )mW + 2~mm(j + 1 )2W = 2~mm(j + 2)2W.

Из A(fci,m) = r](ki,m) получаем искомое неравенство. □

Из леммы 4 следует, что достаточно взять m такое, чтобы выполнялось

т ^ (п + 1) + |~2 log(n + 1) + 2 log(r) + log(H^)], (10)

для вычисления Р(г)* с точностью 2_(-га+1).

Обозначим алгоритм расчёта гипергеометрического ряда, использующий схему (6), через LinSpaceBinSplit (linear space binary splitting).

Алгоритм LinSpaceBinSplit. Приближённое значение ряда (1).

Вход: Запись точности вычисления 2~п.

Выход: Приближённое значение ряда (1) с точностью 2~п.

Описание:

1) вычисляем г := ц(п + 1); подбираем к\ так, чтобы 2kl ^ г; вычисляем Г\ := \r/k\~\] рассчитываем W;

2) рассчитываем m по формуле (10);

3) h := a*ki (с помощью обычного алгоритма двоичного деления с точностью 2~т)]

4) выполняем цикл по г от 2 до к\:

а) рассчитываем V\ := aki_i+l с точностью 2~т с помощью

обычного алгоритма двоичного деления и v2 := i+2 с точ-

ностью 2~т,

б) вычисляем выражение h := v 1 + г^Л.,

в) h присваиваем величину h, округлённую в соответствии с (7);

5) на выход записываем h.

Оценим временную вычислительную сложность данного алгоритма, учитывая, что гит линейно зависят от п:

- O(log(n)) вычислений at дают 0(M(n) log(n)2);

- O(log(n)) вычислений Tt дают 0(M(n) log(n)2);

- O(log(n)) умножений чисел длины 0(п) дают 0(М(п) log(n));

итого получаем О(М{п) log(n)2) битовых операций. Емкостная сложность алгоритма LinSpaceBinSplit—О(п), так как во всех вычислениях в данном алгоритме фигурируют числа длины О (п).

Теорема. Модифицированный алгоритм двоичного деления для расчёта гипергеометрических рядов LinSpaceBinSplit принадлежит классу сложности Sch(FQLINTIME//LINS РАСЕ).

Заключение. В табл. 1, 2 приведены формулы и ряды для расчёта некоторых часто используемых на практике констант математического анализа. Данные ряды линейно сходятся (см. [2-4]), и они удовлетворяют условиям (8); следовательно, для расчёта их приближённых значений можно использовать алгоритм LinSpaceBinSplit, то есть перечисленные константы принадлежат множеству конструктивных чисел Sch{FQLINTIME//IINSPACE)cf- Алгоритм LinSpaceBinSplit можно также использовать для вычисления приближений многих других констант и приближений элементарных функций в рациональных точках.

Если для умножения использовать алгоритм Шёнхаге—Штрассена с временной сложностью O(nlog(n) log log(п)), то временная сложность алгоритма LinSpaceBinSplit будет O(nlog(n)3 log log(n)); при использовании для умножения простого рекурсивного метода с временной сложностью О (nlog(-3)) временная сложность алгоритма LinSpaceBinSplit будет O(nlog(-3) log(n)2).

Отметим, что ряд константы е сходится со скоростью 2-0(ralog(ra)), поэтому временная сложность расчёта е в соответствии с алгоритмом LinSpaceBinSplit будет 0(М(п) log(n)), а ёмкостная — О (n/log(n)).

Таблица 1

Формулы для расчёта констант

Константы Формулы Ряды

е е = ехр(1) ехр(1) = 2£~оМ

7Г 7г = 16а — 4/3 а — arctg (5) — 25 £*=о ( 1)г2(2i+\)52., /3 — arctg (239) — 2239 £i=o ( 1) 2(2i+l)2392*

С(з) V^oo (-l)*(205i2+250i+77)((i+l)!)5(i!)5 2^i=0 2((2i+2)!)5

Таблица 2

Ряды для расчёта констант

Константы a(i) b(£) pU) ч{з)

e 1 2 1 j

7Г 1 2(2* + 1) -1 52

1 2(2* + 1) -1 2392

C(3) 205*2 + 250* + 77 2 p(0) = 1, p(j) = -f 32(2.7 + l)5

Из дальнейших исследований можно отметить построение алгоритмов расчёта гипергеометрических рядов с рациональными коэффициентами с временной сложностью O(nlog(n)fc), к ^ 3, и при этом имеющих линейную ёмкостную вычислительную сложность.

БИБЛИОГРАФИЧЕСКИЙ СПИСОК

1. Schdnhage A., Grotefeld A.F. W, Vetter Е. Fast algorithms. A multitape Turing machine implementation. Mannheim, Germany: Bl-Wissenschaftsverlag, 1994. 298 pp.

2. Haible B., Papanikolaou, T. Fast multiprecision evaluation of series of rational numbers / In: Algorithmic number theory: Proceedings of 3rd international symposium (Portland, OR, 1998) / Lect. Notes Comput. Sci., 1423. Berlin: Springer, 1998. Pp. 338-350.

3. Cheng H., Gergel B., Kim E., Zima E. Space-efficient evaluation of hypergeometric series // SIGSAM Bull., 2005. Vol. 39, no. 2. Pp. 41-52.

4. Yakhontov S. V. FLINSPACE constructive real numbers and functions (in Russian). Saarbriicken, Germany: Lambert Academic Publishing, 2010. 167 pp.

5. Gourdon X., Sebah P. Binary splitting method: http://numbers.computation.free.fr/ Constants/Algorithms/splitting.ps.

6. Ко К.-I Complexity theory of real functions. Progress in Theoretical Computer Science. Boston, MA: Birkhauser Boston, Inc., 1991. 310 pp.

Поступила в редакцию 01/11/2011; в окончательном варианте — 24/VIII/2011.

MSC: 03F60; 68Q17

CALCULATION OF HYPERGEOMETRIC SERIES WITH QUASI-LINEAR TIME AND LINEAR SPACE COMPLEXITY

S. V. Yakhontov

St. Petersburg State University, Mathematics and Mechanics Faculty,

28, Universitetskiy prosp., Stariy Petergoff, St. Petersburg, 198504.

E-mails: sergey_home_mail@inbox. ru

A simple for practical implementation algorithm with the time complexity 0(М(те) log(n)2) and space complexity О (те) for the evaluation of h/ypergeom,etric series with rational coefficients on the Schdnhage machine is constructed (here M(n) is the complexity of integer multiplicatio'ri). It is shown that this algorithm is suitable in practical informatics for constructive analogues of often used constants of analysis.

Key words: constructive real numbers, hypergeometric series, quasi-linear time, linear space complexity.

Original article submitted 01/11/2011; revision submitted 24/VIII/2011.

Sergey V. Yakhontov (Ph.D. (Phys. & Math.)), Associate Professor, Dept, of Informatics.