diff options
Diffstat (limited to 'src/div.c')
| -rw-r--r-- | src/div.c | 136 |
1 files changed, 72 insertions, 64 deletions
@@ -68,28 +68,28 @@ mpc_div_inf_fin (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w) if (a == 1) if (b == 1) { - mpfr_add (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN); + mpfr_add (sign, mpc_realref (w), mpc_imagref (w), MPFR_RNDN); x = MPC_MPFR_SIGN (sign); - mpfr_sub (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN); + mpfr_sub (sign, mpc_realref (w), mpc_imagref (w), MPFR_RNDN); y = MPC_MPFR_SIGN (sign); } else { /* b == -1 */ - mpfr_sub (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN); + mpfr_sub (sign, mpc_realref (w), mpc_imagref (w), MPFR_RNDN); x = MPC_MPFR_SIGN (sign); - mpfr_add (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN); + mpfr_add (sign, mpc_realref (w), mpc_imagref (w), MPFR_RNDN); y = -MPC_MPFR_SIGN (sign); } else /* a == -1 */ if (b == 1) { - mpfr_sub (sign, mpc_imagref (w), mpc_realref (w), GMP_RNDN); + mpfr_sub (sign, mpc_imagref (w), mpc_realref (w), MPFR_RNDN); x = MPC_MPFR_SIGN (sign); - mpfr_add (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN); + mpfr_add (sign, mpc_realref (w), mpc_imagref (w), MPFR_RNDN); y = MPC_MPFR_SIGN (sign); } else { /* b == -1 */ - mpfr_add (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN); + mpfr_add (sign, mpc_realref (w), mpc_imagref (w), MPFR_RNDN); x = -MPC_MPFR_SIGN (sign); - mpfr_sub (sign, mpc_imagref (w), mpc_realref (w), GMP_RNDN); + mpfr_sub (sign, mpc_imagref (w), mpc_realref (w), MPFR_RNDN); y = MPC_MPFR_SIGN (sign); } mpfr_clear (sign); @@ -121,25 +121,25 @@ mpc_div_fin_inf (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w) mpfr_init2 (x, 2); mpfr_init2 (y, 2); mpfr_init2 (zero, 2); - mpfr_set_ui (zero, 0ul, GMP_RNDN); + mpfr_set_ui (zero, 0ul, MPFR_RNDN); mpfr_init2 (a, mpfr_get_prec (mpc_realref (z))); mpfr_init2 (b, mpfr_get_prec (mpc_imagref (z))); - mpfr_set_ui (c, (mpfr_inf_p (mpc_realref (w)) ? 1 : 0), GMP_RNDN); - MPFR_COPYSIGN (c, c, mpc_realref (w), GMP_RNDN); - mpfr_set_ui (d, (mpfr_inf_p (mpc_imagref (w)) ? 1 : 0), GMP_RNDN); - MPFR_COPYSIGN (d, d, mpc_imagref (w), GMP_RNDN); + mpfr_set_ui (c, (mpfr_inf_p (mpc_realref (w)) ? 1 : 0), MPFR_RNDN); + MPFR_COPYSIGN (c, c, mpc_realref (w), MPFR_RNDN); + mpfr_set_ui (d, (mpfr_inf_p (mpc_imagref (w)) ? 1 : 0), MPFR_RNDN); + MPFR_COPYSIGN (d, d, mpc_imagref (w), MPFR_RNDN); - mpfr_mul (a, mpc_realref (z), c, GMP_RNDN); /* exact */ - mpfr_mul (b, mpc_imagref (z), d, GMP_RNDN); - mpfr_add (x, a, b, GMP_RNDN); + mpfr_mul (a, mpc_realref (z), c, MPFR_RNDN); /* exact */ + mpfr_mul (b, mpc_imagref (z), d, MPFR_RNDN); + mpfr_add (x, a, b, MPFR_RNDN); - mpfr_mul (b, mpc_imagref (z), c, GMP_RNDN); - mpfr_mul (a, mpc_realref (z), d, GMP_RNDN); - mpfr_sub (y, b, a, GMP_RNDN); + mpfr_mul (b, mpc_imagref (z), c, MPFR_RNDN); + mpfr_mul (a, mpc_realref (z), d, MPFR_RNDN); + mpfr_sub (y, b, a, MPFR_RNDN); - MPFR_COPYSIGN (mpc_realref (rop), zero, x, GMP_RNDN); - MPFR_COPYSIGN (mpc_imagref (rop), zero, y, GMP_RNDN); + MPFR_COPYSIGN (mpc_realref (rop), zero, x, MPFR_RNDN); + MPFR_COPYSIGN (mpc_imagref (rop), zero, y, MPFR_RNDN); mpfr_clear (c); mpfr_clear (d); @@ -173,10 +173,10 @@ mpc_div_real (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w, mpc_rnd_t rnd) inexact flags */ if (mpfr_zero_p (mpc_realref (rop))) mpfr_setsign (mpc_realref (rop), mpc_realref (rop), (zrs != wrs && zis != wis), - GMP_RNDN); /* exact */ + MPFR_RNDN); /* exact */ if (mpfr_zero_p (mpc_imagref (rop))) mpfr_setsign (mpc_imagref (rop), mpc_imagref (rop), (zis != wrs && zrs == wis), - GMP_RNDN); + MPFR_RNDN); return MPC_INEX(inex_re, inex_im); } @@ -204,7 +204,7 @@ mpc_div_imag (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w, mpc_rnd_t rnd) wloc[0] = mpc_imagref(w)[0]; /* copies mpfr struct IM(w) into wloc */ inex_re = mpfr_div (mpc_realref(dest), mpc_imagref(z), wloc, MPC_RND_RE(rnd)); - mpfr_neg (wloc, wloc, GMP_RNDN); + mpfr_neg (wloc, wloc, MPFR_RNDN); /* changes the sign only in wloc, not in w; no need to correct later */ inex_im = mpfr_div (mpc_imagref(dest), mpc_realref(z), wloc, MPC_RND_IM(rnd)); @@ -221,10 +221,10 @@ mpc_div_imag (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w, mpc_rnd_t rnd) inexact flags */ if (mpfr_zero_p (mpc_realref (rop))) mpfr_setsign (mpc_realref (rop), mpc_realref (rop), (zrs != wrs && zis != wis), - GMP_RNDN); /* exact */ + MPFR_RNDN); /* exact */ if (imag_z) mpfr_setsign (mpc_imagref (rop), mpc_imagref (rop), (zis != wrs && zrs == wis), - GMP_RNDN); + MPFR_RNDN); return MPC_INEX(inex_re, inex_im); } @@ -292,11 +292,11 @@ mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd) /* first compute norm(c) */ mpfr_clear_underflow (); mpfr_clear_overflow (); - inexact_norm = mpc_norm (q, c, GMP_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, GMP_RNDN); + mpfr_set_ui (q, 0ul, MPFR_RNDN); /* to obtain divisions by 0 later on */ /* now compute b*conjugate(c) */ @@ -312,69 +312,77 @@ mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd) hopefully, the side-effects of mpc_mul do indeed raise the mpfr exceptions */ if (overflow_prod) { + int isinf = 0; tmpsgn = mpfr_sgn (mpc_realref(res)); if (tmpsgn > 0) - mpfr_nextabove (mpc_realref(res)); + { + mpfr_nextabove (mpc_realref(res)); + isinf = mpfr_inf_p (mpc_realref(res)); + mpfr_nextbelow (mpc_realref(res)); + } else if (tmpsgn < 0) - mpfr_nextbelow (mpc_realref(res)); - if (mpfr_inf_p (mpc_realref(res))) + { + mpfr_nextbelow (mpc_realref(res)); + isinf = mpfr_inf_p (mpc_realref(res)); + mpfr_nextabove (mpc_realref(res)); + } + if (isinf) { mpfr_set_inf (mpc_realref(res), tmpsgn); overflow_re = 1; } - else - if (tmpsgn > 0) - mpfr_nextbelow (mpc_realref(res)); - else if (tmpsgn < 0) - mpfr_nextabove (mpc_realref(res)); tmpsgn = mpfr_sgn (mpc_imagref(res)); + isinf = 0; if (tmpsgn > 0) - mpfr_nextabove (mpc_imagref(res)); + { + mpfr_nextabove (mpc_imagref(res)); + isinf = mpfr_inf_p (mpc_imagref(res)); + mpfr_nextbelow (mpc_imagref(res)); + } else if (tmpsgn < 0) - mpfr_nextbelow (mpc_imagref(res)); - if (mpfr_inf_p (mpc_imagref(res))) + { + mpfr_nextbelow (mpc_imagref(res)); + isinf = mpfr_inf_p (mpc_imagref(res)); + mpfr_nextabove (mpc_imagref(res)); + } + if (isinf) { mpfr_set_inf (mpc_imagref(res), tmpsgn); overflow_im = 1; } - else - if (tmpsgn > 0) - mpfr_nextbelow (mpc_imagref(res)); - else if (tmpsgn < 0) - mpfr_nextabove (mpc_imagref(res)); mpc_set (a, res, rnd); goto end; } /* divide the product by the norm */ if (inexact_norm == 0 && (inexact_re == 0 || inexact_im == 0)) { - /* The division has good chances to be exact in at least one part. */ - /* Since this can cause problems when not rounding to the nearest, */ - /* we use the division code of mpfr, which handles the situation. */ + /* The division has good chances to be exact in at least one part. */ + /* Since this can cause problems when not rounding to the nearest, */ + /* we use the division code of mpfr, which handles the situation. */ mpfr_clear_underflow (); mpfr_clear_overflow (); - inexact_re |= mpfr_div (mpc_realref (res), mpc_realref (res), q, GMP_RNDZ); + inexact_re |= mpfr_div (mpc_realref (res), mpc_realref (res), q, MPFR_RNDZ); underflow_re = mpfr_underflow_p (); overflow_re = mpfr_overflow_p (); ok_re = !inexact_re || underflow_re || overflow_re - || mpfr_can_round (mpc_realref (res), prec - 4, GMP_RNDN, - GMP_RNDZ, MPC_PREC_RE(a) + (rnd_re == GMP_RNDN)); + || mpfr_can_round (mpc_realref (res), prec - 4, MPFR_RNDN, + MPFR_RNDZ, MPC_PREC_RE(a) + (rnd_re == MPFR_RNDN)); if (ok_re) /* compute imaginary part */ { mpfr_clear_underflow (); mpfr_clear_overflow (); - inexact_im |= mpfr_div (mpc_imagref (res), mpc_imagref (res), q, GMP_RNDZ); + inexact_im |= mpfr_div (mpc_imagref (res), mpc_imagref (res), q, MPFR_RNDZ); underflow_im = mpfr_underflow_p (); overflow_im = mpfr_overflow_p (); ok_im = !inexact_im || underflow_im || overflow_im - || mpfr_can_round (mpc_imagref (res), prec - 4, GMP_RNDN, - GMP_RNDZ, MPC_PREC_IM(a) + (rnd_im == GMP_RNDN)); + || mpfr_can_round (mpc_imagref (res), prec - 4, MPFR_RNDN, + MPFR_RNDZ, MPC_PREC_IM(a) + (rnd_im == MPFR_RNDN)); } } else { /* The division is inexact, so for efficiency reasons we invert q */ /* only once and multiply by the inverse. */ - if (mpfr_ui_div (q, 1ul, q, GMP_RNDZ) || inexact_norm) { + if (mpfr_ui_div (q, 1ul, q, MPFR_RNDZ) || inexact_norm) { /* if 1/q is inexact, the approximations of the real and imaginary part below will be inexact, unless RE(res) or IM(res) is zero */ @@ -383,22 +391,22 @@ mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd) } mpfr_clear_underflow (); mpfr_clear_overflow (); - inexact_re |= mpfr_mul (mpc_realref (res), mpc_realref (res), q, GMP_RNDZ); + inexact_re |= mpfr_mul (mpc_realref (res), mpc_realref (res), q, MPFR_RNDZ); underflow_re = mpfr_underflow_p (); overflow_re = mpfr_overflow_p (); ok_re = !inexact_re || underflow_re || overflow_re - || mpfr_can_round (mpc_realref (res), prec - 4, GMP_RNDN, - GMP_RNDZ, MPC_PREC_RE(a) + (rnd_re == GMP_RNDN)); + || mpfr_can_round (mpc_realref (res), prec - 4, MPFR_RNDN, + MPFR_RNDZ, MPC_PREC_RE(a) + (rnd_re == MPFR_RNDN)); if (ok_re) /* compute imaginary part */ { mpfr_clear_underflow (); mpfr_clear_overflow (); - inexact_im |= mpfr_mul (mpc_imagref (res), mpc_imagref (res), q, GMP_RNDZ); + inexact_im |= mpfr_mul (mpc_imagref (res), mpc_imagref (res), q, MPFR_RNDZ); underflow_im = mpfr_underflow_p (); overflow_im = mpfr_overflow_p (); ok_im = !inexact_im || underflow_im || overflow_im - || mpfr_can_round (mpc_imagref (res), prec - 4, GMP_RNDN, - GMP_RNDZ, MPC_PREC_IM(a) + (rnd_im == GMP_RNDN)); + || mpfr_can_round (mpc_imagref (res), prec - 4, MPFR_RNDN, + MPFR_RNDZ, MPC_PREC_IM(a) + (rnd_im == MPFR_RNDN)); } } } while ((!ok_re || !ok_im) && !underflow_norm && !overflow_norm @@ -416,16 +424,16 @@ mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd) inexact_re = mpfr_sgn (mpc_realref (res)); } else if (underflow_re || (overflow_norm && !overflow_prod)) { - mpfr_set_zero (mpc_realref (a), mpfr_sgn (mpc_realref (res))); - inexact_re = -mpfr_sgn (mpc_realref (res)); + inexact_re = mpfr_signbit (mpc_realref (res)) ? 1 : -1; + mpfr_set_zero (mpc_realref (a), -inexact_re); } if (overflow_im || (underflow_norm && !underflow_prod)) { mpfr_set_inf (mpc_imagref (a), mpfr_sgn (mpc_imagref (res))); inexact_im = mpfr_sgn (mpc_imagref (res)); } else if (underflow_im || (overflow_norm && !overflow_prod)) { - mpfr_set_zero (mpc_imagref (a), mpfr_sgn (mpc_imagref (res))); - inexact_im = -mpfr_sgn (mpc_imagref (res)); + inexact_im = mpfr_signbit (mpc_imagref (res)) ? 1 : -1; + mpfr_set_zero (mpc_imagref (a), -inexact_im); } mpc_clear (res); |