Here's the error analysis, for the record. It's crude, but effective. Assume
our decimal context precision is p (so p = 20 for the above code).
1. d is represented exactly with value in [1.0, 2.0). (Conversions from float
to Decimal are exact.)
2. c.ln(d) < 1.0, and Decimal's ln operation is correctly rounded, so the
*absolute* error in c.ln(d) is at most 0.5 * 10**-p.
3. Similarly, LOG2 has an absolute error of at most 0.5 * 10**-p.
4. Let q be the unique integer satisfying 10**q <= n < 10**(q+1).
5. r * LOG2 < n < 10**(q+1), so ulp(r * LOG2) <= 10**(q + 1 - p),
and the rounding error in computing r * LOG2 is at most
0.5 * 10**(q + 1 - p). We also inherit the error from LOG2, scaled
up by r (which is less than 10**(q+1)), so this inherited error
is also bounded by 0.5 * 10**(q + 1 - p). Our total absolute error
is now bounded by 10**(q + 1 - p).
6. The result of the addition of d to r * LOG2 is also bounded by
n < 10**(q+1), so introduces a new rounding error of up to
0.5 * 10**(q + 1 - p) on top of the errors in d and r * LOG2
(steps 2 and 5). Our total error is now bounded by 1.5*10**(q + 1 - p).
7. Division by n produces a result <= 1.0. It divides our absolute error so
far by n (>= 10**q), giving a new error of at most 15 * 10**-p, and
introduces a new rounding error of at most 0.5 * 10**-p. Bound
so far is 15.5 * 10**-p.
8. The exponential operation gives a result in the range [1.0, 2.0],
and converts our absolute error of 15.5 * 10**-p into a relative
error of at most 15.5 * 10**-p (let's call it 16 * 10**-p to be
safe), which since our result is bounded by 2.0 amounts to an absolute
error of at most 32 * 10**-p. We get a rounding error on the result
of at most 5 * 10**-p (like ln, the decimal module's "exp" operation
is also correctly rounded), making our total error bounded by 37 * 10**-p.
Finally, we want to express the error bound above in terms of ulps of the
original IEEE 754 binary64 format that float (almost certainly) uses.
Since our result (before the ldexp) is in the range [1.0, 2.0], one ulp
for the binary format is 2**-52. So the error bound in step 8, expressed
in binary64 ulps, is 37 / 2 **-52 * 10**-p < 2 * 10**(17-p) ulps. Rounding
from the final Decimal value to float also incurs an error of at most
0.5 ulps, so we end up with an error bound of
final_error <= 0.5 + 2 * 10**(17-p) ulps.
For p = 20, this is 0.502 ulps. |