A related case, now that I've actually looked at the implementation:
> (2251799813685577e0 * 2e0**-51) - 1
1.46105350040671e-13
> sprintf("%.13f", 2251799813685577e0 * 2e0**-51)
1.0000000000002
The Num value in this case is exactly 2251799813685577/2251799813685248,
or 1.000000000000146105350040670600719749927520751953125. (Subtracting
1 correctly shows a few of these digits.) The correct rounding to 13
decimal places is 1.0000000000001. There is plenty of precision left
in the float format after these decimal places: one unit in that last
output place corresponds to over 450 ulp in the float format.
What's going on is that the implementation is first converting to
decimal with its default of 15 total digits, getting a correctly-rounded
1.00000000000015, and then rounding *that* to the output length,
getting 1.0000000000002. This is double rounding, a classic type of
mistake in numerical computation. You need to perform rounding only
once, straight to the output length: either during the binary->decimal
conversion, or after a binary->decimal conversion that produced all the
digits necessary to represent the value exactly.
-zefram
