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

Сер. 10. 2011. Вып. 4

ВЕСТНИК САНКТ-ПЕТЕРБУРГСКОГО УНИВЕРСИТЕТА

УДК 519.688+510.52 С. В. Яхонтов

ЭФФЕКТИВНОЕ ПО ВРЕМЕНИ И ПО ПАМЯТИ ВЫЧИСЛЕНИЕ ЭКСПОНЕНЦИАЛЬНОЙ ФУНКЦИИ КОМПЛЕКСНОГО АРГУМЕНТА НА МАШИНЕ ШЁНХАГЕ

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

Основные сведения о двоично-рациональных числах и конструктивных числах и функциях можно почерпнуть в [2]. Посредством Sch(FQLINTIME//LINSPACE) будем обозначать класс алгоритмов, квазилинейных по времени и линейных по памяти при вычислении на машине Шёнхаге. Под квазилинейностью понимается ограниченность сверху функцией вида O(n log(n)k) при некотором к.

Везде далее n будет обозначать длину записи точности вычисления 2-n двоичнорациональных приближений; x будем использовать для вещественного аргумента, z -для комплексного аргумента. Под функцией log(k) будем понимать логарифм по основанию 2. Через M(n) обозначим сложность умножения целых чисел.

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

В [3] показано, что многие вещественные элементарные функции можно вычислять алгоритмами класса FPTIME//LINSPACE, под которым понимается класс алгоритмов, вычислимых на машинах Тьюринга с полиномиальным временем и линейной памятью. Алгоритмы из [3] можно использовать для расчета элементарных функций на машине Шёнхаге также с полиномиальным временем и линейной памятью. В [4] предложены алгоритмы для быстрого вычисления экспоненциальной и некоторых других элементарных функций с временной сложностью O(nlog(n)3 log(log(n))); количество памяти, используемой данными алгоритмами, ограничено сверху, как будет показано ниже, O(n log(n)). И, далее, в [5] рассматриваются квазилинейные по времени и по памяти алгоритмы вычисления элементарных функций, основанные на разложениях в ряд Тейлора; количество памяти, используемой алгоритмами из [5], ограничено сверху квазилинейной функцией, так как здесь используется метод двоичного деления (binary splitting), требующий в классическом варианте O(nlog(n)) памяти для записи промежуточных результатов на машине Шёнхаге.

Яхонтов Сергей Викторович — кандидат физико-математических наук, доцент кафедры информатики математико-механического факультета Санкт-Петербургского государственного университета. Количество опубликованных работ: 9. Научные направления: конструктивный математический анализ, формальные спецификации, теория вычислительной сложности. E-mail: [email protected].

© С. В. Яхонтов, 2011

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

В статье дается положительный ответ на данный вопрос для экспоненциальной функции комплексного аргумента и некоторых других комплексных элементарных функций. При этом для построения Sch(FQLINTIME//LINSPACE) алгоритма вычисления на машине Шёнхаге экспоненциальной функции используется комбинация двух алгоритмов: алгоритма класса Sch(FQLINTIME//LINSPACE) расчета гипер-геометрических рядов и модификация метода быстрого вычисления экспоненциальной функции [4].

2. Описание вычислительной модели (машины Шёнхаге). Данная машина, введенная в [1], оперирует символами из алфавита Я = {0,12Л — 1} и с последовательностями таких символов (константу Л можно взять, например, равной 32). Машина состоит из одномерных массивов То, ..., Тт для чтения и записи символов из Я, регистров А, В, С, М для арифметических операций над символами (которые интерпретируются как натуральные числа) и управляющего устройства СРи. Массивы потенциально бесконечны в обе стороны; для каждого массива есть указатель р^ на текущий символ, записанный в массиве. Запись < р + ] > означает символ, на который ссылается указатель р + ]. Есть также дополнительный регистр У, который является указателем на текущую выполняемую инструкцию, и дополнительный массив S, потенциально бесконечный в одну сторону, который служит в качестве стека рекурсивных вызовов. Битовый регистр E служит как регистр переполнения при арифметических операциях.

Программа для машины Шёнхаге состоит из нескольких модулей-процедур на языке TPAL, который аналогичен языку ассемблера для МБС процессоров. В TPAL имеются команды загрузки в регистр символа, записанного в массиве, чтения из регистра и запись в массив, арифметические операции над символами, команды увеличения/уменьшения содержимого регистров, команды сдвига, вызов и возврат из процедуры, переход по метке и условный переход. Целые числа, которыми оперирует машина, кодируются как последовательности символов в алфавите Я:

а = ао + а,12х + а2(2х)2 + ... + а^-1(2х)к 1, 0 ^ а ^ 2х — 1.

Бит знака размещается, например, в символе, следующим за старшим символом а^~ 1.

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

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

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

приближения аргументов; по таким приближениям машина вычисляет приближения функции. Условимся, что запрос к оракулу в виде записи точности вычисления записывается в массиве То, в котором также дается приближение аргумента. При оценке временной вычислительной сложности на оракульной машине Шёнхаге запрос к оракулу учитывается как одна операция.

Определение 1. Емкостную вычислительную сложность алгоритма при расчете на оракульной машине Шёнхаге определим как сумму длин используемой памяти по всем массивам плюс максимум длины памяти, занятой стеком.

Отметим, что машина Шёнхаге существенно отличается от машины Тьюринга, РАМ и РАСП [6] возможностью использования рекурсивных процедур. Поэтому из верхних оценок вычислительной сложности (временной и емкостной) для машины Шёнхаге непосредственно не следуют какие-либо верхние оценки сложности для реализации тех же алгоритмов на машинах Тьюринга, РАМ или РАСП.

В то же время машина Шенхаге может рассматриваться как паскалевидная функция с теми же верхними оценками сложности. Из основного вывода в [7] и результатов этой статьи получаем, что есть алгоритм вычисления экспоненциальной функции комплексного аргумента в каждом круге, принадлежащий ЕP.

3. Конструктивные комплексные числа и функции. Под приближением с точностью 2-п комплексного числа г будем понимать комплексное число г' такое, что \г — г'\ ^ 2-п (здесь фигурирует модуль комплексного числа).

Пусть имеются комплексные числа ш = х + іу, ш' = х' + іу' такие, что \х — х'\ ^ 2-(п+1) и \у — у'\ ^ 2-(п+1). Тогда

- ш'\ = л/(х - X1)2 + (у- у')2 < \І2 • 2-2(п+і) < 2~п. (1)

То есть для вычисления приближения комплексного числа ш с точностью 2-п достаточно рассчитывать приближения действительной и мнимой частей этого комплексного числа с точностью 2-(п+1).

Будем говорить, что последовательность ф : N ^ Б х Б, ф(п) = (фх(п), фу(п)), Б -множество двоично-рациональных чисел, двоично-рационально сходится к комплексному числу г, если для любого п Є N выполняется ргес(фх(п)) = п + 2, ргес(фу(п)) = п + 2 и \г(п) — г\ ^ 2-п, где г(п) = фх(п) + іфу(п), ртес(т) - точность представления числа г, под которой понимается число битов справа от двоичной точки. Множество всех функций ф, двоично-рационально сходящихся к числу г, обозначим СЕг. Комплексное число г будем называть СЕ конструктивным, если СЕг содержит вычислимую функцию ф.

Определение 2. Число г Є Є называ,ется Sch(ЕQLINTIME//LINSPACE) конструктивным комплексным числом, если существует функция ф Є СЕг, вычислимая в пределах класса вычислительной сложности Sch(ЕQLINTIME//LINSPACE).

Пусть / - функция /(г) : А ^ С, где А = (г Є С : \г\ ^ К).

Определение 3. Функция /(г) называется Sch(ЕQLINTIME//LINSPACE) конструктивной комплексной функцией в круге А, если для любого г из этого круга существует функция ф из СЕ}(г), вычислимая в пределах класса вычислительной сложности Sch(ЕQLINTIME//LINSPACE).

Отметим, что для расчета комплексной конструктивной функции от аргумента г = х + іу требуется задать функции и(х, у) = Ие (/(г)) и у(х, у) = 1т (/(г)), для вычисления которых, в свою очередь, нужны две оракульные функции, соответствующие действительной и мнимой частям аргумента.

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

ехр(^) = ехр(ж + гу) = ехр(ж) • ехр(гу), (2)

а ряд ехр(гу) сводится к рядам с вещественными аргументами.

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

а(ъ) т-г рС?)

^ Ыг) И д(Л 7

г=0 У 1 3=0

где а, Ь, р, ^ - полиномы с целыми коэффициентами. Линейно сходящиеся гипергеомет-рические ряды используются для расчета многих констант математического анализа и элементарных функций в рациональных точках; ряд линейно сходится, если его частичная сумма

М(к) /• \ г ( '\

*м*» = Е$П^. (»

ъ=0 у 1 3 = 0

где (л(к) - линейная функция от к, отличается от точного значения не более чем на 2-к:

\Б - Б(и(к))\ < 2-к.

В классическом варианте метод двоичного деления состоит в следующем. Обозначим к\ = 1л(к). Рассмотрим частичную сумму (3) для некоторых чисел %\ и *2,

0 ^ *1 ^ к\, 0 ^ *2 ^ к\, *1 ^ %2’.

г а(г)р(%1) ...р(г)

ас -\ а\ЧР\г 1

П) =

Будем определять величины

Р(іі, І2) = _р(іі) ■ ■ .р(І2), Q(Іl, І2) = ч(Іі) ■ ■ ■ Я(І2), В(іі, І2) = Ь(іі) ■ ■ ■ Ь(І2),

Т(іі, І2) = В(іі, І2)Я(Іі, І2^(іі, І2)■

Если іі = І2, то они вычисляются напрямую. Иначе ряд делится на две части, левую и правую, и Р(іі,і2), д(іьі2), В(іі,і2) рассчитываются для каждой из частей рекурсивно. Затем полученные величины комбинируются:

Р (Іі ,І2 ) = Рі Рг , ^(Іі,І2)= ЯїЯг, В (Іі,І2)= ВіВг, Т (Іі,І2)= Вг Qr Ті + ВіРі Т ■ (4)

Алгоритм начинает свою работу с іі = 0, І2 = кі. После вычисления Т(0,кі), В(0,кі), Q(0, кі) осуществляется деление Т(0, кі) на В(0, кі^(0, кі), чтобы получить результат с заданной точностью. Выпишем алгоритм двоичного деления в явном виде.

Алгоритм ВтБрШ. Приближенное значение частичной суммы (3) с точностью 2-к.

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

Выход: Приближенное значение (3) с точностью 2-к.

Описание:

1) к1 := ц(к);

2) [Р,д,Б,Т] := ВтБрЩКвсит8(0,к1);

3) осуществляем деление г := -щ, с точностью 2~к;

4) возвращаем результат г.

В данном алгоритме используется подалгоритм, рассчитывающий рекурсивно Р,

д, б, т.

Алгоритм ВтЯрШКесигн. Вычисление величин Р, д, Б, Т.

Вход: Границы интервала *1, *2.

Выход: Кортеж [Р(гь *2), д(*1,*2), В(н, *2), Т(*1, *2)].

Описание:

1) если *1 = *2, то Т := а^^н) и возвращаем кортеж [р(*1), д(н), Ь(*1), Т];

2) гтЫ:= |>±*|;

3) вычисляем [Р;, Б;,Т] := БinSplitRecuгs(il,imгd);

4) вычисляем [Рг,дг,БГ,ТГ] := БinSplitRecuгs(imгd,І2);

5) Т := БгдгТ + Б;Р;ТГ и возвращаем кортеж [Р;Рг, д;дг, БВГ, Т].

Длины чисел Т(0, к1) и Б(0,к1)д(0,к1) пропорциональны к 1с^(к), т. е. алгоритм двоичного деления является квазилинейным по памяти; временная сложность данного алгоритма - 0(М(к)^(к)2) [5].

5. Метод быстрого вычисления ехр(х). Рассмотрим ряд Тейлора вещественной экспоненциальной функции

г 2 п

, . х '*> х X X X , ч

ехр(х) = ^7Г = 1 + 1Т + ^Г + -- - + ^Г + -- -

г=0

в точке хо, — \ + 2~т < жо < | — 2~т. Будем сначала вычислять значение этого ряда с точностью 2-Пз, пз > 7, в двоично-рациональной точке хт, \хт — хо\ ^ 2-т,

— т < % < т- Пусть /г - наименьшее такое, что

4

к

а т

п3 + 1 ^ 2к, а т = 2к+1. (6)

Представим хт = ±0.00«з«4 ... атат+1 в виде

хт = ±0.00аза4 + ±0.0000а5ава7а8 + ... +

+ ±°.°° . . . 0ат-2к + 1ат-2к + 2 . . . атат+1 =

_ /?2 Рз , /?4 Рк + 1 _

-7^ + ^+у1 + --- + - 72 + 73 + • • • + 7й+Ь

где (З2 = ±аза4, вз = ±а5ава7а8, ... вк+1 = ±ат_2к+1ат_2ь+2 ...атат+1; =

в^2-2 , 2 ^ V ^ к +1; - 2^-1-значное целое число. Значение ехр(хт) запишем как

произведение

ехр(хт) = ехр(72)ехр(7з) . ..ехр(7к+1). (7)

В методе быстрого вычисления экспоненты (БВЕ, или методе Карацубы; данный метод известен так же как Brent trick [5]) величины exp(7v) далее рассчитываются с помощью ряда Тейлора (5):

в в2 вг

ехр(7,) = 1 + ^+2!^ + ... + ^г + ад=6 + ад, (8)

здесь r = m2-v+1. Так как для остаточной суммы выполняется [4] неравенство

в |r+1

|Д„(г)| < 2-

Jv I

'(r + 1)!2(r+1)2v ’

то Rv(r)| < 2-m. Величины получаются из формул

av

:2

Си = т~, аи = £иЪи, Ьи = г! 2

Ьи

в которых целые аи вычисляются с помощью некоторого последовательного процесса группировки членов ряда (8). Отметим, что в формуле для Ьи присутствует факториал г, а г меняется от т2-1 до т2-к. Следовательно, здесь фигурирует величина, пропорциональная пз!, и поэтому длина промежуточных данных вычислений имеет порядок пз ^(пз).

6. Модификация метода двоичного деления. Модифицируем метод двоичного деления для гипергеометрических рядов так, чтобы вычисления находились в пределах класса Sch(FQLINTIME//ЫКБРЛОЕ), приспособив данный алгоритм для расчета серии рядов (8), определяющих величины ехр().

Возьмем г = т2_1/+2, к\ = 1с^(г), г\ = {к\ - натуральное число) и запишем

частичную сумму ряда (8) в виде

Р Ьи) = <71 + Т2[ст2 + тз[аз + ... + г^-1[а^-1 + тк1 ак1 ]]], (9)

где

т 1 I ^ I ^ I I ^

а1 = 1 + 777^ + г> |г>9. 9^ + • • • +

1!22" 2!22 2" (r1 - 1)!2(п-1)2"

д

Т2

в:1

r1!2r i2v 1

Pv Pi IZ1 1

a2 1 + (n + 1)22" + (n + l)(n + 2)22 2" + " ' + (n + 1). . . (2n - l)2(n-i)-2>' ’ тз

в:1

(r 1 + 1)... 2r12ri 2v ’

т 1 I ^ I ^ II ^

a3 = 1 + 777-, + 77;---, ,Wrt--, n„ +... + ■

(2r1 + 1)22v (2r1 + 1)(2r1 + 2)22 2v (2r1 + 1)... (3r1 - 1)2(ri-1)'2v

Величины а будем вычислять классическим методом двоичного деления для суммы

(3), где

(л (к) = Г1 — 1, а(г) = 1, Ь(г) = 1,

3 =0, ,л_ /1, 3 = 0, (10)

Р3 \ ри, 3 =0, 43 {((£ — 1)г1 + 3 )22", 3 =0.

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

Лемма 1. Временная сложность алгоритма двоичного деления для вычисления а на машине Шёнхаге ограничена сверху 0(М(m)log(m)); емкостная сложность -ограничена сверху О(т), где т находится по формуле (6).

Доказательство. Рассмотрим произвольную максимальную цепочку рекурсивных вызовов, полученных в результате вычислений в соответствии с алгоритмом ВтБрШЯесигв, а также пары (Р, г), ^, г), (В, г), (Т, г), в которых г - номер элемента цепочки рекурсивных вызовов, а величины Р, Q, В, Т находятся в г-м элементе цепочки. Условимся, что нумерация в цепочке начинается с самого глубокого элемента цепочки, г = 1,...,^, где ? - длина цепочки, <; ^ [^(2^)].

Пользуясь методом математической индукции по г, покажем, что длина представления числа Т в паре (Т, г) удовлетворяет соотношению 1(Т) < 2*/(и)+2г; здесь и = г22 . При этом заметим, что /(@„) < /(и), так как [3„ - 2^-1-значное целое число. Далее,

1(Р) < 2г/(и), /(Q) < 2г/(и) для пар (Р,г), (Q,г), так как при увеличении г происходит удвоение длины числа; В всегда равно 1.

База индукции: г =1, /(Т) ^ 2/(п) < 21/(и) + 21. Это неравенство следует из формул (10) и равенства Т = а(г1)_р(г!) для (Т, г). Индукционный переход:

/(Т) < 2*/(и) + 2*/(и) + 2* + 1 < 2*+1/(и) + 2*+1.

Это неравенство следует из формул (4): здесь будет два умножения чисел длины 2*/(и) и 2*/(и) + 2* и одно сложение.

Отметим, что, исходя из данной оценки, число Т для пары (Т, ?) имеет длину О(т), так как

Г

1{Т) < С&1{и) < + г) < С2

т2-^+12^

г +

log(r)

Оценим временную сложность вычисления аг. При этом будем учитывать свойство функции сложности умножения: 2М(2-1т) ^ М(т) (квазилинейные и полиномиальные функции удовлетворяют этому свойству).Так как в узле дерева вызовов на уровне г производится С4 умножений 2я-* чисел, длины которых не превосходят 2*/(и) + 2*, а <; ^ С log(m), то получаем следующую оценку количества операций, требуемых для вычисления аг:

Тгте(аг) < С4 ^ 2я-*М(2*/(и) + 2*) < С5 ^ 2я-*М(2*+1/(и)) *=1 *=1

Я Я

(и)) < Се *=1 *=1

^ С7^М(т) ^ С8 log(m)M(т)

С5^ 2я-*М(2я-(я-(т))/(и)) < С^ 2я-*2-я+(т)М(2я/(и)) <

(здесь используется неравенство для 2я/(и) из оценки для /(Т) для пары (Т, ?)). Заключительное деление дает 0(М(т)) операций. Если для умножения целых чисел использовать алгоритм Шёнхаге-Штрассена, то Тгте(аг) ^ т к^(т)2 loglog(m).

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

переменные, - это С8(2г/(и) + 2г), получаем, что количество памяти во всех одновременно существующих рекурсивных вызовах в фиксированной цепочке оценивается таким образом:

я

Брасе(аг) С9(2*/(и)+2*) < С10(2я/(и) + 2я) < С11 т = О(т).

*=1

Лемма доказана.

Лемма 2. Временная сложность алгоритма двоичного деления для вычисления тг на машине Шёнхаге ограничена сверху 0(М(m)log(m)); емкостная сложность -ограничена сверху О(т), где т определяется по формуле (6).

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

Лемма доказана.

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

к1(т1) = а1г,

к*(т1) = а*к1 -*+1 + Тк1-*+2^-Ь г = 1,...,kl, (11)

к(т\) = Ы(т1) + е*;

при г = к1 полагаем Р(уи)* = кк1 (т1). Здесь т1 ^ т, а* - приближения а* с точ-

ностью 2-т1, т* - приближения т* с точностью 2-т1. Величины к* (т{) получаются отбрасыванием битов дт1 + 1?т1+2 ...Цт1+з чисел к*(т1) после двоичной точки, начиная с (т1 + 1)-го, т. е.

Ы = |к*(т1) - к*(т1)| = 0.0 . .. Одт^дт^ . .. Чтщ+з, (12)

а знак е* совпадает со знаком Н*(т1) (ясно, что 1е*1 < 2-т1).

Отметим, что точность расчета ехр(7^) будет 2-т, так как

| ехр(7„)* - ехр(7„)| < 1С - &I + Щ(г)| < 2-(т+1) + 2-(т+1) = 2-т.

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

к*(т1) < 2. (13)

Доказательство. Применим математическую индукцию по 3 для к^ (т1), используя оценки

4 1

\щ\ < ехр(7„) < -, п < 7/1 <

База индукции при 3, равном 1: |к1(т1)| ^ ак1 + 2-т1 < 2. Индукционный переход для

(3 + 1) > 2:

|к5 + 1(т1)| = |ак1-(5'+1) + 1 + тк1 -а+1)+2к0 + еj + 1| ^

4

^3 +

}_ I o-mi

4

2 + 2—mi < 2.

Лемма доказана.

Лемма 4. Погрешность вычисления hki (mi) по схеме (11) оценивается как

Д(&1, mi) < 2—mi+ki.

Доказательство. Обозначим

Hi = aki, Hi = CTfci-j+i + rfci-i+2Hi-i, n(i,mi) = \hi(mi) - Щ.

Воспользуемся методом математической индукции для n(j,mi) по j. База индукции при j, равном 1:

n(l,mi) = \hi(mi) - Hi\ = \aki - aki \ < 2-mi + i.

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

V(j + 1,mi) = Ki-(j+)+i + Tki —(j+i)+2 hj (mi)+£j+i-

- aki — (j+i) + i - Tki — (j+i)+2Hj \ <

< \rUjhj (mi) - Tu hj (mi) + tv hj (mi) - tv Hj \ + 2 • 2—mi <

< 2—mi hj (mi ) + 2—2n( j, mi ) + 2 • 2—mi.

Так как из (13) \hj(mi)\ < 2 и, по индукционному предположению, n(j,mi) < 2—mi+j, то

n(j + 1, mi) < 2 • 2—mi + 2—22—mi+j + 2 • 2—mi < 2—mi+(j+i).

Из Д(^, mi) = n(ki, mi) получаем искомое неравенство.

Лемма доказана.

Из леммы 4 следует, что достаточно взять mi = 2m + 1, чтобы вычислять P(yv) с точностью 2— (m+i). Обозначим алгоритм расчета гипергеометрического ряда, использующий схему (11), через ModBinSplit (modified binary splitting).

Алгоритм ModBinSplit. Приближенное значение гипергеометрического ряда.

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

Выход: Приближенное значение ряда (8) с точностью 2—m.

Описание:

1) mi := 2m + 1;

2) h := a*k (с помощью обычного алгоритма двоичного деления с точностью 2—mi);

3) выполняем цикл по i от 2 до ki:

а) рассчитываем vi := aki+i с точностью 2—mi с помощью обычного алгоритма двоичного деления и V2 := tJ_^_i+2 с точностью 2—mi,

б) вычисляем выражение h := vi + V2h,

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

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

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

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

• O(log(n3)) вычислений ^ дают O(M(пз) log^)2);

• O(log(nз)) умножений чисел длины O(^) дают O(M(пз) log(nз));

итого получаем O(M(пз) log(nз)2). Емкостная сложность модифицированного алгоритма двоичного деления ModBinSplit - это O(^), так как во всех вычислениях в данном алгоритме фигурируют числа длины O(^).

Утверждение 1. Модифицированный алгоритм двоичного деления ModBinSplit принадлежит классу сложности Sch(FQLINTIME//LINSPACE).

7. Модификация метода БВЕ. Нам понадобится некоторая модификация в методе вычисления экспоненты БВЕ с тем, чтобы алгоритм находился в классе вычислительной сложности Sch(FQLINTIME//LINSP ACE).

В классическом алгоритме БВЕ значения exp(xm)J находятся по формуле (Т) с помощью попарного суммирования величин exp^)*; при этом одновременно нужно держать в памяти O(log(log(nз))) чисел длины O(^), и, как уже отмечалось, в данном алгоритме фигурирует величина пз!, т. е. емкостная сложность такого алгоритма хуже линейной. Поэтому рассчитывать exp(xm)J с точностью 2-Пз будем в соответствии со следующим итеративным процессом:

h2(m) = exp(Y2)J,

hi(m) = hi-l(m)exp(Yi)J, i = 2,...,k + 1, (14)

hi(m) = hi(m) + Єі;

при i = k + 1 полагаем exp(xm)* = hk+l(m). Здесь exp^)* - это приближения величин exp(Yi) с точностью 2-m, определенные по формуле (8). Величины hi(m) получаются отбрасыванием битов qm+lqm+2 ...qm+t чисел hi(m) после двоичной точки, начиная с (m + 1)-го, т. е.

= |hi(m) - hi(m)| = G.G . ..Gqm+lqm+2 ...qm+t, (15)

а знак єі совпадает со знаком hi(m) (ясно, что ^j\ < 2-m).

Лемма 5. Для любого i Є 2..k +1 справедлива оценка

h (m) < 2і-1. (16)

Доказательство. Применим математическую индукцию по j для hj(m). База индукции при j, равном 2:

з

\h2(m)\ < I exp(Y2)| + 2~m < - + 2~16 < 22-1 (здесь учитывается то, что І7*| < ^, то 16). Индукционный переход для (j + 1) 3:

|hj+l(m)| = h(m) exp(Yj+l)J + Є'+і| < hj(m)(exp(Yj+l) + 2-m) + 2-m <

< 2j

-l

з

_ + 2~m 2

< 2j.

Лемма доказана.

Лемма 6. Погрешность вычисления кк+1(т) по схеме (14) оценивается как

m

Д(к + 1, m) < 2—m+2(k+i).

Доказательство. Обозначим

H2 =exp(Y2), Hi = exp(Y2)exp(Y3) ...exp(Yi), n(i,m) = \hi(m) - Hi\.

Воспользуемся методом математической индукции для n(j, m) по j. База индукции при j, равном 2:

n(2,m) = \h2(m) - H2\ = \ exp(Y2)J - exp(Y2)\ < 2—m < 2—m+2(i+i). Индукционный переход для (j + 1) ^ 3:

V(j + 1 m) = \hj(m) exp(7j+i)j + £j+i - Hj exp(Yj+i)\ <

< \hj(m)exp(Yj+i)J - hj(m)exp(Yj+i) +

+ hj(m)exp(Yj+i) - Hj exp(Yj+i)\ + 2—m <

< hj(m)2—m + n(j, m) exp(Yj+i) + 2—m.

Так как из (16) \hj(m)\ < 2j—i и, по индукционному предположению, n(j,m) < 2—m+2j, то

3

r){j + 1 ,m)< + [2-m+2j] - < 2-m+2W+1).

Из Д(к + 1, m) = п(к + 1, m) получаем искомое неравенство.

Лемма доказана.

Из леммы 6 следует, что точность 2—m вычисления exp(Yi)* достаточна для расчета exp(xm)j с точностью 2—Пз, так как

- m + 2(к + 1) ^ -пз ^ 2k+i - 2(к + 1) ^ пз,

2k+i - 2(к +1) > 2k+i2—i = 2k > n3 + 1

(напомним, что пз + 1 ^ 2k).

Обозначим алгоритм расчета вещественной экспоненциальной функции, использующий схему (14), через ModFEE (modified fast exponenta evaluation).

Алгоритм ModFEE. Приближенное значение вещественной экспоненциальной функции.

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

Оракульные функции: фх.

Выход: Приближенное значение exp(x) с точностью 2—Пз.

Описание:

1) h := exp(Y2)* (с помощью алгоритма ModBinSplit с точностью 2—m);

2) выполняем цикл по i от 3 до к + 1:

а) рассчитываем vi := exp(Yi)J с помощью алгоритма ModBinSplit с точностью

2—m,

б) вычисляем выражение h := h • vi,

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

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

Вычислительная сложность данного алгоритма при расчете на машине Шёнхаге будет следующей: временная сложность - 0(М(пз)1<о;(пз)3), так как алгоритм расчета гипергеометрических рядов МойЕтБрШ применяет 0(М(пз) 1og(nз)2) операций, а в схеме (14) используется 0(к^(пз)) таких вычислений и 0(1og(nз)) умножений чисел длины 0(пз); емкостная сложность - 0(пз), так как во всех вычислениях в данном алгоритме фигурируют числа длины 0(пз).

Утверждение 2. Модифицированный алгоритм быстрого вычисления экспоненты МоёЕЕЕ принадлежит классу сложности Sch(FQLINTIME//LINSP АСЕ).

8. Вычисление комплексной функции ехр(г). Пусть везде далее величина р - натуральное число, р ^ 0. Функцию ехр(г) будем рассматривать с точностью 2-п в круге \г\ ^ 2Р; тогда будет выполняться \х\ ^ 2Р, \у\ ^ 2Р.

Проведем мультипликативную редукцию интервала для комплексного аргумента. А именно, возьмем целое число в = 2Р+3 и г' = тогда \г'\ = т. е. г' будет в круге \г'\ ^ 2~з и соответственно будет выполняться \х\ ^ 2~з, \у\ ^ 2~з. Тем самым вычисление ехр(г) сводится к расчету ехр(г'), после чего результат возводится в степень в. Легко увидеть, что для показателя пі точности вычисления ехр(г') функция зависимости будет пі = L(n) + С(р), где L(n) - линейная функция от п, С(р) - константа,

зависящая от р (константа в том смысле, что она не зависит от п).

Будем рассматривать функцию ехр(г) в круге \г\ ^ 2~з. Обозначим (х = ехр(х), (у = ехр(іу) и перепишем формулу (2) в виде

ехр(г) = Сх(йе(Су) + *1ш(Су)) = Схйе(Су) + *Сх1ш(Су) = Сх + іСу ■

Согласно (1), для того чтобы получить значение ехр(г) с точностью 2-П1, величины £х и Су достаточно вычислять с точностью 2-(п1+1). Далее имеем

(у = Е

=0

(гуУ

з'-

= 1 , {У , (ІУ) 1! 2!

2

+

(*у)э

3!

+ ■■■ +

(іу)к к\

+ ■■■ =

2 з 4

, . -у у у . у .

1 + гЇЇ“¥“гЗ! + 4!+-"

24 у2 у4

1 . 2! 4!

+ і

з

У У . 1! “ 3! +

Отсюда видно, что |Ке(£у)| < 2, |1ш(£у)| < 2, и, следовательно, для расчета £х и £у с точностью 2-(”1+2) и последующим округлением до точности 2-(”1+1) нужно вычислять значения Сх, К-е(Су), !ш(^у) с точностью 2-”2, П2 = П1 + 5 (см. [3]). Так как все данные величины - ряды с вещественным аргументом, то для их расчета воспользуемся алгоритмом ЫвйЕЕЕ (некоторые члены ряда (8) для Ие(^у) и 1ш(£у) будут нулевыми). Пусть т ^ п2 + 3. Имеем следующее: хт = хо + в12-т, |01| ^ 1, при этом |хо| ^ 2-3

и

| ехр(хо) - ехр(хт)| = | ехр(хо) - ехр(хо)ехр(^2-т) | =ехр(хо)|1 - ехр(^2-т)|.

Так как ехр(012_т) < 1_^_Т71, то | ехр(#12~т) —1| < 122_Т71 < 2~т+1, и поэтому получаем оценку

| ехр(хо) - ехр(хт)| < 2 • 2-т+1 = 2-т+2 < 2-("2+1),

из которой видно, что точность вычисления хт должна быть не ниже 2-(”2+3) для достижения точности 2-(”2+1) значения функции ехр(хо). Поскольку т = 2к+1, то такое условие выполняется.

Теперь нужно учесть, что мы вычисляем приближенное значение ехр(хт)*. Если его определять с точностью 2—”3, пз = П2 + 1, то

| ехр(хо) - ехр(хт)* | < | ехр(хо) - ехр(хт)| + | ехр(хт) - ехр(хт)* | <

< 2 — (”2 + 1) + 2 —(”2 + 1) = 2 —”2

Отсюда следует, что в качестве пз для алгоритма МоСЕЕЕ нужно брать величину П2 + 1. Теперь все готово для описания основного алгоритма.

Алгоритм ЕхрУаЫв. Приближенное значение комплексной экспоненциальной функции.

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

Оракульные функции: фх и фу для аргумента г = х + гу.

Параметры: Константа р.

Выход: Приближенное значение ехр(г) с точностью 2—”.

Описание:

1) П1 := Ь(п) + С(р);

2) П2 := П1 + 5;

3) пз := П2 + 1;

4) рассчитываем к, т так, что выполняется (6);

5) Р1 := р + 3;

6) в := 2Р1;

7) вычисляем х* := фх(шах(1,т- Р1)), у* := фу(шах(1,т -Р1));

8) производим редукцию интервала (ж*)' = —, (у*)' = — (точность аргументов станет 2—т);

9) с помощью алгоритма МоСЕЕЕ рассчитываем «1 := СХ, «2 := К-е(£у)*, «з := 1ш(£у)* с точностью 2—”3, где аргументами служат (х*)', (у*)';

10) вычисляем двоично-рациональные числа «4 := «1«2 («4 = И,е(ехр(г))*),

:= «1«3 («5 = 1ш(ехр(г))*) с точностью 2—(”1+2) и округляем до точности 2—(”1+1); длина битового представления этих чисел будет п1 + 2;

11) на выход пишем комплексное число («4 + г • «5)я.

Свойства алгоритмов МоСБтБрШ и МоСЕЕЕ позволяют сделать выводы в виде следующих утверждений:

Утверждение 3. Алгоритм ЕхрУа1пе расчета комплексной экспоненциальной функции принадлежит классу сложности Sch(FQLINTIME//LINSP АСЕ).

Оценки вычислительной сложности алгоритма ЕхрУаЫв при расчете на машине Шёнхаге будут такие же, как у алгоритма МоСЕЕЕ, т. е. временная сложность -0(М(пз)ко?(пз)3), емкостная - 0(пз). Если для умножения целых чисел использовать алгоритм Шёнхаге-Штрассена, то временная сложность алгоритма ЕхрУаЫв будет ограничена сверху 0(п3 ^(п3)4 ^^(п3)).

Утверждение 4. Функция ехр(г) является Sch(FQLINTIME//LINSPACE) конструктивной в любом круге |z| ^ 2Р.

9. Вычисление комплексных функций вт(г), еов(г), вЬ(г), еЬ(г). Исходя из формул для тригонометрических функций

егг - е—гг егг - е—гг егг + е—гг

зт(г) = ---------- = г-------, соБ(г) = ---------,

^ у 2г -2 ’ ^ У 2 ’

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

Утверждение 5. Функция 8ш(г) является Sch(FQLINTIME//LINSPACE) конструктивной в любом круге |z| ^ 2Р.

Утверждение 6. Функция еов(г) является Sch(FQLINTIME//LINSPACE) конструктивной в любом круге |z| ^ 2Р.

Следующие два утверждения также вытекают непосредственно из формул для гиперболического синуса и косинуса:

pz _ p — z pZ I p — z

sh(z) =----------------, ch(x) =

2

Утверждение 7. Функция sh(z) является Sch(FQLINTIME//LINSPACE) конструктивной в любом круге \z\ ^ 2p.

Утверждение 8. Функция ch(z) является Sch(FQLINTIME//LINSPACE) конструктивной в любом круге \z\ ^ 2p.

10. Заключение. Алгоритм ExpValue расчета экспоненциальной функции можно применять в информатике (вместе с алгоритмами ModBinSplit и ModFEE) как основу Sch(FQLINTIME//LINSPACE) конструктивных комплексных функций exp(z), sin(z), cos(z), sh(z), ch(z), заданных на множестве Sch(FQLINTIME//LINSPACE) конструктивных комплексных чисел.

Отметим также, что если для умножения использовать простой рекурсивный метод с временной сложностью O(nlog(3)), то временная сложность алгоритма ExpValue будет O(nlog(3) log(n)3).

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

Литература

1. Schonhage A., Grotefeld A.F. W., Vetter E. Fast algorithms. A multitape Turing machine implementation. Leipzig: Brockhaus AG, 1994. 298 p.

2. Ko K. Complexity theory of real functions. Boston: Birkhauser, 1991. 310 p.

3. Яхонтов С. В. FLINSPACE конструктивные вещественные числа и функции. Saarbrucken: LAMBERT Academic Publ., 2010. 167 с.

4. Карацуба Е. А. Быстрые вычисления трансцендентных функций // Проблемы передачи информации. 1991. Т. 27, вып. 4. С. 76—99.

5. Haible B., Papanikolaou T. Fast multiple-presicion evaluation of series of rational numbers // Proc. of the Third Intern. Symposium on Algorithmic Number Theory. June 21-25, 1998. P. 338-350.

6. Ахо А., Хопкрофт Дж., Ульман Дж. Построение и анализ вычислительных алгоритмов / пер. с англ. А. О. Слисенко; под ред. Ю. В. Матиясевича. М.: Мир, 1979. 536 c. (Aho A., Hopcroft J., Ull-man J. The design and analysis of computer algorithms.)

7. Косовская Т. М., Косовский Н. К. Принадлежность классу FP дважды полиномиальных паскалевидных функций над подпрограммами из FP // Компьютерные инструменты в образовании. 2010. № 3. С. 3-7.

Статья рекомендована к печати проф. Л. А. Петросяном.

Статья принята к печати 19 мая 2011 г.