CC BY
58
УДК 519.683
С. Ю. Соловьев1
АЛГОРИТМ ВЫЧИСЛЕНИЯ ЛОГАРИФМОВ МЕТОДОМ ВЫТЕСНЕНИЯ
Предлагается метод разложения вещественных чисел на сомножители специального вида, каждый из которых в свою очередь допускает аналогичное разложение. Описывается связь метода с двоичным представлением вещественных чисел. Обсуждаются вопросы программной реализации метода для задачи вычисления натурального логарифма. Рассматриваются подходы к обобщению метода.
Ключевые слова: логарифм, метод, редукция, погрешность, алгоритм.
1. Введение. Любое положительное число X можно однозначно представить в нормализованном виде 2Р?7, где Р целое, а и вещественное из полуинтервала [0.5,1); Р называется порядком, и — мантиссой. Применительно к логарифму это означает существование представления 1п(Х) = 1п(?7) — Р1п(0.5). Иначе говоря, для вычисления любого логарифма достаточно знать алгоритм его вычисления на [0.5,1).
2. Метод вытеснения. Приведем цепочку формальных свойств, обосновывающих способ редукции (см. [1]) аргумента — формулу (1), применимую для вычисления значений некоторых элементарных функций. Свойства относятся к числам из полуинтервала [0.5,1) и в силу своей простоты приводятся без доказательств. Прежде всего введем специальные обозначения для трех числовых последовательностей:
An. —
Br, =
1
— Ai, Сп —
■1'
где п = 1,2,... .
Последовательность Ап порождает семейство нумерованных непересекающихся полуинтервалов • • • вида 1п = [Ап-1, Ап). В совокупности эти полуинтервалы образуют разбиение (рис. 1) множества [0.5,1), т. е. [0.5,1) = /2 и /3 и ... . Для полуинтервала [Ап, 1) будем использовать обозначение:
1п+ = 1п +1 и 1п+2 и... .
Рис. 1. Семейство полуинтервалов /2,/з,
Свойство 1. Ап_ 1 < Вп < Сп < Ап для всех п ^ 2 (рис. 2, а).
Любое число х из [0.5,1) принадлежит ровно одному полуинтервалу /п, а значит, существует однозначная функция г: [0.5,1) —)► {2, 3,...}, принимающая значение п, если х Е 1п• Например, ¿(0.625) = 2, ¿(0.828125) = 3.
1 Факультет ВМК МГУ, проф., д.ф.-м.н., e-mail: [email protected]
Рис. 2. Свойства 1, 2 и 3 как основания свойства 4 и формулы (1)
Пусть х Е 1г = Дг) и [Вг,Аг). Рассмотрим число х/Аг.
Свойство 2. Если х Е то х/Ах Е [С^, А^) С (рис. 2, б).
Свойство 3. Если х Е то х/А% Е 1) = (рис. 2, б).
Как следует из свойств 2 и 3, преобразование полуинтервала Iг посредством деления его чисел на приводит к двойственному результату:
— числа из [Вг,Аг) вытесняются в полуинтервал а
— числа из Востаются в но свойство 1 гарантирует, что повторное деление на все-таки вытеснит их в
Свойство 4. Если х Е то справедливо равенство
^ = {(х/А1)АхАх, если ж < ](х/Аг)Аг — иначе,
причем Ах Е 1, х/А2х Е для х < Вг и х/Ах Е — иначе. Логарифмируя правую и левую части равенства (1), имеем
Свойство 5. Если х Е то справедливы следующие равенство и связанные с ним соотношения:
1п(ж) = / 1п(жМ*) + 21п(Аж), если ж < Да, 1 + — иначе;
Аг Е
если х < Вг, то х/А^ Е > . (3)
если х ^ то х/Аг Е 1г+ J
Другими словами, формула (2) сводит вычисление 1п(ж) для либо к паре подзадач 1п(Аг) (для Е 1) и (для х/А\ Е
либо к паре подзадач 1п(Аг) (для Ах Е /2+1) и (для х/Ах Е
С одной стороны, при таком подходе количество задач удваивается, а с другой стороны, все аргументы вытесняются в Вытеснение аргументов заводит их в левую окрестность единицы, где поведение функции 1п(ж) хорошо изучено.
Свойство 6. 1п(ж) ~ х — 1 при х 1.
Свойство 7. Если 0.5 ^ 1 — 5 ^ ж < 1, то 0 < ж - 1 - \п(х) ^ щз^у ~~ I)2-
Не ограничивая общности, будем рассматривать случай 5 = 2_7?, где 77 — целочисленный параметр точности, 77 = 2, 3,... . При этом
Свойство 8. Условия 1 — 2~Г] ^ х < 1, А^ ^ х < 1, х Е и 2(ж) > 77 эквивалентны. Для аргумента ж из [0.5,1) и положительного параметра точности 77 рекурсивно определим функцию Шп(ж; 77) следующим образом:
ж — 1, если 2 > 77, где 2 =
Шп(ж;т/) = + 2Шп(Лг; 7/), если (4)
Шп(ж/А^;т7) +Шп(А^;т7) — иначе.
Идея вычисления уо = 1п(жо) методом вытеснения состоит в том, что уо « Шп(жо;т7). Параметр 77 считается известным. Алгоритмическую корректность, т.е. конечность процесса вычисления Шп(жо;т7), гарантируют свойства 5 и 6.
"Врожденным пороком" метода является замена 1п(ж) на х — 1 при г(х) > 77. Общая погрешность вычисления складывается из ошибок округления при выполнении арифметических операций
Да и из погрешности замен Де, т.е. |Шп(жо;7?) — уо| < Да + Де. Для оценки Де определим (из самых общих соображений) верхнюю границу числа замен. Наибольшее количество замен возможно, если:
— исходный аргумент .r(, е- / >. и
— одна задача для аргумента из I2 порождает две задачи для аргументов из I3,
— две задачи для аргументов из I3 порождают четыре задачи для аргументов из /4 и т. д.
В конечном итоге на полуинтервале 1п+ заменами решаются 24-1 задач. С учетом свойства 7 погрешность Де оценивается так:
д ^ -1-/2-т2 . 2„-1 = -1-2~(т1+2) < 2~(ч+1).
2(1 — 2"Ч)v ! 1^2 '>
Предложенный в [2] метод редукции аргумента функции log2(a;) позволит улучшить оценку на два порядка: Де < Вместе с тем привлечение дополнительного свойства Ап/А^п_£ /•_>,< • •_> позво-
ляет привести полученную оценку к виду Де < где m(r]) — рекуррентно определенная
функция поправок, принимающая значения 0, 0, 1, 1, 1, 2, 2, 3, 3, 4, 4, 5, 5, 6, 7, 7, 8, 9, 9, 10, 11, 11, 12, 13, 13, 14, 15, 16, 16, 17, 18, 18, 19, 20, 21, 22, 22, 23 для г] = 2, 3,..., 38,40 - 1. Так, при г] = 15, т(г]) = 6 и Де < 2"162"6 = 2"22.
3. Программирование метода вытеснения. Прежде всего отметим три обстоятельства, важные для практической реализации метода.
1. Задача нормализации чисел требует знания стандартов представления вещественных чисел в компьютере [3]. Не отвлекаясь на эти детали, будем считать известной процедуру FEX(X, Р, U), которая по заданному положительному X вычисляет его порядок Р и мантиссу U. Результатом FEX(1, Р, U) являются Р = 1 и U = 0.5, результатом FEX(5.5, Р, U) — Р = 3и11 = 0.6875 и т.д. Кроме того, будем считать фиксированным количество разрядов в двоичном представлении мантиссы. Для определенности будем считать мантиссу сорокаразрядной, полагая, что в настоящем изложении число 40 играет роль параметра.
2. Функция z(x) удобна для записи формул, но необходимость вычисления z(x) определяет переборный характер метода вытеснения. Безусловно, наилучший способ вычисления z(x) состоит в исследовании разрядов двоичного представления х. Фактически z(x) есть номер самого левого нуля в дробной части аргумента: z(0.625) = z(0.10l2) = 2, z(0.828125) = z(0.11010l2) = 3 и т.д. Однако из методологических соображений для вычисления z(x) будем использовать заранее подготовленный массив-справочник вещественных чисел ABD[1,..., 40] = (Ai,..., Л40). В этом случае z(X) = Z, если ABD[Z - 1] < X < ABD[Z].
Так или иначе, вычисление z(x) сводится к перебору номеров 2, 3,... до тех пор, пока не встретится подходящий номер z, а после его обнаружения возникают две новые задачи z(Az) и z{x/Akz) (для к = 1 или к = 2), но, как следует из свойства 5, для их решения процесс перебора номеров можно начинать с z + 1.
3. Окончание вычислительного процесса определяется условием Z > ЕТА, где ЕТА — целочисленная константа. Например, ЕТА = 15.
Буквальная реализация формулы (4) на языке Паскаль состоит из вспомогательной рекурсивной функции RCR и из основной функции LNR:
function RCR(X : real; Z : integer) : integer;
begin Z:=Z+1;
if Z > ETA
then RCR:=X-1
else if X < sqr(ABD[Z])
then RCR:=RCR(X/sqr(ABD[Z]),Z) + 2*RCR(ABD[Z],Z) else if X < ABD LZ]
then RCR:=RCR(X/ABD[Z],Z) + RCR(ABD[Z],Z)
(* next Z *) else RCR:=RCR(X,Z) end;
function LNR(X : real) : real; var P : integer;
begin FEX(X,P,X); LNR:=RCR(X,1)-P*RCR(0.5,1) end;
Первая из них выполняет всевозможные проверки и преобразования для аргумента X на полуинтервале Iz+1, а вторая обеспечивает корректные обращения к RCR и вырабатывает искомое приближение 1п(Х).
Даже при разовом вычислении LNR(X) многократно возникают и решаются "с нуля" совершенно одинаковые задачи вида RCR(ABD[Z], Z). В связи с этим имеет смысл заранее запасти результаты решения этих задач в справочнике ALN[ 1,..., 40]. Кроме того, для чисел Bz также имеет смысл завести справочник BBD[1,...,40], т.е. ALN[Z] = RCR(ABD[Z], Z), BBD[Z] = Bz. В этом случае формула (4) принимает следующий вид:
Г] Г]
Rln (х-Г]) =xY[A~ki -1 + ^kiX ALN[i], где к{ Е {0,1,2}.
г=2 г=2
Соответствующий алгоритм выглядит так:
function LNT(X : real) : real; var T : real;
Z : integer; begin FEX(X,Z,X);
T:=-Z*LNA[1]; for Z:=2 to ETA do
if X < BBD [Z] then begin X:=X/BBD[Z]; T:=T+2*ALN [Z] end else if X < ABD[Z] then begin X:=X/ABD[Z]; T:=T+ ALN[Z] end; LNT:=X-1+T end;
Временная сложность алгоритма LNT есть О (77), а для вычисления первых 77 значений справочника ALN можно предложить алгоритм с оценкой сложности 0(г]2). При переходе от LNR к LNT точность вычислений остается неизменной. Заметим, что для фиксированного параметра точности 77 погрешность вычисления LNT(X) можно уменьшить, если внести в справочник ALN вместо Kln(Az;r]) существенно более точные значения In(Az), полученные, скажем, иным методом вычисления логарифмов
[1, 4].
В приведенном алгоритме источниками погрешности Аа являются данные из справочников ALN и BBD, а также операции деления; остальные действия допускают точное выполнение. Общая проблема всех реализаций связана с вычислением значений Bz- Из-за ограниченности разрядов для записи мантиссы при Z > 40/2 числа Bz теряют младшую единицу и становятся неотличимыми от Az~ 1, поэтому условия X < BBD[Z] оказываются ложными, и целая ветвь вычислений оказывается заблокированной. Таким образом, при ЕТА >40/2 алгоритмы начинают существенно зависеть от погрешностей округления, что ограничивает их применение. И все-таки ситуация не столь драматична.
4. Обобщение метода. Предложенный в п. 2 метод вытеснения можно рассматривать как представителя целого класса методов, отличающихся друг от друга параметрами. Рассмотрим один из способов параметризации, порождающий обобщенный метод вытеснения, схематично представленный на рис. 3, а. Формально обобщенный метод описывается системой ограничений на значения параметров Еп и Яп :
' Ап_ 1 ^ Еп < АП1 что соответствует Еп Е /п;
An-i
Ап ^ П А ^ Я А
Ап ^ Rn < 1, Еп
An < -7- < 1,
-/in
■ ^ 1, [An-i,En)/(RnAn) С In+; Rn E In+'i [En,An)/An С In+.
Рис. 3. Схема обобщенного метода и область ERr
Эквивалентными преобразованиями эта система сводится к неравенствам
Вп ^ Еп ^ Сп ,
Еп/Ап ^ Rn ^ Сп/Ап. ^
Пары (ЕП1 Дп), удовлетворяющие (5), образуют треугольную область ERn, изображенную на рис. 3,
Свойство 9 (обобщение свойства 5). Если х Е Iz и (EZ,RZ) Е ERZl то справедливы формула (6) и соотношения (7):
1п(ж) = + ln(A*) + 1п(Д*), если ж < Я*,
1 — иначе;
^ G 4+1 и Д* G
если ж < EZ1 то (x/Rz)/Az Е > . (7)
если х ^ EZ1 то x/Az Е J
При фиксированном способе выбора Еп и i?n вычисление логарифма с помощью обобщенного метода вытеснения задается формулой
х — 1, если 2 = 2(ж) > г/;
Шт(ж;?7) = < Мт((а;/Дг)/Лг;г)) + Шт(Дг;?7) + Мт(Дг; 77), если z ^ 77 и ж < Шт(ж/А^;?7) +R\m(Az;r]) — иначе.
При некоторых значениях Еп и Rn обобщенный метод вытеснения приобретает дополнительные свойства, способные повлиять и на область применимости алгоритма, и на выбор справочников. В частности, все ERn содержат пары (Efn,An+1), где Е'п = (Ап_ 1 + Ап)/2 — середина полуинтервала 1п. Использование Е'п и Ап+i в качестве значений Еп и Rn порождает алгоритм, который зависит от ошибок округления при ЕТА >40 — 2, что вполне приемлемо для реальных задач.
В общем случае области ERn обладают сеточным строением, как показано на рис. 4 и сформулировано в виде свойства 10.
Рис. 4. Строение областей ERr
Свойство 10. Для любого п ^ 2 справедливы неравенства
Ап < Ап+1 < ... < Аг < Сп/Ап,
где г = г(Сп/Ап), и
Вп = АпАп < т4гг_|_1т4гг < ... < АгАп < Сп.
Таким образом, узлами сетки для области ЕРп являются:
1) пары (АкАп, Ат), где (/г, т) Е {п, п + 1,..., г}2 и к ^ т (случай 1);
2) пары {А^Ап, Сп/Ап), где ^ = п, п + 1,..., г (случай 2);
3) пара (Сп,Сп/Ап) (случай 3).
Если в формуле (6) на полуинтервале 4 значения и Rz выбраны равными соответственно Az+pAz и А^+д (случай 1, п = г, к = г + р, т = г + д), то соотношения (7) принимают следующий уточненный вид:
Аг £ 4 + 1 ИЙ^Е + если Дз ^ ж, то ж/А* Е > •
если Ег > ж, то (х/11г)/Аг Е + )
Другими словами, задача вычисления 1п(ж) (для х € Iz) сводится к подзадаче \yl(Az) (для Az € Iz+1), которая решается:
— либо совместно с подзадачей In(x/Az) (для x/Az G /(г+р)+);
— либо совместно с подзадачами In((x/Rz)/Az) (для (x/Rz)/Az G /(г+р)+) и In(Rz) (для Rz G /(г+в)+).
Только одна из перечисленных подзадач — подзадача вычисления \yl(Az) — вытесняется в следующий полуинтервал Iz+1, в то время как остальные подзадачи вытесняются вперед нар и q полуинтервалов. На первый взгляд Еп и Rn следует выбирать так, чтобы числа р и q были максимальными, а еще лучше использовать Еп = Cn, Rn = Сп/Ап. Однако при существенных отклонениях параметров Еп и Rn от медианных значений алгоритм начинает зависеть от погрешностей округления.
5. Заключительные замечания. 1. Почти все операции, задействованные в методе вытеснения, реализуются быстрыми командами логического типа [2]. Единственное исключение составляет операция деления вещественных чисел, но и она допускает определенную свободу в выборе делителя, что позволяет контролировать ошибку округления. 2. Вычислительные эксперименты подтверждают работоспособность и эффективность метода. 3. Метод вытеснения допускает замену степеней двойки на степени чисел из [2, 3). Варьирование основаниями степеней способно сократить объем вычислений, но исключает использование команд логического типа. 4. Метод вытеснения можно успешно применять и для вычисления квадратных корней.
СПИСОК ЛИТЕРАТУРЫ
1. Brent R., Zimmermann Р. Modern Computer Arithmetic. Cambridge University Press, 2011.
2. Люстерник JI.A., Абрамов A.A., Шестаков В.И., Шура-Бура М.Р. Решение математических задач на автоматических цифровых машинах. М.: Изд-во АН СССР, 1952.
3. Overton М. L. Numerical Computing with IEEE Floating Point Arithmetic. Philadelphia: SIAM, 2001.
4. Сальников M. С. Рекурсивный алгоритм вычисления логарифма // Информационные процессы. 2012. 12. № 3. С. 248-252.
Поступила в редакцию 12.09.12
ALGORITHM LOGARITHM EVALUATION BY DISPLACEMENT METHOD
Soloviev S. Yu.
We propose a method of expanding the real numbers to the factors of a special type, each of which, in turn, allows for a similar expansion. We describe a method with the binary representation of real numbers. We discuss the software implementation of the method for the problem of calculating the natural logarithm. We consider approaches to the synthesis method.
Keywords: logarithm, method, reduction, precision, algorithm.
