diff options
Diffstat (limited to 'src/log.c')
| -rw-r--r-- | src/log.c | 129 |
1 files changed, 99 insertions, 30 deletions
@@ -23,12 +23,16 @@ 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=0; - mpfr_t w; + int ok, underflow = 0; + mpfr_srcptr x, y; + mpfr_t v, w; mpfr_prec_t prec; - int loops = 0; + int loops; int re_cmp, im_cmp; int inex_re, inex_im; + int err; + mpfr_exp_t expw; + int sgnw; /* special values: NaN and infinities */ if (!mpc_fin_p (op)) { @@ -96,7 +100,7 @@ mpc_log (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd){ inex_re = mpfr_log (mpc_realref (rop), mpc_imagref (op), MPC_RND_RE (rnd)); inex_im = mpfr_const_pi (mpc_imagref (rop), MPC_RND_IM (rnd)); /* division by 2 does not change the ternary flag */ - mpfr_div_2ui (mpc_imagref (rop), mpc_imagref (rop), 1, GMP_RNDN); + mpfr_div_2ui (mpc_imagref (rop), mpc_imagref (rop), 1, MPFR_RNDN); } else { w [0] = *mpc_imagref (op); @@ -104,45 +108,110 @@ mpc_log (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd){ inex_re = mpfr_log (mpc_realref (rop), w, MPC_RND_RE (rnd)); inex_im = mpfr_const_pi (mpc_imagref (rop), INV_RND (MPC_RND_IM (rnd))); /* division by 2 does not change the ternary flag */ - mpfr_div_2ui (mpc_imagref (rop), mpc_imagref (rop), 1, GMP_RNDN); - mpfr_neg (mpc_imagref (rop), mpc_imagref (rop), GMP_RNDN); + mpfr_div_2ui (mpc_imagref (rop), mpc_imagref (rop), 1, MPFR_RNDN); + mpfr_neg (mpc_imagref (rop), mpc_imagref (rop), MPFR_RNDN); inex_im = -inex_im; /* negate the ternary flag */ } return MPC_INEX(inex_re, inex_im); } prec = MPC_PREC_RE(rop); - mpfr_init2 (w, prec); - /* let op = x + iy; log = 1/2 log (x^2 + y^2) + i atan2 (y, x) */ - /* loop for the real part: log (x^2 + y^2) */ - do { - loops ++; - prec += (loops <= 2) ? mpc_ceil_log2 (prec) + 4 : prec / 2; + mpfr_init2 (w, 2); + /* let op = x + iy; log = 1/2 log (x^2 + y^2) + i atan2 (y, x) */ + /* loop for the real part: 1/2 log (x^2 + y^2), fast, but unsafe */ + /* implementation */ + ok = 0; + for (loops = 1; !ok && loops <= 2; loops++) { + prec += mpc_ceil_log2 (prec) + 4; mpfr_set_prec (w, prec); - /* w is rounded down */ - mpc_norm (w, op, GMP_RNDD); - /* error 1 ulp */ - MPC_ASSERT (!mpfr_inf_p (w)); - /* FIXME: intermediate overflow; the logarithm may be representable */ - - mpfr_log (w, w, GMP_RNDD); - /* generic error of log: (2^(2 - exp(w)) + 1) ulp */ - - if (mpfr_get_exp (w) >= 2) - ok = mpfr_can_round (w, prec - 2, GMP_RNDD, - MPC_RND_RE(rnd), MPC_PREC_RE(rop)); - else - ok = mpfr_can_round (w, prec - 3 + mpfr_get_exp (w), GMP_RNDD, - MPC_RND_RE(rnd), MPC_PREC_RE(rop)); - } while (ok == 0); + mpc_abs (w, op, MPFR_RNDN); + /* error 0.5 ulp */ + if (mpfr_inf_p (w)) + /* intermediate overflow; the logarithm may be representable. + Intermediate underflow is impossible. */ + break; + + mpfr_log (w, w, MPFR_RNDN); + /* generic error of log: (2^(- exp(w)) + 0.5) ulp */ + + if (mpfr_zero_p (w)) + /* impossible to round, switch to second algorithm */ + break; + + err = MPC_MAX (-mpfr_get_exp (w), 0) + 1; + /* number of lost digits */ + ok = mpfr_can_round (w, prec - err, MPFR_RNDN, MPFR_RNDZ, + mpfr_get_prec (mpc_realref (rop)) + (MPC_RND_RE (rnd) == MPFR_RNDN)); + } + + if (!ok) { + prec = MPC_PREC_RE(rop); + mpfr_init2 (v, 2); + /* compute 1/2 log (x^2 + y^2) = log |x| + 1/2 * log (1 + (y/x)^2) + if |x| >= |y|; otherwise, exchange x and y */ + if (mpfr_cmpabs (mpc_realref (op), mpc_imagref (op)) >= 0) { + x = mpc_realref (op); + y = mpc_imagref (op); + } + else { + 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, MPFR_RNDD); /* error 1 ulp */ + mpfr_sqr (v, v, MPFR_RNDD); + /* generic error of multiplication: + 1 + 2*1*(2+1*2^(1-prec)) <= 5.0625 since prec >= 6 */ + mpfr_log1p (v, v, MPFR_RNDD); + /* error 1 + 4*5.0625 = 21.25 , see algorithms.tex */ + mpfr_div_2ui (v, v, 1, MPFR_RNDD); + /* If the result is 0, then there has been an underflow somewhere. */ + + mpfr_abs (w, x, MPFR_RNDN); /* exact */ + mpfr_log (w, w, MPFR_RNDN); /* error 0.5 ulp */ + expw = mpfr_get_exp (w); + sgnw = mpfr_signbit (w); + + mpfr_add (w, w, v, MPFR_RNDN); + if (!sgnw) /* v is positive, so no cancellation; + error 22.25 ulp; error counts lost bits */ + err = 5; + else + 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; + + /* 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, MPFR_RNDN, MPFR_RNDZ, + mpfr_get_prec (mpc_realref (rop)) + (MPC_RND_RE (rnd) == MPFR_RNDN))); + mpfr_clear (v); + } /* imaginary part */ inex_im = mpfr_atan2 (mpc_imagref (rop), mpc_imagref (op), mpc_realref (op), MPC_RND_IM (rnd)); - /* set the real part; cannot be done before when rop==op */ - inex_re = mpfr_div_2ui (mpc_realref(rop), w, 1ul, MPC_RND_RE (rnd)); + /* set the real part; cannot be done before if rop==op */ + if (underflow) + /* create underflow in result */ + inex_re = mpfr_set_ui_2exp (mpc_realref (rop), 1, + mpfr_get_emin_min () - 2, 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); } |