diff options
| author | enge <enge@211d60ee-9f03-0410-a15a-8952a2c7a4e4> | 2012-06-30 14:32:11 +0000 |
|---|---|---|
| committer | enge <enge@211d60ee-9f03-0410-a15a-8952a2c7a4e4> | 2012-06-30 14:32:11 +0000 |
| commit | 4fa6ba816d900b6df71065ad4971e3ebb74d3056 (patch) | |
| tree | cb026c80e21b2669a0cae48c6abe2579314541e6 | |
| parent | 22ec468b9cec80aa3ffdf58dc853b8cf70f5236c (diff) | |
| download | mpc-4fa6ba816d900b6df71065ad4971e3ebb74d3056.tar.gz | |
log.c: correctly handle log (+-1 + i*eps)
log.dat: add corresponding test cases
git-svn-id: svn://scm.gforge.inria.fr/svn/mpc/trunk@1213 211d60ee-9f03-0410-a15a-8952a2c7a4e4
| -rw-r--r-- | src/log.c | 38 | ||||
| -rw-r--r-- | tests/log.dat | 12 |
2 files changed, 35 insertions, 15 deletions
@@ -23,7 +23,7 @@ along with this program. If not, see http://www.gnu.org/licenses/ . int mpc_log (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd){ - int ok; + int ok, underflow = 0; mpfr_srcptr x, y; mpfr_t v, w; mpfr_prec_t prec; @@ -170,18 +170,20 @@ do { \ x = mpc_imagref (op); y = mpc_realref (op); } + do { prec += mpc_ceil_log2 (prec) + 4; mpfr_set_prec (v, prec); mpfr_set_prec (w, prec); - mpfr_div (v, y, x, GMP_RNDN); /* error 0.5 ulp */ - mpfr_sqr (v, v, GMP_RNDN); + mpfr_div (v, y, x, GMP_RNDD); /* error 1 ulp */ + mpfr_sqr (v, v, GMP_RNDD); /* generic error of multiplication: - 0.5 + 2*0.5*(2+0.5*2^(1-prec)) <= 2.51... since prec >= 6 */ - mpfr_log1p (v, v, GMP_RNDN); - /* error 0.5 + 4*2.51... = 10.54... , see algorithms.tex */ - mpfr_div_2ui (v, v, 1, GMP_RNDN); + 1 + 2*1*(2+1*2^(1-prec)) <= 5.0625 since prec >= 6 */ + mpfr_log1p (v, v, GMP_RNDD); + /* error 1 + 4*5.0625 = 21.25 , see algorithms.tex */ + mpfr_div_2ui (v, v, 1, GMP_RNDD); + /* If the result is 0, then there has been an underflow somewhere. */ mpfr_abs (w, x, GMP_RNDN); /* exact */ mpfr_log (w, w, GMP_RNDN); /* error 0.5 ulp */ @@ -192,13 +194,20 @@ do { \ error 11.54... ulp; error counts lost bits */ err = 4; else - err = MPC_MAX (4 + mpfr_get_exp (v), - /* 10.54... ulp (v) rewritten in ulp (result, now in w) */ + err = MPC_MAX (5 + mpfr_get_exp (v), + /* 21.25 ulp (v) rewritten in ulp (result, now in w) */ -1 + expw - mpfr_get_exp (w) /* 0.5 ulp (previous w), rewritten in ulp (result) */ ) + 2; - } while (!mpfr_can_round (w, prec - err, GMP_RNDN, GMP_RNDZ, + /* handle one special case: |x|=1, and (y/x)^2 underflows; + then 1/2*log(x^2+y^2) \approx 1/2*y^2 also underflows. */ + if ( (mpfr_cmp_si (x, -1) == 0 || mpfr_cmp_ui (x, 1) == 0) + && mpfr_zero_p (w)) + underflow = 1; + + } while (!underflow && + !mpfr_can_round (w, prec - err, GMP_RNDN, GMP_RNDZ, mpfr_get_prec (mpc_realref (rop)) + (MPC_RND_RE (rnd) == GMP_RNDN))); mpfr_clear (v); } @@ -208,7 +217,14 @@ do { \ MPC_RND_IM (rnd)); /* set the real part; cannot be done before if rop==op */ - inex_re = mpfr_set (mpc_realref (rop), w, MPC_RND_RE (rnd)); + if (underflow) { + /* create underflow in result */ + mpfr_set_ui (mpc_realref (rop), 1, GMP_RNDN); + inex_re = mpfr_mul_2si (mpc_realref (rop), mpc_realref (rop), mpfr_get_emin_min (), + MPC_RND_RE (rnd)); + } + else + inex_re = mpfr_set (mpc_realref (rop), w, MPC_RND_RE (rnd)); mpfr_clear (w); return MPC_INEX(inex_re, inex_im); } diff --git a/tests/log.dat b/tests/log.dat index 854ebf0..9b7aab2 100644 --- a/tests/log.dat +++ b/tests/log.dat @@ -173,7 +173,11 @@ # Example leading to intermediate underflow in x^2+y^2 - - 53 -0x58B90BFB3775A8p-25 2 0x3p-2 2 0x1p-1073741824 2 0x1p-1073741824 D D -# Due to intermediate "underflow" (see FIXME in log.c), the following example -# loops. -# log(-1 + i*eps) -# - - 2 0 2 3 2 -1 2 0x1p-1073741813 N N +# log (-1 + i*eps), looped in previous version +- - 2 0 2 3 2 -1 2 0x1p-1073741813 N N +- - 2 0 2 3 2 -1 2 0x1p-1073741813 D D ++ + 2 0x1p-1073741824 2 4 2 -1 2 0x1p-1073741813 U U +# log(1 + i*eps), could cause similar problems +- + 2 0 2 0x1p-1073741824 2 1 2 0x1p-1073741824 N N +- - 2 0 2 0 2 1 2 0x1p-1073741824 D D ++ + 2 0x1p-1073741824 2 0x1p-1073741824 2 1 2 0x1p-1073741824 U U |