diff options
author | zimmerma <zimmerma@280ebfd0-de03-0410-8827-d642c229c3f4> | 2001-10-26 15:32:23 +0000 |
---|---|---|
committer | zimmerma <zimmerma@280ebfd0-de03-0410-8827-d642c229c3f4> | 2001-10-26 15:32:23 +0000 |
commit | 9452541d31eb789a64f9cb963200eddf57c98d79 (patch) | |
tree | 189b50c87d0631b6a17b767b12f8972c43b09ec8 /pow.c | |
parent | 352cf452bda69e3edbb9b0dfeb6c9767e9095b0c (diff) | |
download | mpfr-9452541d31eb789a64f9cb963200eddf57c98d79.tar.gz |
implemented exact rounding (but no ternary flag)
git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@1434 280ebfd0-de03-0410-8827-d642c229c3f4
Diffstat (limited to 'pow.c')
-rw-r--r-- | pow.c | 152 |
1 files changed, 104 insertions, 48 deletions
@@ -25,7 +25,7 @@ MA 02111-1307, USA. */ #include "mpfr.h" #include "mpfr-impl.h" -/* sets x to y^n, and returns ceil(log2(max ulp error)) */ +/* sets x to y^n, and return 0 if exact, non-zero otherwise */ int #if __STDC__ mpfr_pow_ui (mpfr_ptr x, mpfr_srcptr y, unsigned long int n, mp_rnd_t rnd) @@ -37,56 +37,82 @@ mpfr_pow_ui (x, y, n, rnd) mp_rnd_t rnd; #endif { - long int i; + long int i, err; unsigned long m; - double err; - mpfr_t ycopy; + mpfr_t res; + mp_prec_t prec; + int inexact; + mp_rnd_t rnd1; - if (MPFR_IS_NAN(y)) { MPFR_SET_NAN(x); return 0; } + if (MPFR_IS_NAN(y)) + { + MPFR_SET_NAN(x); + MPFR_RET_NAN; + } MPFR_CLEAR_NAN(x); + if (n == 0) /* x^0 = 1 for any x */ + { + mpfr_set_ui (x, 1, rnd); + return 0; + } + if (MPFR_IS_INF(y)) { - if (n == 0) { MPFR_SET_NAN(x); return 0; } /* Inf^n = Inf, (-Inf)^n = Inf for n even, -Inf for n odd */ - else if ((MPFR_SIGN(y) < 0) && (n % 2 == 1)) { - if (MPFR_SIGN(x) > 0) MPFR_CHANGE_SIGN(x); - } + if ((MPFR_SIGN(y) < 0) && (n % 2 == 1)) + { + if (MPFR_SIGN(x) > 0) + MPFR_CHANGE_SIGN(x); + } else if (MPFR_SIGN(x) < 0) MPFR_CHANGE_SIGN(x); - MPFR_SET_INF(x); return 0; + MPFR_SET_INF(x); + return 0; } MPFR_CLEAR_INF(x); - if (n==0) { mpfr_set_ui(x, 1, rnd); return 0; } - - if (x == y) { - mpfr_init2 (ycopy, MPFR_PREC(y)); - mpfr_set (ycopy, y, GMP_RNDN); /* exact */ - } - - mpfr_set(x, y, rnd); - err = 1.0; - for (m=n, i=0; m; i++, m>>=1); - /* now 2^(i-1) <= n < 2^i */ - for (i-=2; i>=0; i--) { - mpfr_mul(x, x, x, rnd); err = 2.0*err+1.0; - if (n & (1<<i)) { - mpfr_mul(x, x, (x == y) ? ycopy : y, rnd); - err += 1.0; + mpfr_init (res); + + prec = MPFR_PREC(x); + + rnd1 = (MPFR_SIGN(y) > 0) ? GMP_RNDU : GMP_RNDD; /* away */ + + do + { + prec += 3; + for (m=n, i=0; m; i++, m>>=1, prec++); + mpfr_set_prec (res, prec); + inexact = mpfr_set (res, y, rnd1); + err = 1; + /* now 2^(i-1) <= n < 2^i */ + for (i-=2; i>=0; i--) + { + if (mpfr_mul (res, res, res, GMP_RNDU)) + inexact = 1; + err++; + if (n & (1 << i)) + if (mpfr_mul (res, res, y, rnd1)) + inexact = 1; + } + err = prec - err; + if (err < 0) + err = 0; } - } - - if (x == y) mpfr_clear (ycopy); - return _mpfr_ceil_log2 (err); -} + while (inexact && (mpfr_can_round (res, err, + (MPFR_SIGN(res) > 0) ? GMP_RNDU : GMP_RNDD, rnd, MPFR_PREC(x)) == 0)); -/* sets x to y^n, and returns ceil(log2(max ulp error)) */ + if (mpfr_set (x, res, rnd)) + inexact = 1; -/* TODO: gerer l'infini en cas de debordement ? Fait mecaniquement si fait - dans mul ? */ + mpfr_clear (res); + + return inexact; +} + +/* sets x to y^n, and return 0 if exact, non-zero otherwise */ int #if __STDC__ @@ -100,22 +126,52 @@ mpfr_ui_pow_ui (x, y, n, rnd) mp_rnd_t rnd; #endif { - long int i; double err; + long int i, err; + unsigned long m; + mpfr_t res; + mp_prec_t prec; + int inexact; - MPFR_CLEAR_FLAGS(x); - if (n==0) { mpfr_set_ui(x, 1, rnd); return 0; } + MPFR_CLEAR_FLAGS(x); + + if (n == 0) /* x^0 = 1 for any x */ + { + mpfr_set_ui (x, 1, rnd); + return 0; + } - mpfr_set_ui(x, y, rnd); - err = 1.0; + mpfr_init (res); - for (i=0;(1<<i)<=n;i++); - /* now 2^(i-1) <= n < 2^i */ - for (i-=2; i>=0; i--) { - mpfr_mul(x, x, x, rnd); err = 2.0 * err + 1.0; - if (n & (1<<i)) { - mpfr_mul_ui(x, x, y, rnd); - err = err + 1.0; + prec = MPFR_PREC(x); + + do + { + prec += 3; + for (m=n, i=0; m; i++, m>>=1, prec++); + mpfr_set_prec (res, prec); + inexact = mpfr_set_ui (res, y, GMP_RNDU); + err = 1; + /* now 2^(i-1) <= n < 2^i */ + for (i-=2; i>=0; i--) + { + if (mpfr_mul (res, res, res, GMP_RNDU)) + inexact = 1; + err++; + if (n & (1 << i)) + if (mpfr_mul_ui (res, res, y, GMP_RNDU)) + inexact = 1; + } + err = prec - err; + if (err < 0) + err = 0; } - } - return _mpfr_ceil_log2 (err); + while (inexact && (mpfr_can_round (res, err, + (MPFR_SIGN(res) > 0) ? GMP_RNDU : GMP_RNDD, rnd, MPFR_PREC(x)) == 0)); + + if (mpfr_set (x, res, rnd)) + inexact = 1; + + mpfr_clear (res); + + return inexact; } |