diff options
author | zimmerma <zimmerma@280ebfd0-de03-0410-8827-d642c229c3f4> | 1999-12-08 12:12:07 +0000 |
---|---|---|
committer | zimmerma <zimmerma@280ebfd0-de03-0410-8827-d642c229c3f4> | 1999-12-08 12:12:07 +0000 |
commit | 3778f7a4b9ddc6c2685cffbcd348f83818601e1c (patch) | |
tree | 5450e636e8a4fee190c002244bae44e3632913f7 /exp.c | |
parent | 049faf5e4b4e2f7ac9f3f2f85c90b5fcbdeed885 (diff) | |
download | mpfr-3778f7a4b9ddc6c2685cffbcd348f83818601e1c.tar.gz |
fixed bug found by V. Lefe`vre (when n<0, we have to compute
an upper bound of log(2) instead of a lower bound)
git-svn-id: svn://scm.gforge.inria.fr/svn/mpfr/trunk@413 280ebfd0-de03-0410-8827-d642c229c3f4
Diffstat (limited to 'exp.c')
-rw-r--r-- | exp.c | 13 |
1 files changed, 10 insertions, 3 deletions
@@ -71,27 +71,33 @@ mpfr_exp(y, x, rnd_mode) n = (int) floor(mpfr_get_d(x)/LOG2); - /* K = (int) log( (double) precy + 2 ); */ K = (int) sqrt( (double) precy ); l = (precy-1)/K + 1; err = K + (int) ceil(log(2.0*(double)l+18.0)/LOG2); /* add K extra bits, i.e. failure probability <= 1/2^K = O(1/precy) */ q = precy + err + K + 3; mpfr_init2(r, q); mpfr_init2(s, q); mpfr_init2(t, q); + /* the algorithm consists in computing an upper bound of exp(x) using + a precision of q bits, and see if we can round to PREC(y) taking + into account the maximal error. Otherwise we increase q. */ do { #ifdef DEBUG printf("n=%d K=%d l=%d q=%d\n",n,K,l,q); #endif - mpfr_log2(s, GMP_RNDZ); + /* if n<0, we have to get an upper bound of log(2) + in order to get an upper bound of r = x-n*log(2) */ + mpfr_log2(s, (n>=0) ? GMP_RNDZ : GMP_RNDU); #ifdef DEBUG printf("n=%d log(2)=",n); mpfr_print_raw(s); putchar('\n'); #endif - mpfr_mul_ui(r, s, (n<0) ? -n : n, GMP_RNDZ); /* r = n*log(2) */ + mpfr_mul_ui(r, s, (n<0) ? -n : n, (n>=0) ? GMP_RNDZ : GMP_RNDU); if (n<0) mpfr_neg(r, r, GMP_RNDD); + /* r = floor(n*log(2)) */ #ifdef DEBUG printf("x=%1.20e\n",mpfr_get_d(x)); + printf(" ="); mpfr_print_raw(x); putchar('\n'); printf("r=%1.20e\n",mpfr_get_d(r)); printf(" ="); mpfr_print_raw(r); putchar('\n'); #endif @@ -104,6 +110,7 @@ mpfr_exp(y, x, rnd_mode) } #ifdef DEBUG printf("x-r=%1.20e\n",mpfr_get_d(r)); + printf(" ="); mpfr_print_raw(r); putchar('\n'); if (SIGN(r)<0) { fprintf(stderr,"Error in mpfr_exp: r<0\n"); exit(1); } #endif mpfr_div_2exp(r, r, K, GMP_RNDU); /* r = (x-n*log(2))/2^K */ |