From 4102e576454b24ff22aef95821389ba075d85be1 Mon Sep 17 00:00:00 2001 From: Paul Zimmermann Date: Thu, 13 Feb 2020 09:30:20 +0100 Subject: use extended exponent range in mpc_exp --- src/exp.c | 15 ++++++++++++--- 1 file changed, 12 insertions(+), 3 deletions(-) diff --git a/src/exp.c b/src/exp.c index b28ccea..7770a5a 100644 --- a/src/exp.c +++ b/src/exp.c @@ -28,6 +28,7 @@ mpc_exp (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) int ok = 0; int inex_re, inex_im; int saved_underflow, saved_overflow; + mpfr_exp_t saved_emin, saved_emax; /* special values */ if (mpfr_nan_p (mpc_realref (op)) || mpfr_nan_p (mpc_imagref (op))) @@ -58,7 +59,6 @@ mpc_exp (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) return MPC_INEX(0, 0); /* NaN is exact */ } - if (mpfr_zero_p (mpc_imagref(op))) /* special case when the input is real exp(x-i*0) = exp(x) -i*0, even if x is NaN @@ -77,7 +77,6 @@ mpc_exp (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) return MPC_INEX(inex_re, inex_im); } - if (mpfr_inf_p (mpc_realref (op))) /* real part is an infinity, exp(-inf +i*y) = 0*(cos y +i*sin y) @@ -130,6 +129,10 @@ mpc_exp (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) return MPC_INEX(0, 0); /* NaN is exact */ } + saved_emin = mpfr_get_emin (); + saved_emax = mpfr_get_emax (); + mpfr_set_emin (mpfr_get_emin_min ()); + mpfr_set_emax (mpfr_get_emax_max ()); /* from now on, both parts of op are regular numbers */ @@ -150,7 +153,7 @@ mpc_exp (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) do { - prec += mpc_ceil_log2 (prec) + 5; + prec += prec / 2 + mpc_ceil_log2 (prec) + 5; mpfr_set_prec (x, prec); mpfr_set_prec (y, prec); @@ -199,5 +202,11 @@ mpc_exp (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) if (saved_overflow) mpfr_set_overflow (); + /* restore the exponent range, and check the range of results */ + mpfr_set_emin (saved_emin); + mpfr_set_emax (saved_emax); + inex_re = mpfr_check_range (mpc_realref (rop), inex_re, MPC_RND_RE (rnd)); + inex_im = mpfr_check_range (mpc_imagref (rop), inex_im, MPC_RND_IM (rnd)); + return MPC_INEX(inex_re, inex_im); } -- cgit v1.2.1