From fbcde25676931fbdbf707d03fd980c1936a00b02 Mon Sep 17 00:00:00 2001 From: Paul Zimmermann Date: Thu, 6 Feb 2020 13:01:20 +0100 Subject: [div.c] use larger exponent range internally (fixes bug20200206 in tdiv.c) --- src/div.c | 75 ++++++++++++++++----------------------------------------------- 1 file changed, 19 insertions(+), 56 deletions(-) diff --git a/src/div.c b/src/div.c index a84b3c6..0915ee4 100644 --- a/src/div.c +++ b/src/div.c @@ -244,9 +244,7 @@ mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd) mpfr_rnd_t rnd_re = MPC_RND_RE (rnd), rnd_im = MPC_RND_IM (rnd); int saved_underflow, saved_overflow; int tmpsgn; - mpfr_exp_t e, emin, emax, emid; /* for scaling of exponents */ - mpc_t b_scaled, c_scaled; - mpfr_t b_re, b_im, c_re, c_im; + mpfr_exp_t saved_emin, saved_emax; /* According to the C standard G.3, there are three types of numbers: */ /* finite (both parts are usual real numbers; contains 0), infinite */ @@ -278,57 +276,16 @@ mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd) mpc_init2 (res, 2); mpfr_init (q); - /* compute scaling of exponents: none of Re(c) and Im(c) can be zero, - but one of Re(b) or Im(b) could be zero */ - - e = mpfr_get_exp (mpc_realref (c)); - emin = emax = e; - e = mpfr_get_exp (mpc_imagref (c)); - if (e > emax) - emax = e; - else if (e < emin) - emin = e; - if (!mpfr_zero_p (mpc_realref (b))) - { - e = mpfr_get_exp (mpc_realref (b)); - if (e > emax) - emax = e; - else if (e < emin) - emin = e; - } - if (!mpfr_zero_p (mpc_imagref (b))) - { - e = mpfr_get_exp (mpc_imagref (b)); - if (e > emax) - emax = e; - else if (e < emin) - emin = e; - } - - /* all input exponents are in [emin, emax] */ - emid = emin / 2 + emax / 2; - - /* scale the inputs */ - b_re[0] = mpc_realref (b)[0]; - if (!mpfr_zero_p (mpc_realref (b))) - MPFR_EXP(b_re) = MPFR_EXP(mpc_realref (b)) - emid; - b_im[0] = mpc_imagref (b)[0]; - if (!mpfr_zero_p (mpc_imagref (b))) - MPFR_EXP(b_im) = MPFR_EXP(mpc_imagref (b)) - emid; - c_re[0] = mpc_realref (c)[0]; - MPFR_EXP(c_re) = MPFR_EXP(mpc_realref (c)) - emid; - c_im[0] = mpc_imagref (c)[0]; - MPFR_EXP(c_im) = MPFR_EXP(mpc_imagref (c)) - emid; - - /* create the scaled inputs without allocating new memory */ - mpc_realref (b_scaled)[0] = b_re[0]; - mpc_imagref (b_scaled)[0] = b_im[0]; - mpc_realref (c_scaled)[0] = c_re[0]; - mpc_imagref (c_scaled)[0] = c_im[0]; + /* we perform the division in the largest possible exponent range, + to avoid underflow/overflow in intermediate computations */ + saved_emin = mpfr_get_emin (); + saved_emax = mpfr_get_emax (); + mpfr_set_emin (mpfr_get_emin_min ()); + mpfr_set_emax (mpfr_get_emax_max ()); /* create the conjugate of c in c_conj without allocating new memory */ - mpc_realref (c_conj)[0] = mpc_realref (c_scaled)[0]; - mpc_imagref (c_conj)[0] = mpc_imagref (c_scaled)[0]; + mpc_realref (c_conj)[0] = mpc_realref (c)[0]; + mpc_imagref (c_conj)[0] = mpc_imagref (c)[0]; MPFR_CHANGE_SIGN (mpc_imagref (c_conj)); /* save the underflow or overflow flags from MPFR */ @@ -342,20 +299,20 @@ mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd) mpc_set_prec (res, prec); mpfr_set_prec (q, prec); - /* first compute norm(c_scaled) */ + /* first compute norm(c) */ mpfr_clear_underflow (); mpfr_clear_overflow (); - inexact_norm = mpc_norm (q, c_scaled, MPFR_RNDU); + inexact_norm = mpc_norm (q, c, MPFR_RNDU); underflow_norm = mpfr_underflow_p (); overflow_norm = mpfr_overflow_p (); if (underflow_norm) mpfr_set_ui (q, 0ul, MPFR_RNDN); /* to obtain divisions by 0 later on */ - /* now compute b_scaled*conjugate(c_scaled) */ + /* now compute b*conjugate(c) */ mpfr_clear_underflow (); mpfr_clear_overflow (); - inexact_prod = mpc_mul (res, b_scaled, c_conj, MPC_RNDZZ); + inexact_prod = mpc_mul (res, b, c_conj, MPC_RNDZZ); inexact_re = MPC_INEX_RE (inexact_prod); inexact_im = MPC_INEX_IM (inexact_prod); underflow_prod = mpfr_underflow_p (); @@ -505,5 +462,11 @@ mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, 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); + inexact_re = mpfr_check_range (mpc_realref (a), inexact_re, rnd_re); + inexact_im = mpfr_check_range (mpc_imagref (a), inexact_im, rnd_im); + return MPC_INEX (inexact_re, inexact_im); } -- cgit v1.2.1