[Python-checkins] cpython: 1) State the relative errors of the power functions for integer exponents.
stefan.krah
python-checkins at python.org
Sat Jun 16 19:46:31 CEST 2012
http://hg.python.org/cpython/rev/73df491612aa
changeset: 77470:73df491612aa
user: Stefan Krah <skrah at bytereef.org>
date: Sat Jun 16 19:45:35 2012 +0200
summary:
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.
files:
Modules/_decimal/libmpdec/mpdecimal.c | 20 +++++++++++++-
1 files changed, 18 insertions(+), 2 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
@@ -5844,6 +5844,12 @@
/*
* 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 @@
/*
* 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 @@
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 @@
}
/*
- * 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 @@
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);
--
Repository URL: http://hg.python.org/cpython
More information about the Python-checkins
mailing list