Mirror
/
cpython
mirror of https://github.com/python/cpython
4
1
Fork 0
Code Issues Releases Wiki Activity

1) List relative error for _mpd_qln10(). 

2) Add rigorous error analysis to _mpd_qlog10 (ACL2 proofs exist).

3) Use the relative error as a basis for the interval generation in the
   correction loop (same as in _mpd_qln()).
This commit is contained in:
Stefan Krah 2012-06-08 18:41:33 +02:00
parent ed36b2e55b
commit 1cf6dfc8b2
1 changed files with 29 additions and 9 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

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

@ -4298,7 +4298,14 @@ static const mpd_t _mpd_ln10 = {
(mpd_uint_t *)mpd_ln10_data (mpd_uint_t *)mpd_ln10_data
}; };
/* Set 'result' to ln(10). ulp error: abs(result - log(10)) < ulp(log(10)) */ /*
* Set 'result' to log(10).
* Ulp error: abs(result - log(10)) < ulp(log(10))
* Relative error : abs(result - log(10)) < 5 * 10**-prec * log(10)
*
* NOTE: The relative error is not derived from the ulp error, but
* calculated separately using the fact that 23/10 < log(10) < 24/10.
*/
void void
mpd_qln10(mpd_t *result, mpd_ssize_t prec, uint32_t *status) mpd_qln10(mpd_t *result, mpd_ssize_t prec, uint32_t *status)
{ {
@ -4712,21 +4719,34 @@ mpd_qln(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
} }
} }
/* Internal log10() function that does not check for specials, zero, ... */ /*
* Internal log10() function that does not check for specials, zero or one.
* Case SKIP_FINALIZE:
* Relative error: abs(result - log10(a)) < 0.1 * 10**-prec * abs(log10(a))
* Case DO_FINALIZE:
* Ulp error: abs(result - log10(a)) < ulp(log10(a))
*/
enum {SKIP_FINALIZE, DO_FINALIZE};
static void static void
_mpd_qlog10(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx, _mpd_qlog10(int action, mpd_t *result, const mpd_t *a,
uint32_t *status) const mpd_context_t *ctx, uint32_t *status)
{ {
mpd_context_t workctx; mpd_context_t workctx;
MPD_NEW_STATIC(ln10,0,0,0,0); MPD_NEW_STATIC(ln10,0,0,0,0);
mpd_maxcontext(&workctx); mpd_maxcontext(&workctx);
workctx.prec = ctx->prec + 3; workctx.prec = ctx->prec + 3;
/* relative error: 0.1 * 10**(-p-3). The specific underflow shortcut
* in _mpd_qln() does not change the final result. */
_mpd_qln(result, a, &workctx, status); _mpd_qln(result, a, &workctx, status);
/* relative error: 5 * 10**(-p-3) */
mpd_qln10(&ln10, workctx.prec, status); mpd_qln10(&ln10, workctx.prec, status);
workctx = *ctx; if (action == DO_FINALIZE) {
workctx.round = MPD_ROUND_HALF_EVEN; workctx = *ctx;
workctx.round = MPD_ROUND_HALF_EVEN;
}
/* SKIP_FINALIZE: relative error: 5 * 10**(-p-3) */
_mpd_qdiv(NO_IDEAL_EXP, result, result, &ln10, &workctx, status); _mpd_qdiv(NO_IDEAL_EXP, result, result, &ln10, &workctx, status);
mpd_del(&ln10); mpd_del(&ln10);
@ -4807,9 +4827,9 @@ mpd_qlog10(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
prec = ctx->prec + 3; prec = ctx->prec + 3;
while (1) { while (1) {
workctx.prec = prec; workctx.prec = prec;
_mpd_qlog10(result, a, &workctx, status); _mpd_qlog10(SKIP_FINALIZE, result, a, &workctx, status);
_ssettriple(&ulp, MPD_POS, 1, _ssettriple(&ulp, MPD_POS, 1,
result->exp + result->digits-workctx.prec-1); result->exp + result->digits-workctx.prec);
workctx.prec = ctx->prec; workctx.prec = ctx->prec;
mpd_qadd(&t1, result, &ulp, &workctx, &workctx.status); mpd_qadd(&t1, result, &ulp, &workctx, &workctx.status);
@ -4829,7 +4849,7 @@ mpd_qlog10(mpd_t *result, const mpd_t *a, const mpd_context_t *ctx,
mpd_del(&aa); mpd_del(&aa);
} }
else { else {
_mpd_qlog10(result, a, &workctx, status); _mpd_qlog10(DO_FINALIZE, result, a, &workctx, status);
mpd_check_underflow(result, &workctx, status); mpd_check_underflow(result, &workctx, status);
} }
} }