traverseda
/
cpython
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

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.
This commit is contained in:
Stefan Krah 2012-06-16 19:45:35 +02:00
parent f185226244
commit c62bd13cb2
1 changed files with 18 additions and 2 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

20
Modules/_decimal/libmpdec/mpdecimal.c
View File

@ -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);