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

2013 Вычислительные методы в дискретной математике №2(20)

УДК 510.25+510.52+519.688

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

НА МАШИНЕ ШЁНХАГЕ

С. В. Яхонтов

Санкт-Петербургский государственный университет, г. Санкт-Петербург, Россия E-mail: SergeyV.Yakhontov@gmail.com

Строится алгоритм FLE быстрого вычисления логарифмической функции ln(1 + x) вещественного аргумента на полуинтервале [2-5,1 — 2-5) на машине Шёнхаге с оракульной функцией и даётся верхняя оценка его временной и емкостной сложности. Алгоритм FLE строится на основе разложения в ряд Тейлора по аналогии с алгоритмом быстрого вычисления экспоненты FEE, при этом дополнительно строится модифицированный алгоритм двоичного деления ModifBinSplit для гипергеометрических рядов. Для алгоритмов ModifBinSplit и FLE показывается квазилинейность по времени и линейность по памяти при вычислении на машине Шёнхаге, то есть принадлежность классу Sch(FQLIN-TIME//LIN-SPACE).

Для расчёта логарифмической функции на произвольном промежутке используется мультипликативная редукция интервала.

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

Введение

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

Приводится мера временной и емкостной сложности [1-4] вычислений на машине Шёнхаге и даётся верхняя оценка временной и емкостной сложности предлагаемого алгоритма вычисления логарифмической функции ln(1 + x) вещественного аргумента на полуинтервале [2-5,1 — 2-5) на машине Шёнхаге.

Основные сведения о двоично-рациональных числах и алгоритмических числах и функциях можно найти в [5]. Посредством Sch(FQLIN-TIME//LINS-PACE) будем обозначать, как и в [1-3], класс алгоритмов, квазилинейных по времени и линейных по памяти при вычислении на машине Шёнхаге. Под квазилинейностью понимается ограниченность сверху функцией вида O(nlog(n)k) при некотором k.

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

Известно, что многие вещественные элементарные функции, с одной стороны, вычислимы с помощью разложения в ряд полиномиальными по времени и линейными по памяти алгоритмами [3] и, с другой — квазилинейными по времени и квазилинейными по памяти алгоритмами [6-8]. В данных работах рассматривается сложность

алгоритмов при реализации на машине Тьюринга [3] и в рамках модели битовых вычислений [6-8], но оценки вычислительной сложности, приведённые в них, верны и для машины Шёнхаге.

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

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

Метод двоичного деления (англ. binary splitting) [9], предназначенный для вычисления гипергеометрических рядов, является рекурсивным вариантом метода FEE [6-8,10- 13] вычисления гипергеометрических рядов: в методе из [9] дерево вычислений обходится сверху вниз с помощью рекурсии, в отличие от итеративного метода FEE, в котором, в частности, дерево вычислений обходится снизу вверх. Будем использовать рекурсивный метод из [9]; это связано с удобством его реализации на машине Шёнхаге, в которой есть рекурсивные вызовы. К тому же при использовании рекурсии как сам алгоритм, так и оценки его вычислительной сложности получаются проще.

В данной работе подход, использованный в [2, 3] для экспоненциальной функции комплексного аргумента, применяется для логарифмической функции вещественного аргумента: алгоритм класса Sch(FQLIN-TIME//LIN-SPACE) вычисления логарифмической функции на машине Шёнхаге основан на комбинации алгоритма расчёта ги-пергеометрических рядов из [1 - 3] и алгоритма быстрого вычисления логарифма, который строится для логарифмической функции по схеме метода быстрого вычисления экспоненты.

Везде далее под функцией log(k) будем понимать логарифм по основанию 2; через Msch(n) обозначим верхнюю оценку временной сложности алгоритма Шёнхаге — Штрассена [4] умножения целых чисел на машине Шёнхаге.

1. Описание вычислительной модели (машины Шёнхаге)

Данная машина, введенная в [4], оперирует символами из алфавита Е = {0,1,... , 2Л — 1} и с последовательностями таких символов, где Л — некторая константа (данную константу можно взять, например, равной 32). Машина состоит из одномерных массивов Т0, ..., Тт для чтения и записи символов из Е, регистров A, B, C, M для арифметических операций над символами (которые интерпретируются в данном случае как натуральные числа) и управляющего устройства CPU. Массивы потенциально бесконечны в обе стороны; для каждого массива есть указатель p на текущий символ, записанный в массиве. Запись (p + j) означает символ, на который ссылается указатель p + j. Есть также дополнительный регистр Y — указатель на текущую выполняемую инструкцию, и дополнительный массив S, потенциально бесконечный в одну сторону, который служит в качестве стека рекурсивных вызовов. Битовый регистр E служит как регистр переполнения при арифметических операциях.

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

а = а0 + а12Л + а2(2Л)2 + ... + ак-1(2Л)к-1, 0 ^ а, ^ 2Л - 1.

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

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

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

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

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

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

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

2. Алгоритмические вещественные числа и функции

В данной работе основу представления конструктивных объектов составляет понятие алгоритмической последовательности <^, сходящейся по Коши [5], при этом в качестве вычислительной модели берётся машина Шёнхаге. Такая последовательность определяется на множестве всех натуральных чисел N включая 0, а областью аппроксимирующих значений является всюду плотное в К естественное подмножество

множества рациональных чисел. Для последовательности, сходящейся по Коши и задающей вещественное число х, требуют, чтобы выполнялось |<^(п) — х| ^ 2-п для любого натурального п.

В качестве множества аппроксимирующих значений берётся множество двоичнорациональных чисел О [5]. Рациональное число й называется двоично-рациональным, если й = т/2п для некоторого целого т и натурального п. Двоично-рациональные числа имеют конечное двоичное представление: строка в, равная ±ирир-1.. .щ.у1у2 ... vr, обозначает число

й = ± ( Е и*2* + т Е V 2-

\г=0 ]=1

Длина представления двоично-рационального числа определяется как количество символов в строке в, равное, с учётом знака и двоичной точки, р+г + 3, и обозначается / (в). Под точностью представления ргес(в) понимается число битов справа от двоичной точки, то есть г. С точки зрения изучения вычислительной сложности двоично-рациональные числа удобны тем, что для любого п двоично-рациональные числа с точностью п равномерно распределены на вещественной прямой [5].

Последовательность ^ ^ О двоично-рационально сходится к вещественному

числу х, если для любого п € N выполняется ргес(^(п)) = п +1 и |<^(п) — х| ^ 2-п. Множество всех функций, двоично-рационально сходящихся к вещественному числу х, обозначается СЕХ. Вещественное число х называется СЕ-алгоритмическим [5], если СЕХ содержит вычислимую функцию <^.

Вещественная функция f, заданная на отрезке [а,Ь] (или на любом другом промежутке), называется алгоритмической функцией [5] на этом отрезке, если существует машина Шёнхаге М с оракульной функцией, такая, что для любого х € [а, Ь] и любой вычислимой функции ^ € СЕХ функция ф, вычисляемая М с оракульной функцией <^, принадлежит CЕf (Х). Фактически, это означает, что для любой вычислимой функции ^ € СЕХ и любого п € N машина М последовательно вычисляет т € N и й € О, такие, что |<^(т) — х| ^ 2-т, |й — f (х)| ^ 2-п.

Сложность расчёта двоично-рациональных приближений алгоритмических чисел и функций определяется в [5] на основе длины двоичного представления точности вычисления. Память ленты запроса и ленты ответа оракульной функции при оценке емкостной вычислительной сложности алгоритма не учитывается. Обращение к ора-кульной функции ^ € СЕХ аргумента х алгоритмической функции осуществляется следующим образом:

— на ленту запроса оракульной функции записывается точность вычисления аргумента 2-т в виде 0т (унарная запись);

— рассчитывается значение ^(т) оракульной функции и результат записывается на ленту ответа;

— значение ^(т) считывается с ленты ответа в промежуточную память.

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

Определение 2 [1-3]. Число х € К назовём Sch(FQLIN-TIME//LIN-SPACE)-алгоритмическим вещественным числом, если существует функция ^ € СЕХ, вычислимая в пределах Sch(FQLIN-TIME//LIN-SPACE).

Определение 3 [2, 3]. Вещественную функцию f, заданную на отрезке [а,Ь], назовём Sch(FQLIN-TIME//LIN-SPACE)-алгорumмuческой вещественной функцией на отрезке [а, Ь], если для любого х € [а,Ь] существует Sch(FQLIN-TIME//LIN-SPACE)-вычислимая функция ф из CЕf(Х).

Множества Sch(FQLIN-TIME//LIN-SPACE) алгоритмических вещественных чисел и функций будем обозначать Sch(FQLIN-TIME//LIN-SPACE)CF и Sch(FQ-LIN-TIME//LIN-SPACE)C [а,6] соответственно. Здесь индекс С [а, Ь] обусловлен тем, что алгоритмические функции являются непрерывными на всей области определения [5]. Построение алгоритмического аналога вещественной функции f на отрезке [а, Ь] означает описание алгоритма, вычисляющего двоично-рациональные приближения с произвольной точностью значений f (х) для х € [а, Ь].

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

Под приближением £ с точностью 2-п будем понимать двоично-рациональное число ¿*, такое, что |£* — £| ^ 2-п. Для модуля выполняются неравенства

|£*| — |£| ^ |£* — £| ^ 2-п, то есть |£*| ^ |£| + 2-п,

|£| — |£*| ^ |£ — ¿*| ^ 2-п, то есть |£*| ^ |£| — 2-п.

Далее, пусть двоично-рациональное число £** таково, что |£** — £| ^ 2-(п+1). Отбросим все биты £** после двоичной точки, начиная с (п+2)-го, а новую величину обозначим ¿*. Тогда рге^*) = п +1 и |е| = |£* — ¿**| < 2-(п+1). Отсюда получаем

|£* — £| ^ |£* — ¿**| + |Г* — £| < 2-(п+1) + 2-(п+1) = 2-п.

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

3. Метод двоичного деления

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

^=е Ш п р(4. (1)

г=0 Ь(*) ]=0 )

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

5 ("№)) = 1 щ п § ■ (2)

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

|5 — 5(^(к))| ^ 2-к.

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

0 ^ ¿і ^ &і, 0 ^ ¿2 ^ кі, ¿і ^ ¿2:

^ а(г)р(гі) ...р(г)

Ь(.¿і, ¿2) = Е ил (■ \------Тл .

і=іі о(г)д(гі)... д(г)

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

Р(¿і, ¿2) = р(іі) . . .р(*2), ^(¿і, ¿2) = ^(¿і) . . . <?(*2), В(іі, ¿2) = Ь(іі) . . . Ь(І2),

Т(¿і, ¿2) = в(гі,г2)ф(гі,г2)£(¿і, ¿2).

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

Р (¿і ,¿2) = Рг Рг, ^(¿і, ¿2) = ОгО-, В (¿і, ¿2 ) = ВіВг, Т (¿і, ¿2) = В- О- Ті + Ві РіТ.

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

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

Вход: Запись 0к точности вычисления 2-к Выход: Кортеж [Р (¿і, ¿2), ^(¿і, ¿2), В (¿і, ¿2), Т (¿і, ¿2)]

1: кі := ^(к);

2: [Р, О,В,Т] := B¿nSp/¿íRecмrs(0, кі);

Т _к

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

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

— к

В алгоритме 1 используется подалгоритм БгпБрШКеспгз (алгоритм 2), рассчитывающий рекурсивно Р, О, В, Т.

Алгоритм 2. БтБрШКесптв. Вычисление величин Р, О, В, Т Вход: Границы отрезка ¿і, ¿2

Выход: Кортеж [Р (¿і, ¿2), ^(¿і, ¿2), В (¿і, ¿2), Т (¿і, ¿2)]

1: Если ¿і = ¿2, то

2: Т := а^)р^) и возвращаем кортеж [р^), ^(¿і), ^(¿0, Т];

¿і + ¿2 _ 2 _

4: вычисляем [Рг,Ог, Вг,Тг] := B¿nS,p/¿tЯecurs(¿1, ¿тіа); вычисляем [Рг, ОГ, Вг, Тг] := B¿nS'p/¿íЯecиrs(¿mid, ¿2);

Т := ВгОгТг + ВгРгТг и возвращаем кортеж [РгРг, ОО, ВгВг, Т].

3: ¿тіїі : =

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

4. Метод быстрого вычисления exp(x)

Рассмотрим ряд Тейлора вещественной экспоненциальной функции [16]

СО 'X <x2

exP(x) = О -тг = 1 + тт + -27 + ... + -7 + r„(x). (3)

i=0 1- 2 —•

Будем вычислять значение этого ряда с точностью 2-n в двоично-рациональной точке xm, |xm — x0| ^ 2-m, -1/4 < < 1/4. Пусть k — наименьшее число, такое, что

- +1 ^ 2k, и m = 2fc+1. Представим = ±0,00а3а4 ... amam+1 в виде

= ±0,00а3а4 + ±0,0000а5«ба7а8 + ... +

+ ±0,00 . . . 0am-2k+1am-2k +2 . . . amam+1 =

= в2 + в3 + в4 + + ^fc+1 = + + +

= 24 + 28 + 2!6 + ... + ~2^ = 72 + 73 + ... + 7fc+b

где в2 ±a3a4, e3 ±а5аба7а8; . . .; ^fc+1 ±am-2k + 1am-2k+2 . . . amam+1; Y в2 ,

2 ^ с ^ k + 1; в — 2?-1-значное целое число. Значение exp(xm) запишем как произведение

exp(xm) = exp(Y2) exp(Y3) ■ ... ■ exp(Yfc+1).

В методе быстрого вычисления экспоненты (методе Карацубы; англ. FEE [7]) величины exp(Y?) далее рассчитываются с помощью ряда Тейлора (3):

в в2 вГ

ехр(т?) = 1 + + 212^2 + • • • + Г22^ + (г) = & + (г) (4)

Здесь г = т2-?+1. Так как для остаточного члена выполняется неравенство [7]

|в?|г+1

IR(r)| < 2-

I

(г + 1)!22? (Г+1)’ то |Я?(г)| < 2-т. Величины получаются из формул

& = ат, = г!22? г,

к

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

5. Модификация метода двоичного деления для ряда Тейлора

логарифмической функции

Рассмотрим ряд Тейлора логарифмической функции [16]:

ОО х2 х3 хП

1п(1 + ж) = £ (-1)(г-1) - = х - — + — - ••• + (-1)(п-1) — + г„(х). (5)

i=1

i 2 3 n

Этот ряд сходится для значений х на полуинтервале (-1,1]. Остаточный член ряда (5) для |х| < 1 в форме Коши имеет вид [16]

ln+1 / 1 _ Q N n

|r,,(x)K 1-^ , где 0 <«< 1.

При х ^ 0 имеем 1 + 0х > 1 — и последний множитель меньше единицы. Поэтому для остаточного члена ряда на полуинтервале [0,1) получаем следующую оценку:

|х| -+1

Мх)| « 1—-Щ.

Пусть т = 2к+1, т ^ 25, где к — натуральное число. Будем вычислять ряд (5) с точностью 2-т в точке 7(?) = в(?)2—«, где £ = 2? — 1:

1п (1 + 7(?)) = — (Уг + • • • ± + д(?)(г), (6)

где с и г — натуральные числа, такие, что с ^ 1, г = т2(-?+2), в(?) — 2(?-1)-значное целое число. Ряд (6) линейно сходится относительно т, так как для остаточного члена на полуинтервале [0,1) выполняется следующая оценка:

І в (?)Г+1 22(?-1) (г+1) 2 2

|Д«(г)| < 2^________ < 22_________ = _____2____ = — < 2-(т+1)

1 ( )| 2?(г+1) 2?(г+1) 2?т2(-с+2) 22^ ,

то есть сходимость ряда (6) не зависит от с.

Модифицируем метод двоичного деления для частичной суммы ряда (6) так, чтобы вычисления находились в пределах класса Sch(FQLIN-TIME//LIN-SPACE).

Возьмём к1 = 1с^(г), г1 = [г/к1] (к1 —натуральное число) и запишем частичную сумму ряда (6) в виде

где

Р (7(?)) = т + Г2[^2 + тз[тз + • • • + 7кі_1[0А;і_1 + Ткіт^]]], (7

в(?) (в(?))2 , , (в(?))Г1-1

Т1 = “ТЗ---- —ТТТЗ + • • • ±

Т2

2« 2 • 2«2 ' (г 1 — 1)2«(гі-1) ’

(в «г

2«гі ’

Гі —1

т ,і вы , (вЬ)) (вы)

Т2 = ±— Т 7-----------, ±7-------, Т • • • ±

тз = Т2

г1 (г1 + 1)2« (г1 + 2)2«2 ’’ (2г1 — 1)2«(гі—1) ’

(в^ )Гі 2«гі ,

Гі —1

т і вм , (вй) т . (вы)

Тз = ± — Т 72------------------------------, -.4^ ± 72-, Т • • • ±

'2гі (2г1 + 1)2« (2г1 + 2)2«2 ’’ (3г1 — 1)2«(гі—1) ’

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

^(к) = г 1 - 1; м ил /(1 + ^ * = 0,

а(г) = ±1’ ВД=\№ - 1)Г1 + •■). *> 0; (8)

р(3 )= I1, 3 = 0’ 9(3 )= I1, 3 = 0’

) I вЧ 3 = 0; 9(3) 12«, з> 0.

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

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

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

Доказательство лемм 1 и 2 аналогично доказательству подобных лемм из [2, 3]. Рассчитывать приближённые значения Р (7(?)) с точностью 2-(т+1) по формуле (7) будем в соответствии со следующим итеративным процессом:

^1(Ш1) = ^,

^¿(Ш1) = окх-г+1 + Тс1-г+2^-Ъ г = 1, . . . , к1, (9)

^г(ш1) = ^г(ш1) + £»;

при г = к1 полагаем Р (7(?))* = Лк1 (т1). Здесь т1 ^ т, а* — приближения а с точностью 2-т1, Т*—приближения Т с точностью 2-т1. Величины ^г(т1) получаются отбрасыванием битов 5ш1+15ш1+2 ... 5т1+^ чисел Л,^(т1) после двоичной точки, начиная с (т1 + 1)-го, то есть

Ы = |^г(т1) - ^г(т1)| = 0,0 ... 0^т1+1^т1+2 . . . , (10)

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

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

|Л,Дт,1)| < 4. (11)

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

Ы < 1, Ы = (7(Т < 2.

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

(3 + 1) ^ 2:

|^^'+1(т1)| = К-0'+1)+1 + Т/с1-(^+1)+2^^ + е^'+1| < 1 +

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

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

Д(кьт1) < 2-т1+к1+3.

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

Н1 аЙ1, Нг аЙ1-г+1 + ТЙ1-г+2Нг—1, п(^ т1) |Лг(т1) Нг!.

_ I П-ТО1 2 + 2

4 + 2-т1 < 4.

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

n(1,mi) = |hi(mi) - Hi| = |о^ - afcl| < 2-mi+4.

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

n(j + 1,mi) = |^ki-(j+)+i + rfci-(j+i)+2hj(mi) + ei+i - °fci-(j+i) + i - rfci-(j+i)+2Hj 1 <

< |tUhj(mi) — tuhj(mi) + hj(mi) — Hj| + 2 ■ 2-mi ^ 2-mihj(mi) + 2-in(j,mi) + 2 ■ 2-mi.

Так как из (11) |hj(mi)| < 4 и, по индукционному предположению, n(j, mi) < 2-mi+j+3, то

n(j + 1,mi) < 4 ■ 2-mi + 2-i2-mi+j+3 + 2 ■ 2-mi < 2-mi+(j+i)+3.

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

Из леммы 4 следует, что достаточно взять mi = 2m + 3, чтобы вычислять P (y(?)) с точностью 2-(m+i).

Обозначим алгоритм расчета частичной суммы ряда (6), использующий схему (9), через ModifBinSplit (modified binary splitting, алгоритм 3).

Алгоритм 3. ModifBinSplit. Приближенное значение ряда Тейлора логарифмической

функции

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

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

1 mi := 2m + 3

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

3 Для всех г = 2, 3,..., к1

4 рассчитываем а := а£ _^+1 с точностью 2-т1 с помощью двоичного деления и Ь := т* ¿+2 с точностью 2-т1; обычного алгоритма

5 вычисляем выражение Л := а + Ь • Л;

6 Л присваиваем величину Л, округленную в соответствии с ;10);

7 на выход записываем Л.

Оценим временную вычислительную сложность алгоритма 3 при расчёте на машине Шёнхаге:

— 0(^(т)) вычислений а* дают 0(Мэсь(т) ^(т,)2);

— 0(^(т)) вычислений т дают 0(М8сь(т) к^(т)2);

— 0(log(m)) умножений чисел длины 0(т) дают 0(М8сь(т) log(m));

итого получаем 0(Мэсь(m)log(m)2). Емкостная сложность модифицированного алгоритма двоичного деления ModifBinSplit — это 0 (т), так как во всех вычислениях в данном алгоритме фигурируют числа длины 0(т).

Утверждение 1. Модифицированный алгоритм ModifBinSplit двоичного деления для вычисления ряда Тейлора логарифмической функции принадлежит классу алгоритмов 8сЬ(гдЬЩ-Т1МЕ//ЬЩ-8РАСЕ).

6. Быстрое вычисление функции 1п(1 + ж)

Пусть и = 1 + ж, а(х)(5) и а(и)(5) —полуинтервалы [2-5,1 — 2-5) и [1 + 2-5, 2 — 2-5) на вещественной прямой соответственно.

Построим метод вычисления функции 1п(1 + ж) на полуинтервале а(х)(5), аналогичный быстрому методу вычисления экспоненциальной функции, и назовём его, по аналогии с методом быстрого вычисления экспоненты, быстрым методом вычисления логарифмической функции.

Будем вычислять приближённое значение 1п(—0 )* с точностью 2-п, где ио = 1 + ж0, ио Е а(и)(5), п ^ 24. Возьмём натуральные к и т, такие, что

к ^ ш1п{3 : п ^ 2^}, т = 2к+1. (12)

Пусть ит — двоично-рациональное приближение аргумента и0 с точностью 2-т, то есть

|ит — и0| ^ 2-т. Так как ит Е [1, 2), то двоично-рациональная запись числа ит есть

_ 1 (1) (1) (1) (1) ит — 1,а1 а2 . . . ат+1.

(1) и(1)

Обозначим ит через и(1); возьмём г>(1) = 1,а1 ) и вычислим и(2) = —-т- с точно-

^(1)

и(1)

стью 2-т, то есть и(2) = —-т- + —(2), где |—(2) | ^ 2-т. Легко проверить, что первый

-у(1)

символ в двоичной записи величины и(2) —это 0, то есть и(2) = 1,0а22)а32) ... аетг)а(^+1. В силу соотношения 1п(и) = 1п ^ ^ = 1п ^^ + 1п(г>) имеем следующее равенство:

1п (и(1)) = 1п ^и^(1)^ = 1п (и(2) — —(2)) + 1п (^(1)) = 1п ^1 — —^л)^ + 1п (и(2)) + 1п (^(1)) .

Далее, возьмём ^(2) = 1,0а22)а32) и поделим и(2) на ^(2) с точностью 2-т. В результате получим соотношение

1п (и(1)) = 1п ^ — и2)^ + 1п ^ — и(3)^ + 1п (и(3)) + 1п (^(1)) + 1п (^(2)) .

Затем берём ^(3) = 1,000а43)а53)а63)а73) и выводим аналогичную формулу для 1п (и(1)). Возьмём натуральное число р = к + 2. На р-м шаге данного процесса, начиная с шага 2, получим следующую формулу для 1п (и(1)):

р / —(*) \ р-1 р р-1

1п (и(1)) = £ 1п ( 1 + —Г) ) + 1п (—(р)) + £ 1п (^(г)) = £ Е* + С + £ Я*.

г=2 \ и(г)/ г=1 г=2 г=1

Оценим сверху величины и С, используя неравенство | 1п(1 + ж)| < |ж| (данное

неравенство следует из свойств знакопеременных рядов [10]). Так как |и(г)| ^ 1, то —(г)

—— ^ 2-т и < 2-т. Далее, так как —(р) = 1 + ((р), где |£(р)| ^ 2-т, то С < 2-т.

и(г)

Получаем, что если вычислять величины Яг с точностью 2-т, то погрешность вычисления 1п (и(1)) оценивается следующим образом:

е1(т) = 1п (и(1))* — 1п (и(1)) < 4р • 2-т = (^(т) + 2)2-т+2.

Если взять т ^ 2п (что согласуется с (12)), то е1 (т) < 2-(га+2). Далее, так как

е2 = | 1п (и(1)) — 1п (и0)| = | 1п (и0 + С(1)) — 1п (и0)|

1М 1 + —

и0

< 2"

где |£(1)| ^ 2 т, то получаем оценку погрешности вычисления функции 1п(—) в точке и0:

1п (и(1))* — 1п(и0) ^ е1 + е2 < 2-(га+1).

Вычислять приближения величин Я?, то есть величины 1п ('У(?^ , будем с помощью алгоритма ModifBinSplit для ряда (6); для этого запишем ^(?) для 1 ^ с ^ (р — 1) в виде 1 + 7(?), где

7(?) = в М2-(2?-1),

(1)

в(i) = «1 , в(2)

(2) (2) «2 ^3 ),

в (з) = Л,(з)Л,(з)Л1 (3)^(3)

Здесь в(?) — 2(? 1) -значное целое число.

Алгоритм 4. FLE (fast logarithm evaluation). Быстрое вычисление логарифмической функции вещественного аргумента

Вход: Запись 0n точности вычисления 2-n. Оракульные функции: аргумента x

Выход: Приближённое значение ln(1 + x) с точностью 2-n на полуинтервале [2-5,1 — 2-5)

1: Возьмём к и m так, чтобы выполнялось (12);

2: p : = к + 2;

3: u(i) := <£x(m);

4: S := 0 (S — двоично-рациональное число).

5: Для всех i = 2, 3,... ,p

(i) u(i-i)

6: вычисляем u(i) = ;

v(i-i) ’

7: рассчитываем H(i) := ln (1 + 7(i)) с помощью алгоритма ModifBinSplit с точно-

стью 2-m;

8: S := S + H(i);

9: на выход записываем значение S, округленное до точности 2-n.

Вычислительная сложность алгоритма 4 при расчёте на машине Шёнхаге следующая: временная сложность — 0(Мэсь(п) log(n)3), так как алгоритм расчёта гипер-геометрических рядов ModifBinSplit использует 0(Мэсь(п) к^(п)2) операций, а в алгоритме ЕЬЕ используется 0(log(n)) таких вычислений и 0(log(n)) сложений чисел длины 0(п); емкостная сложность — 0(п), так как во всех вычислениях фигурируют числа такой длины.

Утверждение 2. Алгоритм ЕЬЕ быстрого вычисления логарифмической функции вещественного аргумента принадлежит классу вычислительной сложности на машине Шёнхаге 8сКЕдЬШ-Т1МЕ//ЬЩ-8РАСЕ).

Заключение

Алгоритм ЕЬЕ расчета логарифмической функции можно применять в информатике как основу Sch(FQLIN-TIME//ЬШ-8РАСЕ)-алгоритмической функции 1п(1 + ж),

заданной на подмножестве множества Sch(FQLIN-TIME//LIN-SPACE) алгоритмических вещественных чисел.

Отметим, что если для умножения использовать рекурсивный метод Карацубы [17] с временной сложностью O(nlog3), то временная сложность алгоритма расчёта логарифмической функции FLE равна O(nlog3 log n), то есть она полиномиальна степени меньше чем 2.

Для построения алгоритма расчёта логарифмической функции на произвольном промежутке можно использовать мультипликативную редукцию интервала [18] по аналогии с тем, как осуществляется мультипликативная редукция интервала для логарифмической функции в [3].

ЛИТЕРАТУРА

1. Яхонтов С. В. Вычисление гипергеометрических рядов с квазилинейной временной и линейной емкостной сложностью // Вестник Самарского государственного технического университета. Сер. Физико-математические науки. 2011. Вып. 2 (17). С. 239-249.

2. Яхонтов С. В. Эффективное по времени и по памяти вычисление экспоненциальной функции комплексного аргумента на машине Шёнхаге // Вестник С.-Петербург. унта. Сер. 10. Прикладная математика. Информатика. Процессы управления. 2011. Вып. 4. С. 97-110.

3. Яхонтов С. В., Косовский Н. К., Косовская Т. М. Эффективные по времени и памяти алгоритмические приближения чисел и функций: учеб. пособие. СПб., 2012. 256с.

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

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

6. Карацуба Е. А. Быстрое вычисление exp(x) // 17-я Всесоюз. школа по теории информации и её приложениям. Проблемы передачи информации. 1990. Т. 26. №3. С. 109.

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

8. Карацуба Е. А. Fast evaluation of hypergeometric function by FEE // Proc. 3rd CMFT conference on computational methods and function theory, Nicosia, Cyprus, October 13-17, 1997. Singapore: World Scientific, 1999. Ser. Approx. Decompos. V. 11. P. 303-314.

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

10. Карацуба Е. А. Быстрое вычисление ((3) // Проблемы передачи информации. 1993. Т. 29. №1. С. 68-73.

11. Karatsuba C. А. Fast evaluation of Bessel functions // Integral Transforms and Special Functions. 1993. V. 1. No. 4. P. 269-276.

12. Карацуба Е. А. Быстрое вычисление дзета-функции Римана ((s) для целых значений ар-

гумента s // Проблемы передачи информации. 1995. Т. 31. No. 4. С. 69-80.

13. Карацуба Е. А. Быстрое вычисление дзета-функции Гурвица и L-рядов Дирихле // Про-

блемы передачи информации. 1998. Т. 34. №4. С. 342-353.

14. Ахо А., Хопкрофт Дж., Ульман Дж. Построение и анализ вычислительных алгоритмов. М.: Мир, 1979. 536c.

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

16. Фихтенгольц Г. М. Курс дифференциального и интегрального исчисления. М.: Физмат-лит, 2003. Т. 2.

17. Карацуба А. А., Офман Ю. П. Умножение многозначных чисел на автоматах // Доклады АН СССР. 1962. Т. 145. №2. С. 293-294.

18. Muller J.-M. Elementary functions. Algorithms and implementation. Boston: Birkhauser, 1997. 204 p.