382
l
P. T. P. Tang
Step 2. Let T, be the largest working-precision number smaller than
e -l/16 - 1, and let 7’Z be the smallest one that is larger than e”‘6 + 1.
Then the inequality e-“‘6 - 1 < X < e1’16 - 1 is equivalent to
Ti < X < T2. The exact IEEE representations for T, and T2 are as follows:
Precision
T, Tz
single
BD782A03 3D8415AC
double
BFAF0540 438FD5C4 3FB082B5 77D34ED8
Step 3. Because 1 + X is not representable in general and may even overflow
(if the prevalent rounding mode is toward +m), we must calculate m, F,
and f carefully.
-Let T3 be a threshold chosen so that 1 + T3 does not overflow in any
rounding modes and
X = round-to-nearest(X + 1)
whenever X 2 T3. Thus, for example, T3 2 226 for single precision
and T3 2 255 for double. Calculate Y as follows: If X < T3, then
Y : = 1 + X. Otherwise, Y : = X.
-Determine the unique integer m such that 1 I 2-“Y < 2. Since Y is
always normalized, m can be obtained by logb( Y).
-Determine F as follows:
Y := 2-” * Y . . . exact; save Y for later use
F := round-to-nearest integer(27 * Y)
F := 2-7 e F
. . . exact
-Computation off requires much care. Let d be the number of signifi-
cant bits in working precision (24 for single and 53 for double).
Compute
f
by one of the three methods below, depending on m’s
value.
Case 1. m 5 -2. Compute f by
f:=
Y-F.
Case 2. -1 I m I d - 1. Compute
f
by
f :=
(2-” - F) + 2-” * X.
Case 3. m 1 d. Compute
f
by
f :=
(2-” * X - F) + 2-“.
Note that the computations ensure that
2m(F +
f) =
1 + X
except when m 2 d + 8. But when that is the case,
1 + X = 2m(F +
f) + a,
lal 5 1.
Thus,
log(1 + X) = log(2”(F +
f )) +
log(1 + 2-m~/(F +
f ))
= log(2”(F +
f )) + p, ( p 1 5
1.1 x 2-“.
ACM Transactions on Mathematical Software, Vol. 16, No. 4, December 1990.