http://hg.python.org/cpython/rev/fabc593a5822 changeset: 77268:fabc593a5822 user: Stefan Krah <skrah@bytereef.org> date: Thu May 31 20:01:05 2012 +0200 summary: Improve Underflow handling in the correct-rounding loop. The case for Underflow to zero hasn't changed: _mpd_qexp() internally uses MIN_EMIN, so the result would also underflow to zero for all emin > MIN_EMIN. In case digits are left, the informal argument is as follows: Underflow can occur only once in the last multiplication of the power stage (in the Horner stage Underflow provably cannot occur, and if Underflow occurred twice in the power stage, the result would underflow to zero on the second occasion). Since there is no double rounding during Underflow, the effective work precision is now 1 <= result->digits < prec. It can be shown by a somewhat tedious argument that abs(result - e**x) < ulp(result, result->digits). Therefore the correct rounding loop now uses ulp(result, result->digits) to generate the bounds for e**x in case of Underflow. files: Modules/_decimal/libmpdec/mpdecimal.c | 21 ++++++++++++-- 1 files changed, 17 insertions(+), 4 deletions(-) diff --git a/Modules/_decimal/libmpdec/mpdecimal.c b/Modules/_decimal/libmpdec/mpdecimal.c --- a/Modules/_decimal/libmpdec/mpdecimal.c +++ b/Modules/_decimal/libmpdec/mpdecimal.c @@ -4122,6 +4122,8 @@ MPD_NEW_STATIC(ulp, 0,0,0,0); MPD_NEW_STATIC(aa, 0,0,0,0); mpd_ssize_t prec; + mpd_ssize_t ulpexp; + uint32_t workstatus; if (result == a) { if (!mpd_qcopy(&aa, a, status)) { @@ -4135,12 +4137,20 @@ prec = ctx->prec + 3; while (1) { workctx.prec = prec; - _mpd_qexp(result, a, &workctx, status); - _ssettriple(&ulp, MPD_POS, 1, - result->exp + result->digits-workctx.prec); + workstatus = 0; + + _mpd_qexp(result, a, &workctx, &workstatus); + *status |= workstatus; + + ulpexp = result->exp + result->digits - workctx.prec; + if (workstatus & MPD_Underflow) { + /* The effective work precision is result->digits. */ + ulpexp = result->exp; + } + _ssettriple(&ulp, MPD_POS, 1, ulpexp); /* - * At this point: + * At this point [1]: * 1) abs(result - e**x) < 0.5 * 10**(-prec) * e**x * 2) result - ulp < e**x < result + ulp * 3) result - ulp < result < result + ulp @@ -4148,6 +4158,9 @@ * If round(result-ulp)==round(result+ulp), then * round(result)==round(e**x). Therefore the result * is correctly rounded. + * + * [1] If abs(a) <= 9 * 10**(-prec-1), use the absolute + * error for a similar argument. */ workctx.prec = ctx->prec; mpd_qadd(&t1, result, &ulp, &workctx, &workctx.status); -- Repository URL: http://hg.python.org/cpython