From c62bd13cb203f27bc11d36903554f20386b5b8c5 Mon Sep 17 00:00:00 2001 From: Stefan Krah Date: Sat, 16 Jun 2012 19:45:35 +0200 Subject: 1) State the relative errors of the power functions for integer exponents. 2) _mpd_qpow_mpd(): Abort the loop for all specials, not only infinity. 3) _mpd_qpow_mpd(): Make the function more general and distinguish between zero clamping and folding down the exponent. The latter case is currently handled by setting context->clamp to 0 before calling the function. 4) _mpd_qpow_int(): Add one to the work precision in case of a negative exponent. This is to get the same relative error (0.1 * 10**-prec) for both positive and negative exponents. The previous relative error for negative exponents was (0.2 * 10**-prec). Both errors are _before_ the final rounding to the context precision. --- Modules/_decimal/libmpdec/mpdecimal.c | 20 ++++++++++++++++++-- 1 file changed, 18 insertions(+), 2 deletions(-) diff --git a/Modules/_decimal/libmpdec/mpdecimal.c b/Modules/_decimal/libmpdec/mpdecimal.c index 262e834..1f58203 100644 --- a/Modules/_decimal/libmpdec/mpdecimal.c +++ b/Modules/_decimal/libmpdec/mpdecimal.c @@ -5844,6 +5844,12 @@ mpd_qnext_toward(mpd_t *result, const mpd_t *a, const mpd_t *b, /* * Internal function: Integer power with mpd_uint_t exponent. The function * can fail with MPD_Malloc_error. + * + * The error is equal to the error incurred in k-1 multiplications. Assuming + * the upper bound for the relative error in each operation: + * + * abs(err) = 5 * 10**-prec + * result = x**k * (1 + err)**(k-1) */ static inline void _mpd_qpow_uint(mpd_t *result, const mpd_t *base, mpd_uint_t exp, @@ -5880,6 +5886,12 @@ _mpd_qpow_uint(mpd_t *result, const mpd_t *base, mpd_uint_t exp, /* * Internal function: Integer power with mpd_t exponent, tbase and texp * are modified!! Function can fail with MPD_Malloc_error. + * + * The error is equal to the error incurred in k multiplications. Assuming + * the upper bound for the relative error in each operation: + * + * abs(err) = 5 * 10**-prec + * result = x**k * (1 + err)**k */ static inline void _mpd_qpow_mpd(mpd_t *result, mpd_t *tbase, mpd_t *texp, uint8_t resultsign, @@ -5899,7 +5911,8 @@ _mpd_qpow_mpd(mpd_t *result, mpd_t *tbase, mpd_t *texp, uint8_t resultsign, if (mpd_isodd(texp)) { mpd_qmul(result, result, tbase, ctx, &workstatus); *status |= workstatus; - if (workstatus & (MPD_Overflow|MPD_Clamped)) { + if (mpd_isspecial(result) || + (mpd_iszerocoeff(result) && (workstatus & MPD_Clamped))) { break; } } @@ -5914,7 +5927,9 @@ _mpd_qpow_mpd(mpd_t *result, mpd_t *tbase, mpd_t *texp, uint8_t resultsign, } /* - * The power function for integer exponents. + * The power function for integer exponents. Relative error _before_ the + * final rounding to prec: + * abs(result - base**exp) < 0.1 * 10**-prec * abs(base**exp) */ static void _mpd_qpow_int(mpd_t *result, const mpd_t *base, const mpd_t *exp, @@ -5932,6 +5947,7 @@ _mpd_qpow_int(mpd_t *result, const mpd_t *base, const mpd_t *exp, workctx.round = MPD_ROUND_HALF_EVEN; workctx.clamp = 0; if (mpd_isnegative(exp)) { + workctx.prec += 1; mpd_qdiv(&tbase, &one, base, &workctx, status); if (*status&MPD_Errors) { mpd_setspecial(result, MPD_POS, MPD_NAN); -- cgit v0.12