From aea5ef33f384277e827a0fa8bf530d6378ba08b3 Mon Sep 17 00:00:00 2001 From: enge Date: Mon, 1 Oct 2012 12:15:45 +0000 Subject: log10.c: shortened code without changing functionality git-svn-id: svn://scm.gforge.inria.fr/svn/mpc/trunk@1278 211d60ee-9f03-0410-a15a-8952a2c7a4e4 --- src/log10.c | 127 ++++++++++++++++++------------------------------------------ 1 file changed, 37 insertions(+), 90 deletions(-) diff --git a/src/log10.c b/src/log10.c index a4ff9af..0d35c76 100644 --- a/src/log10.c +++ b/src/log10.c @@ -21,19 +21,17 @@ along with this program. If not, see http://www.gnu.org/licenses/ . #include /* for CHAR_BIT */ #include "mpc-impl.h" -/* Auxiliary functions which implement Ziv's strategy for special cases. - if flag = 0: compute only real part - if flag = 1: compute only imaginary +/* Auxiliary function which implements Ziv's strategy for special cases. Exact cases should be dealt with separately. */ static int -mpc_log10_aux (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd, int flag, int nb) +mpc_log10_aux (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd, int nb) { mp_prec_t prec = (MPFR_PREC_MIN > 4) ? MPFR_PREC_MIN : 4; mpc_t tmp; mpfr_t log10; int ok = 0, ret; - prec = mpfr_get_prec ((flag == 0) ? mpc_realref (rop) : mpc_imagref (rop)); + prec = mpfr_get_prec (mpc_imagref (rop)); prec += 10; mpc_init2 (tmp, prec); mpfr_init2 (log10, prec); @@ -44,40 +42,17 @@ mpc_log10_aux (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd, int flag, int nb) /* In each case we have two roundings, thus the final value is x * (1+u)^2 where x is the exact value, and |u| <= 2^(-prec-1). Thus the error is always less than 3 ulps. */ - switch (nb) - { - case 0: /* imag <- atan2(y/x) */ - mpfr_atan2 (mpc_imagref (tmp), mpc_imagref (op), mpc_realref (op), - MPC_RND_IM (rnd)); - mpfr_div (mpc_imagref (tmp), mpc_imagref (tmp), log10, MPFR_RNDN); - ok = mpfr_can_round (mpc_imagref (tmp), prec - 2, MPFR_RNDN, - MPFR_RNDZ, MPC_PREC_IM(rop) + - (MPC_RND_IM (rnd) == MPFR_RNDN)); - if (ok) - ret = mpfr_set (mpc_imagref (rop), mpc_imagref (tmp), - MPC_RND_IM (rnd)); - break; - case 1: /* real <- log(x) */ - mpfr_log (mpc_realref (tmp), mpc_realref (op), MPC_RND_RE (rnd)); - mpfr_div (mpc_realref (tmp), mpc_realref (tmp), log10, MPFR_RNDN); - ok = mpfr_can_round (mpc_realref (tmp), prec - 2, MPFR_RNDN, - MPFR_RNDZ, MPC_PREC_RE(rop) + - (MPC_RND_RE (rnd) == MPFR_RNDN)); - if (ok) - ret = mpfr_set (mpc_realref (rop), mpc_realref (tmp), - MPC_RND_RE (rnd)); - break; - case 2: /* imag <- pi */ - mpfr_const_pi (mpc_imagref (tmp), MPC_RND_IM (rnd)); - mpfr_div (mpc_imagref (tmp), mpc_imagref (tmp), log10, MPFR_RNDN); - ok = mpfr_can_round (mpc_imagref (tmp), prec - 2, MPFR_RNDN, - MPFR_RNDZ, MPC_PREC_IM(rop) + - (MPC_RND_IM (rnd) == MPFR_RNDN)); - if (ok) - ret = mpfr_set (mpc_imagref (rop), mpc_imagref (tmp), - MPC_RND_IM (rnd)); - break; - } + if (nb == 0) + mpfr_atan2 (mpc_imagref (tmp), mpc_imagref (op), mpc_realref (op), + MPC_RND_IM (rnd)); + else + mpfr_const_pi (mpc_imagref (tmp), MPC_RND_IM (rnd)); + mpfr_div (mpc_imagref (tmp), mpc_imagref (tmp), log10, MPFR_RNDN); + ok = mpfr_can_round (mpc_imagref (tmp), prec - 2, MPFR_RNDN, + MPFR_RNDZ, MPC_PREC_IM(rop) + (MPC_RND_IM (rnd) == MPFR_RNDN)); + if (ok) + ret = mpfr_set (mpc_imagref (rop), mpc_imagref (tmp), + MPC_RND_IM (rnd)); prec += prec / 2; mpc_set_prec (tmp, prec); mpfr_set_prec (log10, prec); @@ -97,46 +72,22 @@ mpc_log10 (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) mpc_t ww; mpc_rnd_t invrnd; - /* special values: NaN and infinities: same as mpc_log */ - if (!mpc_fin_p (op)) /* real or imaginary parts are NaN or Inf */ - { - if (mpfr_nan_p (mpc_realref (op))) - { - if (mpfr_inf_p (mpc_imagref (op))) - /* (NaN, Inf) -> (+Inf, NaN) */ - mpfr_set_inf (mpc_realref (rop), +1); - else - /* (NaN, xxx) -> (NaN, NaN) */ - mpfr_set_nan (mpc_realref (rop)); - mpfr_set_nan (mpc_imagref (rop)); - inex_im = 0; /* Inf/NaN is exact */ - } - else if (mpfr_nan_p (mpc_imagref (op))) - { - if (mpfr_inf_p (mpc_realref (op))) - /* (Inf, NaN) -> (+Inf, NaN) */ - mpfr_set_inf (mpc_realref (rop), +1); - else - /* (xxx, NaN) -> (NaN, NaN) */ - mpfr_set_nan (mpc_realref (rop)); - mpfr_set_nan (mpc_imagref (rop)); - inex_im = 0; /* Inf/NaN is exact */ - } - else /* We have an infinity in at least one part. */ - { - /* (+Inf, y) -> (+Inf, 0) for finite positive-signed y */ - if (mpfr_inf_p (mpc_realref (op)) && mpfr_signbit (mpc_realref (op)) - == 0 && mpfr_number_p (mpc_imagref (op))) - inex_im = mpfr_atan2 (mpc_imagref (rop), mpc_imagref (op), - mpc_realref (op), MPC_RND_IM (rnd)); - else - /* (xxx, Inf) -> (+Inf, atan2(Inf/xxx)) - (Inf, yyy) -> (+Inf, atan2(yyy/Inf)) */ - inex_im = mpc_log10_aux (rop, op, rnd, 1, 0); - mpfr_set_inf (mpc_realref (rop), +1); - } - return MPC_INEX(0, inex_im); - } + /* special values: NaN and infinities: same as mpc_log */ + if (!mpc_fin_p (op)) { + if ( mpfr_nan_p (mpc_realref (op)) + || mpfr_nan_p (mpc_imagref (op)) + || ( mpfr_inf_p (mpc_realref (op)) + && mpfr_signbit (mpc_realref (op))== 0 + && mpfr_number_p (mpc_imagref (op)))) + return mpc_log (rop, op, rnd); + else { + /* (xxx, Inf) -> (+Inf, atan2(Inf/xxx)) + (Inf, yyy) -> (+Inf, atan2(yyy/Inf)) */ + inex_im = mpc_log10_aux (rop, op, rnd, 0); + mpfr_set_inf (mpc_realref (rop), +1); + return MPC_INEX(0, inex_im); + } + } /* special cases: real and purely imaginary numbers */ re_cmp = mpfr_cmp_ui (mpc_realref (op), 0); @@ -149,7 +100,7 @@ mpc_log10 (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) inex_im = mpfr_atan2 (mpc_imagref (rop), mpc_imagref (op), mpc_realref (op), MPC_RND_IM (rnd)); else - inex_im = mpc_log10_aux (rop, op, rnd, 1, 0); + inex_im = mpc_log10_aux (rop, op, rnd, 0); mpfr_set_inf (mpc_realref (rop), -1); inex_re = 0; /* -Inf is exact */ } @@ -167,14 +118,10 @@ mpc_log10 (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) rnd_im = INV_RND (MPC_RND_IM (rnd)); else rnd_im = MPC_RND_IM (rnd); - ww->re[0] = *mpc_realref (op); - MPFR_CHANGE_SIGN (ww->re); - ww->im[0] = *mpc_imagref (op); - if (mpfr_cmp_ui (ww->re, 1) == 0) - inex_re = mpfr_set_ui (mpc_realref (rop), 0, MPC_RND_RE (rnd)); - else - inex_re = mpc_log10_aux (rop, ww, rnd, 0, 1); - inex_im = mpc_log10_aux (rop, op, MPC_RND (0,rnd_im), 1, 2); + w [0] = *mpc_realref (op); + MPFR_CHANGE_SIGN (w); + inex_re = mpfr_log10 (mpc_realref (rop), w, MPC_RND_RE (rnd)); + inex_im = mpc_log10_aux (rop, op, MPC_RND (0,rnd_im), 2); if (negative_zero) { mpc_conj (rop, rop, MPC_RNDNN); @@ -188,7 +135,7 @@ mpc_log10 (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) if (im_cmp > 0) { inex_re = mpfr_log10 (mpc_realref (rop), mpc_imagref (op), MPC_RND_RE (rnd)); - inex_im = mpc_log10_aux (rop, op, rnd, 1, 2); + inex_im = mpc_log10_aux (rop, op, rnd, 2); /* division by 2 does not change the ternary flag */ mpfr_div_2ui (mpc_imagref (rop), mpc_imagref (rop), 1, MPFR_RNDN); } @@ -198,7 +145,7 @@ mpc_log10 (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) MPFR_CHANGE_SIGN (w); inex_re = mpfr_log10 (mpc_realref (rop), w, MPC_RND_RE (rnd)); invrnd = MPC_RND (0, INV_RND (MPC_RND_IM (rnd))); - inex_im = mpc_log10_aux (rop, op, invrnd, 1, 2); + inex_im = mpc_log10_aux (rop, op, invrnd, 2); /* division by 2 does not change the ternary flag */ mpfr_div_2ui (mpc_imagref (rop), mpc_imagref (rop), 1, MPFR_RNDN); mpfr_neg (mpc_imagref (rop), mpc_imagref (rop), MPFR_RNDN); -- cgit v1.2.1