LINSPACE конструктивный аналог логарифмической функции Текст научной статьи по специальности «Математика»

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

Вестник СПбГУ. Сер. 10, 2008, вып. 2

LINSPACE КОНСТРУКТИВНЫЙ АНАЛОГ ЛОГАРИФМИЧЕСКОЙ ФУНКЦИИ

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

Основные сведения о классах вычислительной сложности, в частности о соотношениях классов предикатов Р и LINSPACE, можно найти в [2, 3]. FLINSPACE является стандартным обозначением для класса алгоритмов, использующих в процессе вычислений промежуточную память, объем которой ограничен сверху линейной функцией от длины входных данных (т. е. С-п). Под вычислениями в пределах LINSPACE будем понимать те, которые применяют алгоритмы из класса FLINSPACE. Через М(п) будем обозначать количество битовых операций, достаточных для умножения n-битовых целых чисел (для простого алгоритма это С ■ п2).

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

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

К алгоритмам, основанным на арифметико-геометрическом среднем, относятся алгоритмы, изложенные в работе [4], в которой показано, что некоторые константы математического анализа и большинство элементарных функций можно рассчитывать с точностью 2-п с временной сложностью 0(М(п) log2(n)) (возможно, что данная верхняя оценка является и нижней оценкой).

Алгоритмы, основанные на вычислении значений линейно сходящихся рядов с рациональными коэффициентами, предложены, например, в [5]. Данные ряды, в частности линейно сходящиеся гипергеометрические, используются для расчета констант математического анализа и элементарных функций в рациональных точках. Временная сложность алгоритмов из [5] оценивается сверху 0(M(n)(log2(n))2) (через п, как

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

ив [4], обозначается длина записи точности вычисления 2-п), а объем требуемой памяти пропорционален 0(nlog2(n)), так как быстрые алгоритмы вычисления рядов из [5] используют метод двоичного деления (binary splitting), требующий в классическом варианте 0(nlog2(n)) памяти для промежуточных результатов. В последнее время появились работы, например [6], в которых описываются алгоритмы (на основе модифицированного метода двоичного деления) вычисления значений линейно сходящихся гипергеометрических рядов с временной сложностью 0(М(n)(log2(n))2) и емкостной сложностью О(п) при определенных ограничениях на коэффициенты ряда. Имеется множество других способов расчета значений степенных рядов, например метод Кара-цубы, но, по всей видимости, они не являются вычислениями в пределах LIN SPACE.

В конструктивном математическом анализе аргументы конструктивных вещественных функций - конструктивные вещественные числа и задаются алгоритмами, определяющими рациональные приближения, а функции вычисляются на основе приближенных значений аргументов. Поэтому при оценке емкостной вычислительной сложности конструктивных функций нужно учитывать длины записей рациональных приближений аргументов [1]. На данный момент нами не найдены публикации других авторов, в которых бы содержался подобного рода анализ емкостной вычислительной сложности элементарных функций, в частности логарифмической. Например, в [4] указывается, что для расчета элементарных функций достаточно знать п + 0(log2(n)) битов аргумента, но оценки объема промежуточной памяти, необходимой для вычислений, не приводятся. В [5] содержатся общие доводы, каким образом можно использовать приведенные алгоритмы для расчета элементарных функций для вещественных аргументов, но нет анализа возникающих при этом погрешностей и не даны оценки объема промежуточной памяти, нужной для вычислений с приближенными значениями вещественных аргументов. Для функций ех, sin(a;), cos(a;) известен способ, с помощью которого двоичная запись приближения вещественного аргумента разбивается на группы по 2", = 1,2,..., битов. Затем ряд Тейлора функции вычисляется для каждой группы битов с помощью подходящего быстрого алгоритма, и общий результат получается как комбинация промежуточных результатов. Но для логарифмической функции такое разбиение аргумента не подходит, так как не ясно, как получить общий результат.

Для практических приложений важно уметь вычислять значения конструктивных элементарных функций на множестве конструктивных вещественных чисел с помощью алгоритмов из некоторого емкостного класса сложности, которые одновременно были бы полиномиальными по времени. В качестве класса емкостной сложности, важного для практических приложений, представляется разумным взять класс LIN SPACE. По аналогии с классом временной сложности Р, его можно назвать классом осуществимых вычислений с точки зрения требуемой для вычислений памяти. Подтверждением этой неформальной характеристики может служить то, что в настоящее время в работах по теоретической информатике для вычислений в пределах LINSPACE утвердился термин «space-efficient evaluation» [6].

В [7] предлагается алгоритм из класса FLINSPACE для построения конструктивных аналогов элементарных функций на основе разложений в ряд Тейлора. В [8] алгоритм из [7] используется для вычисления в пределах LIN SPACE значений логарифмической функции на всей области определения. В данной работе на основе модификации алгоритма из [7] и вычислительной схемы из [8] приведены в явном виде алгоритмы для вычисления с произвольной точностью значений логарифмической функции. Предлагается программная реализация этих алгоритмов на языке программирования C# в среде программирования Visual Studio 7.1.

2. Конструктивные вещественные числа и функции. В монографии [1] в качестве вычислительной используется модель, основанная на понятии машины Тьюринга, а основу представления конструктивных объектов составляет понятие вычислимой последовательности, сходящейся по Коши, ф : N —¥ I) (здесь и далее N - множество всех натуральных чисел, включая 0). Область значений I) такой последовательности является всюду плотным в К естественным подмножеством множества рациональных чисел. Для последовательности, сходящейся по Коши и задающей вещественное число х, требуют, чтобы выполнялось |ф(п) — х\ ^ 2~п для любого натурального п.

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

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

Последовательность ф : N ^ И двоично-рационально сходится к вещественному числу х, если для любого п £ N выполняется ргес(ф(п)) = п + 1 и |ф(п) — х\ ^ 2-п. Множество всех функций, двоично-рационально сходящихся к вещественному числу х, обозначается Вещественное число х называется СР конструктивным, если

содержит вычислимую функцию ф.

Вещественная функция /, заданная на отрезке [а, 6], называется конструктивной функцией на этом отрезке, если существует машина Тьюринга М с оракульной функцией такая, что для любого х £ [а, 6] и любой вычислимой функции ф £ СРХ функция ф, вычисляемая М с оракульной функцией ф, принадлежит СР^ху Фактически это означает, что для любой вычислимой функции ф £ СРХ и любого п £ N машина М последовательно вычисляет т £ N ш й £ И такие, что

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

• на ленту запроса оракульной функции записывается точность вычисления аргумента 2-,п;

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

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

р

Г

\ф(т)^х\^2-т, И^/(Ж)| ^2-”.

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

Определение 1 [9]. Вещественное число х называется LINSPACE конструктивным вещественным числом, если существует LINSPACE вычислимая функция ф£СРх.

Определение 2 [7]. Функция f : [a, b] ^ R называется LINSPACE конструктивной вещественной функцией на отрезке [а,Ъ], если для любого х £ [а, 6] существует LINSPACE вычислимая функция ф из CFj(x).

Множества LINSPACE конструктивных вещественных чисел и функций будем обозначать LINSPACEcf и LINSPACEc\a,b} соответственно. Индекс С[а, 6] в обозначении LINSPАСЕС{а,ь] обусловлен тем, что конструктивные функции являются непрерывными на всей области определения [1]. Под построением конструктивного аналога вещественной функции / на отрезке [а, 6] будем понимать описание алгоритма, вычисляющего двоично-рациональные приближения с произвольной точностью значений f(x) для х £ [а, Ь].

Под приближением t с точностью 2-п будем понимать двоично-рациональное число t* такое, что 11* — t\ ^ 2~п и prec(t*) = п + 1. Для модуля t* выполняются неравенства

П - 1*1^ Г - 2-'\ т. е. r|^|i|+2-'\

|i|-|i*| ^ ^2-", т. е. |**| ^ |*| — 2-”.

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

11* - t\ «С 11* I + 11** - t\ < 2-(n+1) + 2-(n+1) = 2-n. (2)

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

3. Вычисление значений степенного ряда. Рассмотрим степенной ряд

0G

S(x) = щх1 = a,Q + сi\X + сI2X2 + ■ ■ ■ + ai,xk + ... (3)

i=0

на отрезке [—е, е], в > 0, в ^ Pi гДе Р """"" РаДиУс сходимости. Пусть щ - LINSPACE конструктивные вещественные числа, х - произвольное LINSPACE конструктивное вещественное число из отрезка [—£>, у], а» - LINSPACE вычислимые функции из CFai, ф — LINSPACE вычислимая функция из CFX. Пусть также т — некоторое натуральное число, оц(т) - двоично-рациональные приближения коэффициентов щ с точностью 2~т, ф(т) — двоично-рациональные приближения аргумента х с точностью 2~т. В [7] рассматриваются частичные суммы Sk(x) = ^1; ич,х‘ ряда (3) с подстановкой а»(т) вместо коэффициентов щ и ф(т) вместо аргумента х:

к

Sk(rn) = У^аг(гп)ф(гп)1.

*=0

Для расчета приближенных значений S£ (т) таких выражений используется схема Горнера для полиномов с округлением промежуточных результатов:

h0(m) = ак(т),

hr(m) = a.k-r{'fn) + hr-i(m)<f)(m), r = l,...,k, (4)

hr(m) = hr(m) + er;

при r = к полагаем S%(m) = hk(m). Величины hr(m) получаются отбрасыванием битов Qm+iQm+'i ■ ■ -Qm+p чисел hr(m) после двоичной точки, начиная с (т + 1)-го, т. е.

|er| = \hr(m) — hr(m)\ = 0.0 ... 0qm+1qm+2 ... qm+p, (5)

а знак er совпадает со знаком hr(m) (ясно, что |ег| < 2~т). Затем в [7] оценивается

погрешность вычислений

Д(k,m) = |S*k(m) - Sk(x)\

для случая произвольного отрезка [—£, £>] из области сходимости ряда (3). В отличие от [7], приведем формулировки условий, при которых будет выполняться Д(к,т) < 2~к при определенных ограничениях на модуль коэффициентов и модуль аргумента. Пусть Q - такая константа, 0 < Q ^ 2-1, что для всех i

oiiQ + 2 \ (6)

Для вывода основной оценки понадобится следующее вспомогательное утверждение. Утверждение 1. Для любого г ^ 0, если т ^ 5, справедлива оценка

\hr(m)\ < 2. (7)

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

\ho(rn)\ = \ak(rn)\ «С Ы + 2-™ «С 2-1 + 2-5 < 2.

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

\hj+i(m)\ = |ak_u+1)(m) + hj{т)ф(т) + ej+1\ <

< |a*_a+1)| + 2-™ + \hj(m) \ (|Ж| + 2-™) + 2-™ <

<Q + 2-m + 2(Q + 2-5 + 2-™) + 2~m sC sC 3 ■ Q + 6 ■ 2-5 sC 3 ■ 2-1 + 6 ■ 2-5 < 2.

Теперь можно сформулировать утверждение относительно Д(й,т).

Утверждение 2. При выполнении условия утверждения 1 погрешность вычислений по схеме (4) оценивается как

A(k,rn) <2-т22(к + 1). (8)

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

с = maxk=Q\ai(m) - £ = \ф(т) - х\,

Hr(x) = ak^r + (afc_(r_i) H---h (a*-1 + akx)x ...)x,

H0(x)=ak, r](r,m) = \hr(m) — Hr(x)\.

Ясно, что С si 2_то, £ si 2~m. Воспользуемся методом математической индукции для г](г,т) по г. База индукции при г, равном 0:

7](0,т) = \hQ(m) - Н0(х)\ = \ак(т) - ак| ^ 2~т < 2_то22(0 + 1).

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

г](г + 1,т) = |аЛ_(г+1)(т) + Ьг(т)ф(т) + ег+1 - а*_(г+1) - Нг(х)х\ =

= | (afc_(r+1)(m) - afc_(r+i)) + er+i + {Кг{т)ф{т) - hr(m)x)+

+ (hr(m)x - Hr(x)x)| ^ С + kr+i| + IMm)l? + rj(r, m)g.

Так как из (7) \hr(m)\ < 2 и по индукционному предположению г](г,т) < 2-TO22(r + 1), то

Т](г + 1,т) < С + |er+i| + 2£ + 2_то22(г + 1)р ^

«С 2-т(2 + 2 + 2(r + 1)) < 2_т22(г + 2).

Из Sf,(m) = hk(m), Sk(x) = Hk(x), Д(k,m) = rj(k,m) получаем искомое неравенство.

Пусть требуется вычислить приближенное значение частичной суммы S^h) (х) РяДа (3) с точностью 2~к для некоторого натурального к, где ц - функция от к, ц(к) ^ 0. Положим

m = m(ki, к2) = 4 + |"log2(fci + 1)] + (к2 + 1), (9)

где к\ = /i(fc), к2 = к. Тогда выполняется условие утверждения 1, т. е. пг ^ 5, и,

следовательно, погрешность вычисления приближенного значения (х) будет оце-

ниваться в соответствии с (8):

A(fi(k),rn(kuk2)) < 2-т22{ц{к) + 1) < 2-к2 = 2-'\ (10)

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

Алгоритм Partial Sum. Приближенное значение Sk(x).

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

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

Выход: Приближенное значение S^^x) с точностью 2~к.

Описание:

1) вычисляем кг := ц(к), к2 := к, т по формуле (9);

2) вычисляем ф := ф(гп);

3) инициализируем цикл: h:=a{k\,m)-,

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

а) рассчитываем очередной коэффициент ряда а := а(к\ —г,т),

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

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

5) на выходную ленту записываем двоичное представление к (в этом представлении будет т битов после двоичной точки).

Пусть степенной ряд (3) линейно сходится, т. е. для остаточного члена выполняется |гс(ж)| = |5с(ж) — ^(ж)! ^ 2-(п+1) для любого х £ [—р, р], где ц - линейная функция, С = 1л(п + 1). Тогда

Возьмем т = m(ki, к2) по формуле (9) для к\ = ц(п + 1), к2 = п + 1. Тогда, используя неравенство (10), получаем

Отсюда видно, что линейной зависимости т от п достаточно для вычисления суммы степенного ряда (3) на отрезке [р, р] с требуемой точностью. По условию, коэффициенты ряда (3) являются LINSPACE вычислимыми, поэтому алгоритм Partial Sum, определяющий hr(m) по схеме (4), будет оперировать с числами, длины которых линейно зависят от п, т. е. данный алгоритм будет принадлежать классу FLINSPACE. Если коэффициенты ряда (3) являются также полиномиально вычислимыми по времени, то алгоритм Partial Sum будет в классе FLIN SPACE П FP (через FP обозначается класс алгоритмов, полиномиальных по времени).

4. Вычисление значений логарифмической функции. В [8] рассматривается ряд Тейлора логарифмической функции

Этот ряд сходится для значений х в интервале (—1,1]. С помощью редукции интервала вычисление логарифма для аргумента и = 1 + х можно свести [8] к расчету логарифма для аргумента £■:

здесь г - некоторое целое число. Ввведем следующие обозначения: и' = ^г, х* - приближение х, (и')* - приближение и'. Приближенные значения аргумента конструктивной функции являются двоично-рациональными числами, а такие числа имеют конечное двоичное представление. Поэтому для вычисления (и')* достаточно сдвигать битовое представление и* так, чтобы старший значащий бит оказался сразу после двоичной точки (т. е. (и')* = 0.1,...). Если и* ^ 1, то г будет такое, что г > 0, иначе г ^ 0. При этом (и1)* будет принадлежать интервалу [|-, 1), а величина (х1)*, равная (и1)* — 1, будет находиться в интервале 0).

Рассмотрим функцию 1п(1+ж') = 2-11п(1+ж'), ряд которой получается умножением коэффициентов ряда (11) на 2-1:

|Sc*(m) - S(x)\ «С |Sc*(m) - Sc(x)\ + \Sc(x) - S(x)\ sC

«С Д(С,m) + |rc(z)| sC A((,m) + 2^(n+1).

|S£(m) - S(x)\ < 2-(n+1) + 2-(,l+1) = 2-"

DO

(11)

DO

Так как (х')* £ [—1-,0), то х' £ [—^ — 5, 5) при условии, что х' вычисляется с

точностью 2~т, т ^ 5. Обозначим интервал [—| — 2 г'.2 г') через ст(5).

На основании формул из [8] для |гп(ж')|, х' < 0, можно получить оценку для остаточного члена ряда (12) на интервале [—^ — 2_5,0):

М*')1 «С «С

17'

¥

п+1

об

. ±_ < 2(п+1)(4.1-5)22 _ 2-0.9(п+1)+2

15

Для х' £ [0, 2””5) заведомо будет |гп(ж')| < 2~п. То есть ряд (12) линейно сходится на интервале а(5). Так как коэффициенты данного ряда ограничены сверху 2-1, то можно воспользоваться схемой (4) с Q = 2-1 для вычисления в пределах LIN SPACE данного ряда для х' на интервале ст(5), в частности для расчета ln(2-1) = 2-1 ln(l — 2-1). При этом в качестве /л(п + 1) достаточно взять выражение f 1.2гг. + 3], так как в этом случае будет |rM(n+i)(x')| < 2-°'9(1'2п+3+1)+2 < 2~(п+1\

Из (9), учитывая, что к\ = f 1.2гг. + 3], к2 = п + 1, получаем точность аргумента, достаточную для определения значения 1п(1 + х') с точностью 2~п:

2-т(кг,к2)^ ще т{кък2) = 4+ [~log2 {к\ + 1)] + (к2 + 1).

Легко убедиться, используя простой алгоритм, что операция обращения натурального числа является LINSPACE вычислимой, т. е. коэффициенты ряда (12) являются LINSPACE вычислимыми. Поэтому алгоритм вычисления значения функции 1п(1 + х') по схеме (4) будет принадлежать классу FLINSPACE. Очевидно также, что все вычисления будут полиномиальными по времени.

Приведем из [10] точности вычисления аргументов х, у арифметических операций, достаточные для получения результата с точностью 2-”: для сложения - это 2~1'п+1\ для умножения - это 2~(n+q+1\ где q таково, что \х\ + 1 < 29, |t/| + 1 < 2q. Далее, при умножении на 2 достаточно знать множитель с точностью 2-("+1) для расчета произведения с точностью 2-п, так как

| f — i| = |2 - (t')* - 2-t'\ = 2- | (t')* -t' | sC 2 ■ 2-(n+1) = 2-n,

здесь t = 2 ■ t', (t')* — приближенное значение t'. Следовательно, для вычисления величины исходной функции

ln(l + х) = ln(u) = ln(u') - г ■ ln^-1) = 2 ■ (ln(u') - г ■ ln(2-1))

с точностью 2_(n+1) и последующим округлением в соответствии с (2) до точности 2~п нужно знать значение слагаемых, т. е. ln(u') и г ■ 1п(2-1), с точностью 2-*"+3).

Рассмотрим слагаемое ln(u'). Величину и' нужно вычислять с точностью 2-”', где w находится по формуле (9) для

кг = |~1.2п + 6], к2 = п + 4, w = m(ki,k2) = 4 + [log2(fci +1)] + (к2 + 1). (13)

Введем обозначения:

и< = w + |г|, и> = тах(ш — |г|, 0).

Рассмотрим отдельно случаи, когда г<0иг>0(и' = и для г = 0). Если г < 0, т. е. происходит сдвиг значащих битов по направлению от младших разрядов к старшим, то х достаточно знать с точностью 2~'1><, так как

1 + X* 1+х

и =

2Г

2Г

= трк* -х\ «S 2-

Аналогично, если г > 0, то х достаточно знать с точностью 2~'0>. Значение V = тах(и<,г;>) = ш + |г| подойдет для обоих случаев знака г.

Здесь возникает некоторая сложность, состоящая в том, что для вычисления х* нужна величина сдвига, которую можно узнать только после расчета и* = 1 + х*. Поэтому предположим, что

х 5? 2~р - 1 и X «с 2Р - 2 (14)

для некоторого натурального р ^ 1. Тогда, используя неравенства (1) для х* = ф{г>), получаем следующие соотношения, которые выполняются при V ^ р + 1:

ф(у) > 2-р - 1 - 2-(р+1) = 2-(р+1) - 1 и ф(ь) «С2р-2 + 2-(р+1) <2^-1,

поэтому и* ;> 2_(г>+1) и и* < 2Р. Величина (и1)* получается сдвигом и* на г битов. Если г < 0, то 1 > (и1)* ;> 2-(г>+1)+1г1, и для г будет выполняться (р + 1) — |г| > 0, т. е. |г| < р + 1. В случае г > 0 будет 2-1 ^ (и1)* < 2Р~Г, и тогда р — г > —1, т. е. г < р + 1. Следовательно, всегда г такое, что |г| < р + 1 для р ) 1, и в качестве V достаточно взять т+р.В этом случае количество битов после двоичной точки величины (и1)* не будет превышать ш + 2р.

Далее, множители произведения г ■ 1п(2-1) нужно, в свою очередь, вычислять с точностью 2-("+3+«+1), и так как |1п(2-1)| + 1 < 2 и |г| + 1 < р + 2, то в качестве д можно взять

q = max(1, |"log2(p + 2)]) = |"log2(р + 2)]. (15)

Заметим, что величина q не зависит от х.

Учитывая свойства алгоритма Partial Sum, а также (13), (14), получаем следующее утверждение.

Утверждение 3. Для вычисления логарифмической функции 1п(1 + х) на отрезке [2~р — 1, 2^ — 2] с точностью 2~п в пределах LINSPACE достаточно знать приближенное значение аргумента х с точностью 2-(u'+p), где w = 4+ [log2(fci +1)] +(fc2 + l), ki = |"1.2n + 6], fc2 = n + 4. При этом объем промежуточной памяти будет ограничен сверху величиной С ■ (п + log2 (п) + р).

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

Алгоритм LnValue. Приближенное значение логарифмической функции. Вход: Запись точности вычисления 2~п.

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

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

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

Описание:

1) вычисляем w по формулам (13);

2) вычисляем v := w + р;

3) вычисляем ф := ф(t>);

4) вычисляем а : I + о:

5) определяем г:

а) если и < 1, то г ^ О,

б) если и ^ 1, то г > 0;

6) сдвигая битовое представление и на г позиций, влево или вправо в зависимости от знака г, получаем и'-,

7) вычисляем х' := и' — 1;

8) используем алгоритм Partial Sum для вычисления приближенного значения ряда (12) с точностью 2-("+3), т. е. а\ := 1п(1 + х');

9) находим q по формуле (15);

10) используем алгоритм PartialSum для расчета с точностью 2-("+3+«+1) приближенного значения <22 := 1п(2-1);

11) вычисляем произведение <23 := г ■ 1п(2-1) с точностью 2-*"+3);

12) вычисляем <24 := <2j — <23 с точностью 2-*"+2);

13) вычисляем <25 := 2 ■ <24 с точностью 2-*"+1);

14) на выходную ленту записываем двоичное представление <25, в котором после двоичной точки оставляем только (п + 1) бит.

Если для умножений и делений использовать алгоритмы из класса FLINSPACE с временной сложностью 0(М(п)), то временная сложность алгоритма LnValue будет 0(п ■ М(п + log2(n) + pj).

5. Программная реализация. Набор классов на языке программирования C# для вычисления значений вещественной логарифмической функции с произвольной точностью является частью библиотеки классов, позволяющей работать с LINSPACE конструктивными вещественными числами и функциями. Поэтому некоторое внимание будет уделено описанию общих интерфейсов и классов.

Для записи двоично-рационального числа вида Щ- используется структура DRNumber, состоящая из знака, массива битов числа m и натурального числа п. Все операции над двоично-рациональными числами возвращают результат в нормализованном виде: битовый вектор содержит только значащие биты числа т, при этом п интерпретируется как сдвиг массива битов m относительно двоичной точки вправо на п позиций.

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

public interface IComputableReal {

DRNumber CalcApprox(uint n);

}

абстрактный метод CalcApprox которого вычисляет двоично-рациональные приближения конструктивного числа с точностью 2-". Интерфейс

public interface IComputableFunc {

DRNumber CalcApprox(IComputableReal[] args, uint n);

}

обобщает понятие конструктивной вещественной функции.

Реализацией интерфейса IComputableReal является класс

public abstract class CompRealNumber : IComputableReal { public abstract DRNumber CalcApprox(uint n) public CompRealNumber() {}

}

Предусмотрены классы для конструктивных рациональных чисел, соответствующие целым числам, произвольным рациональным числам и рациональным числам вида 2к. Все такие классы являются производными от класса CompRealNumber. Реализация арифметических операторов, определенных для конструктивных вещественных чисел, находится в отдельных классах:

CompRealAdd, CompRealSubstr, CompRealMult, CompReallnverse,

CompRealDiv.

Абстрактный метод

public abstract DRNumber CalcApprox(IComputableReal[] args, uint n);

класса

public abstract class CompRealFunc : IComputableFunc { protected readonly Rationlnterval Domain; public CompRealFunc(Rationlnterval Domain)

}

возвращает аппроксимацию значения функции с точностью 2-п с набором аргументов в массиве args.

Приближенные значения частичной суммы S(х) степенного ряда (3) в точке х, являющейся конструктивным вещественным числом, вычисляются в соответствии с алгоритмом PartialSum, переопределенным методом

public override DRNumber CalcApprox(CompRealNumber[] args, uint k)

{•••}

класса

public abstract class PartialSum : CompRealFunc { protected int from; protected int to;

}

Величины from, to задаются в производных классах и означают начальный и конечный индексы частичной суммы степенного ряда, т. е. значение to равно ц(к). Величина т, обозначающая показатель точности 2-то, с которой нужно находить коэффициенты и аргумент степенного ряда, вычисляется по формуле (9):

uint m = 4+MiscCalc.UpperLog2(to+l)+(k+l);

Величина S^k) (х) рассчитывается в цикле с помощью методов

protected abstract DRNumber Calc_ai_approx(int i, uint m); protected abstract int Calc_ai_sign(int i);

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

for(int r=l; r<=to; r++) {

DRNumber ai_approx = Calc_ai_approx(to-r,m); int ai_sign = Calc_ai_sign(to-r);

DRNumber b = h*phi;

DRNumber h_hat = null;

if(ai_sign==l)

h_hat = ai_approx+b; else if(ai_sign==-l) h_hat = b-ai_approx; else if(ai_sign==0) h_hat = b;

h = DRNumberCalc.Round(h_hat,m);

}

Здесь ai - значение текущего коэффициента, hJiat - величина h(m) из схемы (4). У полученного значения h берется m битов после двоичной точки.

Класс LnHatValue, производный от абстрактного класса PartialSum, вычисляет функцию 1п(1 + х') = ijln(l + х') на интервале <т(5) с помощью переопределенных методов Calc_ai_approx и Calc_ai_sign. Величины логарифмической функции на отрезке [2~р — 1, 2^ — 2] определяются в соответствии с алгоритмом LnValue с помощью класса

public sealed class LnValue : CompRealFunc { private readonly uint p; private readonly LnHatValue ln_hat;

}

В конструкторе данного класса создается объект класса LnHatValue:

public LnValue(uint p) :base(null) { this.p = p;

ln_hat = new LnHatValue();

}

с помощью которого в методе CalcApprox, после редукции интервала, вычисляются двоично-рациональные приближения значения исходной функции. Для этого сначала определяется приближенное значение и с точностью 2-”', где w вычисляется по формулам (13). Затем с помощью сдвига и получается и':

uint v = 4+MiscCalc.UpperLog2(kl+l)+(k2+l)+p; DRNumber phi = x.CalcApprox(v);

DRNumber u = DRNumber.const1+phi; int r = u.To+1; DRNumber u_ = null;

if (r<0)

u_ = DRNumberCalc.ShiftLeft(u,-r); else if (r>0)

u_ = DRNumberCalc.ShiftRight(u,r); else if (r==0) u_ = u;

Здесь выражение u.To означает позицию в двоичной записи самого старшего ненулевого бита числа и. Значение 1п(1 + х') рассчитывается с помощью объекта lnJiat:

IComputableReal х_ = new DRNumberWrap(u_-DRNumber.const1);

IComputableReal[] ln_hat_args = new IComputableReal[]{x_};

DRNumber ln_hat_x_ = ln_hat.CalcApprox(ln_hat_args,n+3);

Для вычисления выражения ln(2-1) также используется объект lnJiat:

IComputableReal _1_2 = CRNumberOrder.Neg(

new CompRealRational(DRNumberConstr.Whole(1),DRNumberConstr.Whole(2))); ln_hat_args = new IComputableReal[]{_1_2}; uint q =

MiscCalc.UpperLog2(p+2); DRNumber ln_hat_l_2 = ln_hat.CalcApprox(ln_hat_args,n+3+q+l); DRNumber d =

DRNumberConstr.Whole(r)*ln_hat_1_2;

Затем на основе полученных величин находится значение исходной функции:

DRNumber res = ln_hat_x_-d; res = DRNumberConstr.Whole(2)*res; res = DRNumberCalc.Round(res,n+1);

Алгоритм LnValue был опробован на компьютере со следующей конфигурацией: процессор AMD Athlon 64 Х2 Dual Core 3800+, системная плата MSI MS-7250, память 2Gb DDR2-SDRAM 670 MHz. Приложение строилось в конфигурации Release с включенной оптимизацией кода. Для вычисления приближенного значения 1п(5) с 3320 двоичными знаками (приблизительно 1000 десятичных знаков) после точки понадобилось 27 с. При этом для умножения двоично-рациональных чисел использовался простой алгоритм с временной сложностью С ■ п2.

6. Заключение. Алгоритм LnValue (вместе с базовым алгоритмом PartialSum) можно применять в информатике как основу LINSPACE конструктивной вещественной логарифмической функции, заданной на множестве LINSPACE конструктивных вещественных чисел. Предложенные алгоритмы не являются оптимальными по времени выполнения, но для построения LINSPACE вычислимого математического анализа важна принципиальная возможность построения линейных по памяти алгоритмов.

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

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

Summary

Yakhontov S. V. The LINSPACE constructive analogue of the logarithm function.

The algorithms for evaluation of the real logarithm function are given. It is proved the algorithms are in the class FLINSPACE and they are also polynomial-time algorithms. The implementation of the algorithms using C# programming language and Microsoft Visual Studio 7.1 IDE is described.

Литература

1. Ко К. Complexity theory of real functions. Boston: Birkhauser, 1991. 310 p.

2. Du D., Ко К. Theory of computational complexity. New York: John Wiley & Sons, 2000. 491 p.

3. Косовская Т. М., Косовский H. К. Различия между классом LINSPACE и рядом других классов сложности // Тез. докл. XIV Междунар. конференции «Проблемы теоретической кибернетики». Пенза, 23-28 мая 2005 г. М.: Изд-во мех.-матем. факультета Моск. ун-та, 2005. С. 72.

4. Brent R. Fast multiple-precision evaluation of elementary functions // J. of the ACM. 1976. Vol. 23, N 2. P. 242-251.

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. Cheng H., G erg el B., Kim E., Zima E. Space-efficient evaluation of hypergeometric series j j ACM SIGSAM Bull. Communications in Computer Algebra. 2005. Vol. 39~ N 2. P. 41-52.

7. Яхонтов С. В. Вычисление степенных рядов в пределах LINSPACE // Материалы VIII Междунар. семинара «Дискретная математика и ее приложения». Москва, 2-6 февраля 2004 г. М.: Изд-во мех.-матем. факультета Моск. ун-та, 2004. С. 118-122.

8. Яхонтов С. В. Вычисление логарифмической функции в пределах LINSPACE для конструктивных вещественных чисел // Материалы IX Междунар. семинара «Дискретная математика и ее приложения». Москва, 18-23 июня 2007 г. М.: Изд-во мех.-матем. факультета Моск. ун-та, 2007. С. 149-151.

9. Яхонтов С. В. Ь/Уб'РДС'-Е-конструктивность алгебраических чисел // Труды V Междунар. конференции «Дискретные модели в теории управляющих систем». Ратмино, 2629 мая 2003 г. М.: Изд. отдел факультета ВМиК Моск. ун-та, 2003. С. 96-97.

10. Шурыгин В. А. Основы конструктивного математического анализа. М.: Едиториал УРСС, 2004. 328 с.

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

Статья принята к печати 4 декабря 2007 г.