CC BY
88
2014 Вычислительные методы в дискретной математике №3(25)
УДК 510.25+510.52+519.688
ВЫЧИСЛЕНИЕ ВЕЩЕСТВЕННОЙ W-ФУНКЦИИ ЛАМБЕРТА W0
В ПРЕДЕЛАХ FP//LINSPACE
М. А. Старицын, С. В. Яхонтов
Санкт-Петербургский государственный университет, г. Санкт-Петербург, Россия E-mail: [email protected], [email protected]
Строится FP//LINSPACE алгоритмический аналог вещественной W-функции Ламберта W0(x) на отрезке [-(re)-1, (re)-1] FP//LINSPACE алгоритмических вещественных чисел, где r — рациональное, r > 4/3 (в качестве r можно брать любое рациональное с таким условием). Для построения алгоритмического аналога вещественной W-функции Ламберта Wo(x) предлагается алгоритм WLE расчёта двоично-рациональных приближений данной функции на отрезке [-(re)-1, (re)-1] с полиномиальной временной и линейной емкостной сложностью на машине Тьюринга. Алгоритм WLE строится на основе разложения в ряд Тейлора данной функции, при этом показывается и используется в алгоритме линейная сходимость ряда Тейлора W-функции Ламберта W0(x) на отрезке [-(re)-1, (re)-1].
Ключевые слова: вещественная W-функция Ламберта W0, алгоритмические вещественные функции, машина Тьюринга, полиномиальная временная сложность, линейная емкостная сложность.
Введение
В работе предлагается алгоритм расчёта вещественной W-функции Ламберта W0 [1] на отрезке [-(re)-1, (re)-1], где r — рациональное, r > 4/3 (говоря более точно, основной ветви W0 вещественной W-функции Ламберта W) с полиномиальной временной и линейной емкостной сложностью на машине Тьюринга. Алгоритм строится на основе разложения в ряд Тейлора данной функции с использованием алгоритма вычисления линейно сходящихся степенных рядов в пределах FP//LINSPACE из [2] в качестве базового алгоритма.
При построении алгоритмического аналога вещественной W-функции Ламберта W0 берётся модель алгоритмических чисел и функций, изложенная в [3]. Посредством FP//LINSPACE будем обозначать класс алгоритмов, полиномиальных по времени и линейных по памяти при вычислении на машине Тьюринга [4].
W-функция Ламберта W является трансцендентной функцией и относится к классу специальных функций (т. е. как специальная функция W-функция Ламберта W не выражается через элементарные функции). W-функция Ламберта W интересна как с практической точки зрения, так как используется, например, в математических задачах физики, так и с теоретической точки зрения, например при рассмотрении вычислительной сложности констант и функций математического анализа в теоретической информатике.
Имеется достаточно большое количество алгоритмов расчёта констант и функций математического анализа: алгоритмы на основе AGM [5], метод Карацубы быстрого вычисления экспоненциальной функции [6], метод binary splitting [7] как вариант метода Карацубы, алгоритм bit-burst расчёта голономных функций (D-finite на англ.) [8], метод Ньютона вычисления обратных функций и др. Но на данный мо-
мент не известно никаких результатов относительно использования этих методов вычисления констант и функций математического анализа для вычисления вещественной W-функции Ламберта W0 на отрезке [-(re)-1, (re)-1 ] с линейной памятью на машине Тьюринга. Кроме того, некоторые из перечисленных методов не могут быть применены для расчёта W-функции Ламберта W0; например, алгоритм bit-burst расчёта голономных функций [8] неприменим к расчету W-функции Ламберта W0, так как данная функция не принадлежит к классу голономных функций. Поэтому результат, изложенный в данной работе, является новым в области вычислительной сложности алгоритмических чисел и функций.
Частично результат работы представлялся на конференции СПИСОК-14 [9].
1. Алгоритмические вещественные числа и функции
В данной работе основу представления конструктивных объектов (чисел и функций) составляет понятие алгоритмической последовательности <^, сходящейся по Коши [3], при этом в качестве вычислительной модели берётся машина Тьюринга. Такая последовательность определяется на множестве всех натуральных чисел N, включая 0, а областью аппроксимирующих значений является всюду плотное в R естественное подмножество множества рациональных чисел. Для последовательности, сходящейся по Коши и задающей вещественное число х, требуют, чтобы выполнялось
|^(n) — х| ^ 2-n
для любого натурального n.
В качестве множества аппроксимирующих значений берётся множество двоичнорациональных чисел D [3]. Рациональное число d называется двоично-рациональным, если d = m/2n для некоторого целого m и натурального n. Двоично-рациональные числа имеют конечное двоичное представление: строка s, равная
±UpUp-1 ... u0.v1 v2 ... vr,
обозначает число
d = ± ( £ Ui2i + ¿ Vj2-j
y=0 j=1
Длина представления двоично-рационального числа определяется как количество символов в строке s, равное, с учётом знака и двоичной точки, p+r + 3, и обозначается / (s). Под точностью представления prec(s) понимается число битов справа от двоичной точки, то есть r. С точки зрения изучения вычислительной сложности двоично-рациональные числа удобны тем, что для любого n двоично-рациональные числа с точностью n равномерно распределены на вещественной прямой [3].
Последовательность ^ ^ D двоично-рационально сходится к вещественному
числу х, если для любого n G N выполняется prec(^(n)) = n +1 и
|<^(n) — х| ^ 2-n.
Множество всех функций, двоично-рационально сходящихся к вещественному числу х, обозначается CFx.
Вещественное число х называется CF-алгоритмическим [3], если содержит
вычислимую функцию <^.
Вещественная функция f, заданная на отрезке [а, 6], называется алгоритмической функцией [3] на этом отрезке, если существует машина Тьюринга М с оракульной функцией, такая, что для любого х Е [а, 6] и любой вычислимой функции ^ Е СЕх функция ф, вычисляемая М с оракульной функцией <^, принадлежит CFf (х).
Фактически это означает, что для любой вычислимой функции ^ Е СЕж и любого п Е N машина М последовательно вычисляет т Е N и d Е О, такие, что
|^(ш) - х| ^ 2-т, |d - f (х)| ^ 2-п.
Сложность расчёта двоично-рациональных приближений алгоритмических чисел и функций определяется в [3] на основе длины двоичного представления точности вычисления. Память ленты запроса и ленты ответа оракульной функции при оценке емкостной вычислительной сложности алгоритма не учитывается. Обращение к оракульной функции ^ Е СЕх аргумента х алгоритмической функции осуществляется следующим образом:
— на ленту запроса оракульной функции записывается точность вычисления аргумента 2-т в виде 0т (унарная запись);
— рассчитывается значение ^(т) оракульной функции, и результат записывается на
ленту ответа;
— значение ^(т) считывается с ленты ответа в промежуточную память.
Определение 1 [2]. Число х Е К назовём ЕР//ЬШ8РАСЕ алгоритмическим вещественным числом, если существует функция ^ Е СЕх, вычислимая в пределах
ЕР//ЬШ8РАСЕ.
Определение 2 [2]. Вещественную функцию f, заданную на отрезке [а, 6], назовём ЕР//ЬШ8РАСЕ алгоритмической вещественной функцией на отрезке [а, 6], если для любого х Е [а, 6] функция ф (указанная в определении алгоритмической функции) из CFf(х) является РР//ЬШ8РАСЕ вычислимой.
Множества РР//ЬШ8РАСЕ алгоритмических вещественных чисел и функций будем обозначать ЕР//ЬШ8РАСЕс^ и ЕР//ЬЩ8РАСЕС[а6] соответственно. Здесь использование индекса С [а, 6] обусловлено тем, что алгоритмические функции являются непрерывными на всей области определения [3]. Построение алгоритмического аналога вещественной функции f на отрезке [а, 6] означает описание алгоритма, вычисляющего двоично-рациональные приближения с произвольной точностью значений f (х) для х Е [а, 6].
2. Вычисление W-функции Ламберта Ш0
Вещественная Ш-функция Ламберта Ш определяется как решение функционального уравнения
х = Ш (х)е^(х).
Данное решение является функцией, обратной к функции f (х) = хех. Вещественная Ш-функция Ламберта Ш определена на полуинтервале [—е-1, то) и имеет две ветви, верхнюю Ш0 и нижнюю Ш-1.
Будем строить алгоритм расчёта верхней ветви Ш0, который обозначим через ШЬЕ, на основе разложения в ряд Тейлора данной функции с использованием алгоритма вычисления линейно сходящихся степенных рядов в пределах РР//ЬШ8РАСЕ из [2] в качестве базового алгоритма.
ГО
Напомним, что степенной ряд Б = ГО а называется линейно сходящимся степенно
Мк)
ным рядом, если его частичная сумма 5^(&) = ГО а^, такая, что ^(к) —линейная функ-
г=0
ция от к, отличается от точного значения не более чем на 2-к: |Б — 5^)1 ^ 2-к.
Рассмотрим ряд Тейлора функции Ж0 [1] в окрестности точки х = 0:
го го (_к)(к-1)
^о(х) = ГО айХк = ГО----и----Xй, (1)
й=1 й=1 к-
радиусом сходимости данного ряда является величина е-1. Перепишем ряд (1) в виде
Ш01)(х) = Е ак1)хк = Е —-¡—тг-(ех)к , (2)
к=1 к=1 е ' к-
рассмотрим данный ряд для х Е [— (ге)-1, (ге)-1 ] и оценим сверху модуль п-го остатка
Е
дП1)(х) = Е ак)хк (3)
к=п+1
данного ряда. В силу неравенства к! ^ 22 ■ кк+1/2е к, следующего из формулы Стирлинга [10], получаем
кк
|а(1)тк| < ______к______ < 2-с1^к
1 к 1 ек ■ 22 ■ кк+1/2е-к (ге)к
для х Е [—(ге)-1, (ге)-1]. Здесь С1 —константа, зависящая от рационального г. Следовательно,
Е
д11)(х) < Е 2-Сгк = С2 ■ 2-С1(га+1),
к=п+1
то есть ряд (2) линейно сходится.
Далее, покажем ЕР//ЬШ8РАСЕ вычислимость коэффициентов а^ ряда (2) (отметим, что входными данными для алгоритма вычисления коэффициентов а^ является двоичная запись точности вычисления 2-т, что является точностью вычисления аргумента х [2]).
(1) (1) 1 к-1 Для этого запишем коэффициенты ак ) в виде произведения ак ) = ^П(-1)^-,
ек ^—1
к р
где 6^- = —т; обозначим а(р,к) = П (—1)^"-16^-. Будем вычислять величину а(к-1,к) (рав-е^ ^=1
ную ак1)) последовательно в цикле дляр Е {1,..., к—2}, на каждом шаге выполняя произведение а*рк)6р+1 (при этом а(1,к) = 61), где а*рк) —приближённое значение величины а(р,к) с некоторой точностью Ер; 6р+1 —приближённое значение величины 6р+1 с той же точностью Ер = 2-9 < 2-1; q — некоторое натуральное. Произведение ( = а*рк)6р+1 будем округлять с точностью Ер, то есть отбрасывать биты числа £ после двоичной точки, начиная с д-го бита.
Используя метод математической индукции по р Е {1,... , к — 2}, покажем, что если взять е 1 ^ 2-Сзк+(1°§2(к)+2), где Сз — некоторая константа, то
-сзк+ ^ (1°§2 (к)-1о§2 СЛ+2)
Ер ^ 2 ^
для любого р Е {1,... , к — 2}. База индукции: р = 1; в этом случае величина а*1 к), равная 6(, вычисляется с точностью Е1. Индукционный переход: пусть для р Е {1,... , к—3} выполняется |а*рк) — а(р,к)| ^ Ер. Тогда
|а(р+1,к) — а(р+1,к)1 ^ |а(р,к)6р+1 — а(р,к)6р+11 + Ер = |а(р,к)6(+1 — а(р,к)6р+1 + а(р,к)6р+1 — а(р,к)6р+11 + Ер ^
к
^ |а(р,к)(6(+1 — 6р+1)| + |6р+1(а(р,к) — а(р,к))| + Ер < 2ер + р + 1 Ер < 21°й(к) 1°§2(р+1)+2Ер
Р+1
1 -Сзк+ Й (1°§2(к)-1°§2(Л + 2)
(здесь используется оценка а(р,к) < 2-1), то есть Ер+1 < 2 ^'=1 .
Теперь оценим сверху величину
V = Е (1°§2 (к) — 1о§2(^)).
.7=1
В силу неравенства р! ^ 2рР+1/2е-р, следующего из формулы Стирлинга, получаем
применив тот факт, что функция f(х) = х(^2(2к) — ^2(х)) возрастает на отрезке [1, к].
В результате имеем оценку Ер+1 < 2-Сбк. Так как к ^ п (где п берётся из формулы (3)), то для вычислений а(р,к) можно взять точности Ер+1, такие, что Ер+1 < 2-СбП. Это означает ЕР//ЬШ8РАСЕ вычислимость коэффициентов а^ ряда Тейлора (2) функции Ш0(1)(х) (так как т линейно зависит от п для линейно сходящихся степенных рядов, а т такое, что 2-т — точность вычисления аргумента х [2]).
Так как ряд Тейлора (2) линейно сходится на отрезке [—(ге)-1, (ге)-1], где г — рациональное, г > 4/3, то воспользуемся алгоритмом беггезбит^ из [2] для вычисления данного ряда в пределах ЕР//ЬШ8РАСЕ. Алгоритм беггезбит^ позволяет вычислять с полиномиальным временем и линейной памятью на машине Тьюринга линейно схо-
Е
дящиеся степенные ряды вида Б(х) = ^ а»хг на любом отрезке а С [—д, д] при условии,
г=0
что |а^| ^ 1, д ^ 3/4 + 2-5 и величины ак являются ЕР//ЬШ8РАСЕ вычислимыми. Так как все условия, при которых алгоритм беггезбит^ применим, выполняются для ряда Тейлора функции Ш0(1)(х) (а значит, и для ряда Тейлора функции Ш0), то верна следующая теорема.
Теорема 1. Основная ветвь Ш0 вещественной Ш-функции Ламберта является ЕР//ЬШ8РАСЕ алгоритмической вещественной функцией на любом отрезке [—(ге)-1, (ге)-1] ЕР//ЬШ8РАСЕ алгоритмических вещественных чисел, где г — рациональное, г > 4/3.
Заключение
Алгоритм ШЬЕ расчета вещественной Ш-функции Ламберта Ш0 можно применять в информатике как основу ЕР//ЬШ8РАСЕ алгоритмической вещественной Ш-функ-ции Ламберта Ш0, заданной на отрезке [—(ге)-1, (ге)-1] ЕР//ЬШ8РАСЕ алгоритмических вещественных чисел, где г — рациональное, г > 4/3.
Из дальнейших исследований стоит отметить задачу построения алгоритмов, основанных на разложении в ряды, для FP//LINSPACE алгоритмических аналогов других вещественных функций, не выражающихся через элементарные функции.
ЛИТЕРАТУРА
1. Дубинов А. Е., Дубинова И. Д., Сайков С. К. W-функция Ламберта и ее применение в математических задачах физики. Саров: Изд-во ФГУП <РФЯЦ-ВНИИЭФ>, 2006. 160 с.
2. Яхонтов С. В., Косовский Н. К., Косовская Т. М. Эффективные по времени и памяти алгоритмические приближения чисел и функций. Учеб. пособие. СПб.: Изд-во СПбГУ, 2012. 256с.
3. Ko K. Complexity Theory of Real Functions. Boston: Birkhauser, 1991. 310 p.
4. Du D. and Ko K. Theory of Computational Complexity. N.Y.: John Wiley & Sons, 2000. 491 p.
5. Brent R. P. Fast multiple-precision evaluation of elementary functions // J. ACM. 1976. V. 23. No. 2. P. 242-251.
6. Карацуба Е. А. Быстрые вычисления трансцендентных функций // Проблемы передачи информации. 1991. Т. 27. Вып. 4. С. 76-99.
7. Haible B. and Papanikolaou T. Fast multiprecision evaluation of series of rational numbers // Proc. Third Intern. Symposium on Algorithmic Number Theory, Portland, Orgeon, USA, June 21-25, 1998. P. 338-350.
8. Mezzarobba M. A note on the space complexity of fast D-finite function evaluation // Computer Algebra in Scientific Comput. 2012. V. 7442. P. 212-223.
9. Старицын М. А., Яхонтов С. В. Эффективное по времени и памяти вычисление W-функции Ламберта // Четвертая Всерос. науч. конф. по проблемам информатики СПИСОК-14. СПб., 2014 (в печати).
10. Фихтенгольц Г. М. Курс дифференциального и интегрального исчисления. Т. 2. М.: Физ-матлит, 2003. 680 с.
