/* mpc_div -- Divide two complex numbers.
-Copyright (C) 2002, 2003, 2004, 2005, 2008, 2009, 2010, 2011, 2012 INRIA
+Copyright (C) 2002, 2003, 2004, 2005, 2008, 2009 Andreas Enge, Paul Zimmermann, Philippe Th\'eveny
-This file is part of GNU MPC.
+This file is part of the MPC Library.
-GNU MPC is free software; you can redistribute it and/or modify it under
-the terms of the GNU Lesser General Public License as published by the
-Free Software Foundation; either version 3 of the License, or (at your
+The MPC Library is free software; you can redistribute it and/or modify
+it under the terms of the GNU Lesser General Public License as published by
+the Free Software Foundation; either version 2.1 of the License, or (at your
option) any later version.
-GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY
-WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
-FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
-more details.
+The MPC Library is distributed in the hope that it will be useful, but
+WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
+or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
+License for more details.
You should have received a copy of the GNU Lesser General Public License
-along with this program. If not, see http://www.gnu.org/licenses/ .
-*/
+along with the MPC Library; see the file COPYING.LIB. If not, write to
+the Free Software Foundation, Inc., 59 Temple Place - Suite 330, Boston,
+MA 02111-1307, USA. */
#include "mpc-impl.h"
-/* this routine deals with the case where w is zero */
static int
mpc_div_zero (mpc_ptr a, mpc_srcptr z, mpc_srcptr w, mpc_rnd_t rnd)
-/* Assumes w==0, implementation according to C99 G.5.1.8 */
{
- int sign = MPFR_SIGNBIT (mpc_realref (w));
+ /* Assumes w==0, implementation according to C99 G.5.1.8 */
+ int sign = MPFR_SIGNBIT (MPC_RE (w));
mpfr_t infty;
-
- mpfr_init2 (infty, MPFR_PREC_MIN);
mpfr_set_inf (infty, sign);
- mpfr_mul (mpc_realref (a), infty, mpc_realref (z), MPC_RND_RE (rnd));
- mpfr_mul (mpc_imagref (a), infty, mpc_imagref (z), MPC_RND_IM (rnd));
- mpfr_clear (infty);
+ mpfr_mul (MPC_RE (a), infty, MPC_RE (z), MPC_RND_RE (rnd));
+ mpfr_mul (MPC_IM (a), infty, MPC_IM (z), MPC_RND_IM (rnd));
return MPC_INEX (0, 0); /* exact */
}
-/* this routine deals with the case where z is infinite and w finite */
static int
mpc_div_inf_fin (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w)
-/* Assumes w finite and non-zero and z infinite; implementation
- according to C99 G.5.1.8 */
{
+ /* Assumes w finite and non-zero and z infinite; implementation
+ according to C99 G.5.1.8 */
int a, b, x, y;
- a = (mpfr_inf_p (mpc_realref (z)) ? MPFR_SIGNBIT (mpc_realref (z)) : 0);
- b = (mpfr_inf_p (mpc_imagref (z)) ? MPFR_SIGNBIT (mpc_imagref (z)) : 0);
-
- /* a is -1 if Re(z) = -Inf, 1 if Re(z) = +Inf, 0 if Re(z) is finite
- b is -1 if Im(z) = -Inf, 1 if Im(z) = +Inf, 0 if Im(z) is finite */
+ a = (mpfr_inf_p (MPC_RE (z)) ? MPFR_SIGNBIT (MPC_RE (z)) : 0);
+ b = (mpfr_inf_p (MPC_IM (z)) ? MPFR_SIGNBIT (MPC_IM (z)) : 0);
- /* x = MPC_MPFR_SIGN (a * mpc_realref (w) + b * mpc_imagref (w)) */
- /* y = MPC_MPFR_SIGN (b * mpc_realref (w) - a * mpc_imagref (w)) */
+ /* x = MPC_MPFR_SIGN (a * MPC_RE (w) + b * MPC_IM (w)) */
+ /* y = MPC_MPFR_SIGN (b * MPC_RE (w) - a * MPC_IM (w)) */
if (a == 0 || b == 0) {
- /* only one of a or b can be zero, since z is infinite */
- x = a * MPC_MPFR_SIGN (mpc_realref (w)) + b * MPC_MPFR_SIGN (mpc_imagref (w));
- y = b * MPC_MPFR_SIGN (mpc_realref (w)) - a * MPC_MPFR_SIGN (mpc_imagref (w));
+ x = a * MPC_MPFR_SIGN (MPC_RE (w)) + b * MPC_MPFR_SIGN (MPC_IM (w));
+ y = b * MPC_MPFR_SIGN (MPC_RE (w)) - a * MPC_MPFR_SIGN (MPC_IM (w));
}
else {
/* Both parts of z are infinite; x could be determined by sign
considerations and comparisons. Since operations with non-finite
numbers are not considered time-critical, we let mpfr do the work. */
mpfr_t sign;
-
mpfr_init2 (sign, 2);
- /* This is enough to determine the sign of sums and differences. */
+ /* This is enough to determine the sign of sums and differences. */
if (a == 1)
if (b == 1) {
- mpfr_add (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN);
+ mpfr_add (sign, MPC_RE (w), MPC_IM (w), GMP_RNDN);
x = MPC_MPFR_SIGN (sign);
- mpfr_sub (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN);
+ mpfr_sub (sign, MPC_RE (w), MPC_IM (w), GMP_RNDN);
y = MPC_MPFR_SIGN (sign);
}
else { /* b == -1 */
- mpfr_sub (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN);
+ mpfr_sub (sign, MPC_RE (w), MPC_IM (w), GMP_RNDN);
x = MPC_MPFR_SIGN (sign);
- mpfr_add (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN);
+ mpfr_add (sign, MPC_RE (w), MPC_IM (w), GMP_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_IM (w), MPC_RE (w), GMP_RNDN);
x = MPC_MPFR_SIGN (sign);
- mpfr_add (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN);
+ mpfr_add (sign, MPC_RE (w), MPC_IM (w), GMP_RNDN);
y = MPC_MPFR_SIGN (sign);
}
else { /* b == -1 */
- mpfr_add (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN);
+ mpfr_add (sign, MPC_RE (w), MPC_IM (w), GMP_RNDN);
x = -MPC_MPFR_SIGN (sign);
- mpfr_sub (sign, mpc_imagref (w), mpc_realref (w), GMP_RNDN);
+ mpfr_sub (sign, MPC_IM (w), MPC_RE (w), GMP_RNDN);
y = MPC_MPFR_SIGN (sign);
}
mpfr_clear (sign);
}
if (x == 0)
- mpfr_set_nan (mpc_realref (rop));
+ mpfr_set_nan (MPC_RE (rop));
else
- mpfr_set_inf (mpc_realref (rop), x);
+ mpfr_set_inf (MPC_RE (rop), x);
if (y == 0)
- mpfr_set_nan (mpc_imagref (rop));
+ mpfr_set_nan (MPC_IM (rop));
else
- mpfr_set_inf (mpc_imagref (rop), y);
+ mpfr_set_inf (MPC_IM (rop), y);
return MPC_INEX (0, 0); /* exact */
}
-/* this routine deals with the case where z if finite and w infinite */
static int
mpc_div_fin_inf (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w)
-/* Assumes z finite and w infinite; implementation according to
- C99 G.5.1.8 */
{
+ /* Assumes z finite and w infinite; implementation according to
+ C99 G.5.1.8 */
mpfr_t c, d, a, b, x, y, zero;
mpfr_init2 (c, 2); /* needed to hold a signed zero, +1 or -1 */
mpfr_init2 (y, 2);
mpfr_init2 (zero, 2);
mpfr_set_ui (zero, 0ul, GMP_RNDN);
- mpfr_init2 (a, mpfr_get_prec (mpc_realref (z)));
- mpfr_init2 (b, mpfr_get_prec (mpc_imagref (z)));
+ mpfr_init2 (a, mpfr_get_prec (MPC_RE (z)));
+ mpfr_init2 (b, mpfr_get_prec (MPC_IM (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_RE (w)) ? 1 : 0), GMP_RNDN);
+ MPFR_COPYSIGN (c, c, MPC_RE (w), GMP_RNDN);
+ mpfr_set_ui (d, (mpfr_inf_p (MPC_IM (w)) ? 1 : 0), GMP_RNDN);
+ MPFR_COPYSIGN (d, d, MPC_IM (w), GMP_RNDN);
- mpfr_mul (a, mpc_realref (z), c, GMP_RNDN); /* exact */
- mpfr_mul (b, mpc_imagref (z), d, GMP_RNDN);
+ mpfr_mul (a, MPC_RE (z), c, GMP_RNDN); /* exact */
+ mpfr_mul (b, MPC_IM (z), d, GMP_RNDN);
mpfr_add (x, a, b, GMP_RNDN);
- mpfr_mul (b, mpc_imagref (z), c, GMP_RNDN);
- mpfr_mul (a, mpc_realref (z), d, GMP_RNDN);
+ mpfr_mul (b, MPC_IM (z), c, GMP_RNDN);
+ mpfr_mul (a, MPC_RE (z), d, GMP_RNDN);
mpfr_sub (y, b, a, GMP_RNDN);
- MPFR_COPYSIGN (mpc_realref (rop), zero, x, GMP_RNDN);
- MPFR_COPYSIGN (mpc_imagref (rop), zero, y, GMP_RNDN);
+ MPFR_COPYSIGN (MPC_RE (rop), zero, x, GMP_RNDN);
+ MPFR_COPYSIGN (MPC_IM (rop), zero, y, GMP_RNDN);
mpfr_clear (c);
mpfr_clear (d);
}
-static int
-mpc_div_real (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w, mpc_rnd_t rnd)
-/* Assumes z finite and w finite and non-zero, with imaginary part
- of w a signed zero. */
-{
- int inex_re, inex_im;
- /* save signs of operands in case there are overlaps */
- int zrs = MPFR_SIGNBIT (mpc_realref (z));
- int zis = MPFR_SIGNBIT (mpc_imagref (z));
- int wrs = MPFR_SIGNBIT (mpc_realref (w));
- int wis = MPFR_SIGNBIT (mpc_imagref (w));
-
- /* warning: rop may overlap with z,w so treat the imaginary part first */
- inex_im = mpfr_div (mpc_imagref(rop), mpc_imagref(z), mpc_realref(w), MPC_RND_IM(rnd));
- inex_re = mpfr_div (mpc_realref(rop), mpc_realref(z), mpc_realref(w), MPC_RND_RE(rnd));
-
- /* correct signs of zeroes if necessary, which does not affect the
- inexact flags */
- if (mpfr_zero_p (mpc_realref (rop)))
- mpfr_setsign (mpc_realref (rop), mpc_realref (rop), (zrs != wrs && zis != wis),
- GMP_RNDN); /* exact */
- if (mpfr_zero_p (mpc_imagref (rop)))
- mpfr_setsign (mpc_imagref (rop), mpc_imagref (rop), (zis != wrs && zrs == wis),
- GMP_RNDN);
-
- return MPC_INEX(inex_re, inex_im);
-}
-
-
-static int
-mpc_div_imag (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w, mpc_rnd_t rnd)
-/* Assumes z finite and w finite and non-zero, with real part
- of w a signed zero. */
-{
- int inex_re, inex_im;
- int overlap = (rop == z) || (rop == w);
- int imag_z = mpfr_zero_p (mpc_realref (z));
- mpfr_t wloc;
- mpc_t tmprop;
- mpc_ptr dest = (overlap) ? tmprop : rop;
- /* save signs of operands in case there are overlaps */
- int zrs = MPFR_SIGNBIT (mpc_realref (z));
- int zis = MPFR_SIGNBIT (mpc_imagref (z));
- int wrs = MPFR_SIGNBIT (mpc_realref (w));
- int wis = MPFR_SIGNBIT (mpc_imagref (w));
-
- if (overlap)
- mpc_init3 (tmprop, MPC_PREC_RE (rop), MPC_PREC_IM (rop));
-
- 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);
- /* 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));
-
- if (overlap) {
- /* Note: we could use mpc_swap here, but this might cause problems
- if rop and tmprop have been allocated using different methods, since
- it will swap the significands of rop and tmprop. See
- http://lists.gforge.inria.fr/pipermail/mpc-discuss/2009-August/000504.html */
- mpc_set (rop, tmprop, MPC_RNDNN); /* exact */
- mpc_clear (tmprop);
- }
-
- /* correct signs of zeroes if necessary, which does not affect the
- inexact flags */
- if (mpfr_zero_p (mpc_realref (rop)))
- mpfr_setsign (mpc_realref (rop), mpc_realref (rop), (zrs != wrs && zis != wis),
- GMP_RNDN); /* exact */
- if (imag_z)
- mpfr_setsign (mpc_imagref (rop), mpc_imagref (rop), (zis != wrs && zrs == wis),
- GMP_RNDN);
-
- return MPC_INEX(inex_re, inex_im);
-}
-
-
int
mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd)
{
int ok_re = 0, ok_im = 0;
mpc_t res, c_conj;
mpfr_t q;
- mpfr_prec_t prec;
- int inex, inexact_prod, inexact_norm, inexact_re, inexact_im, loops = 0;
- int underflow_norm, overflow_norm, underflow_prod, overflow_prod;
- int underflow_re = 0, overflow_re = 0, underflow_im = 0, overflow_im = 0;
- mpfr_rnd_t rnd_re = MPC_RND_RE (rnd), rnd_im = MPC_RND_IM (rnd);
- int saved_underflow, saved_overflow;
- int tmpsgn;
+ mp_prec_t prec;
+ int inexact_prod, inexact_norm, inexact_re, inexact_im, loops = 0;
+
+ /* save signs of operands in case there are overlaps */
+ int brs = MPFR_SIGNBIT (MPC_RE (b));
+ int bis = MPFR_SIGNBIT (MPC_IM (b));
+ int crs = MPFR_SIGNBIT (MPC_RE (c));
+ int cis = MPFR_SIGNBIT (MPC_IM (c));
/* According to the C standard G.3, there are three types of numbers: */
/* finite (both parts are usual real numbers; contains 0), infinite */
/* a real; we handle it separately instead. */
if (mpc_zero_p (c))
return mpc_div_zero (a, b, c, rnd);
- else if (mpc_inf_p (b) && mpc_fin_p (c))
+ else {
+ if (mpc_inf_p (b) && mpc_fin_p (c))
return mpc_div_inf_fin (a, b, c);
- else if (mpc_fin_p (b) && mpc_inf_p (c))
+ else if (mpc_fin_p (b) && mpc_inf_p (c))
return mpc_div_fin_inf (a, b, c);
- else if (!mpc_fin_p (b) || !mpc_fin_p (c)) {
- mpc_set_nan (a);
- return MPC_INEX (0, 0);
+ else if (!mpc_fin_p (b) || !mpc_fin_p (c)) {
+ mpfr_set_nan (MPC_RE (a));
+ mpfr_set_nan (MPC_IM (a));
+ return MPC_INEX (0, 0);
+ }
+ }
+
+ /* check for real divisor */
+ if (mpfr_zero_p(MPC_IM(c))) /* (re_b+i*im_b)/c = re_b/c + i * (im_b/c) */
+ {
+ /* warning: a may overlap with b,c so treat the imaginary part first */
+ inexact_im = mpfr_div (MPC_IM(a), MPC_IM(b), MPC_RE(c), MPC_RND_IM(rnd));
+ inexact_re = mpfr_div (MPC_RE(a), MPC_RE(b), MPC_RE(c), MPC_RND_RE(rnd));
+
+ /* correct signs of zeroes if necessary, which does not affect the
+ inexact flags */
+ if (mpfr_zero_p (MPC_RE (a)))
+ mpfr_setsign (MPC_RE (a), MPC_RE (a), (brs != crs && bis != cis),
+ GMP_RNDN); /* exact */
+ if (mpfr_zero_p (MPC_IM (a)))
+ mpfr_setsign (MPC_IM (a), MPC_IM (a), (bis != crs && brs == cis),
+ GMP_RNDN);
+
+ return MPC_INEX(inexact_re, inexact_im);
}
- else if (mpfr_zero_p(mpc_imagref(c)))
- return mpc_div_real (a, b, c, rnd);
- else if (mpfr_zero_p(mpc_realref(c)))
- return mpc_div_imag (a, b, c, rnd);
-
+
+ /* check for purely imaginary divisor */
+ if (mpfr_zero_p(MPC_RE(c)))
+ {
+ /* (re_b+i*im_b)/(i*c) = im_b/c - i * (re_b/c) */
+ int overlap = (a == b) || (a == c);
+ int imag_b = mpfr_zero_p (MPC_RE (b));
+ mpfr_t cloc;
+ mpc_t tmpa;
+ mpc_ptr dest = (overlap) ? tmpa : a;
+
+ if (overlap)
+ mpc_init3 (tmpa, MPFR_PREC (MPC_RE (a)), MPFR_PREC (MPC_IM (a)));
+
+ cloc[0] = MPC_IM(c)[0]; /* copies mpfr struct IM(c) into cloc */
+ inexact_re = mpfr_div (MPC_RE(dest), MPC_IM(b), cloc, MPC_RND_RE(rnd));
+ mpfr_neg (cloc, cloc, GMP_RNDN);
+ /* changes the sign only in cloc, not in c; no need to correct later */
+ inexact_im = mpfr_div (MPC_IM(dest), MPC_RE(b), cloc, MPC_RND_IM(rnd));
+
+ if (overlap)
+ {
+ /* Note: we could use mpc_swap here, but this might cause problems
+ if a and tmpa have been allocated using different methods, since
+ it will swap the significands of a and tmpa. See http://
+ lists.gforge.inria.fr/pipermail/mpc-discuss/2009-August/000504.html */
+ mpc_set (a, tmpa, MPC_RNDNN); /* exact */
+ mpc_clear (tmpa);
+ }
+
+ /* correct signs of zeroes if necessary, which does not affect the
+ inexact flags */
+ if (mpfr_zero_p (MPC_RE (a)))
+ mpfr_setsign (MPC_RE (a), MPC_RE (a), (brs != crs && bis != cis),
+ GMP_RNDN); /* exact */
+ if (imag_b)
+ mpfr_setsign (MPC_IM (a), MPC_IM (a), (bis != crs && brs == cis),
+ GMP_RNDN);
+
+ return MPC_INEX(inexact_re, inexact_im);
+ }
+
prec = MPC_MAX_PREC(a);
mpc_init2 (res, 2);
mpfr_init (q);
/* create the conjugate of c in c_conj without allocating new memory */
- 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));
+ MPC_RE (c_conj)[0] = MPC_RE (c)[0];
+ MPC_IM (c_conj)[0] = MPC_IM (c)[0];
+ MPFR_CHANGE_SIGN (MPC_IM (c_conj));
- /* save the underflow or overflow flags from MPFR */
- saved_underflow = mpfr_underflow_p ();
- saved_overflow = mpfr_overflow_p ();
-
- do {
+ do
+ {
loops ++;
- prec += loops <= 2 ? mpc_ceil_log2 (prec) + 5 : prec / 2;
+ prec += (loops <= 2) ? mpc_ceil_log2 (prec) + 5 : prec / 2;
mpc_set_prec (res, prec);
mpfr_set_prec (q, prec);
- /* first compute norm(c) */
- mpfr_clear_underflow ();
- mpfr_clear_overflow ();
- inexact_norm = mpc_norm (q, c, GMP_RNDU);
- underflow_norm = mpfr_underflow_p ();
- overflow_norm = mpfr_overflow_p ();
- if (underflow_norm)
- mpfr_set_ui (q, 0ul, GMP_RNDN);
- /* to obtain divisions by 0 later on */
+ /* first compute norm(c)^2 */
+ inexact_norm = mpc_norm (q, c, GMP_RNDD);
/* now compute b*conjugate(c) */
- mpfr_clear_underflow ();
- mpfr_clear_overflow ();
+ /* We need rounding away from zero for both the real and the imagin- */
+ /* ary part; then the final result is also rounded away from zero. */
+ /* The error is less than 1 ulp. Since this is not implemented in */
+ /* mpfr, we round towards zero and add 1 ulp to the absolute values */
+ /* if they are not exact. */
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 ();
- overflow_prod = mpfr_overflow_p ();
- /* unfortunately, does not distinguish between under-/overflow
- in real or imaginary parts
- 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));
- isinf = mpfr_inf_p (mpc_realref(res));
- mpfr_nextbelow (mpc_realref(res));
- }
- else if (tmpsgn < 0)
- {
- 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;
- }
- tmpsgn = mpfr_sgn (mpc_imagref(res));
- isinf = 0;
- if (tmpsgn > 0)
- {
- 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));
- isinf = mpfr_inf_p (mpc_imagref(res));
- mpfr_nextabove (mpc_imagref(res));
- }
- if (isinf)
- {
- mpfr_set_inf (mpc_imagref(res), tmpsgn);
- overflow_im = 1;
- }
- mpc_set (a, res, rnd);
- goto end;
- }
+ if (inexact_re != 0)
+ mpfr_add_one_ulp (MPC_RE (res), GMP_RNDN);
+ if (inexact_im != 0)
+ mpfr_add_one_ulp (MPC_IM (res), GMP_RNDN);
/* 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. */
- mpfr_clear_underflow ();
- mpfr_clear_overflow ();
- inexact_re |= mpfr_div (mpc_realref (res), mpc_realref (res), q, GMP_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));
-
- 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);
- 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));
+ 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. */
+ if (MPFR_SIGN (MPC_RE (res)) > 0)
+ {
+ inexact_re |= mpfr_div (MPC_RE (res), MPC_RE (res), q, GMP_RNDU);
+ ok_re = mpfr_can_round (MPC_RE (res), prec - 4, GMP_RNDU,
+ MPC_RND_RE(rnd), MPFR_PREC(MPC_RE(a)));
+ }
+ else
+ {
+ inexact_re |= mpfr_div (MPC_RE (res), MPC_RE (res), q, GMP_RNDD);
+ ok_re = mpfr_can_round (MPC_RE (res), prec - 4, GMP_RNDD,
+ MPC_RND_RE(rnd), MPFR_PREC(MPC_RE(a)));
+ }
+
+ if (ok_re || !inexact_re) /* compute imaginary part */
+ {
+ if (MPFR_SIGN (MPC_IM (res)) > 0)
+ {
+ inexact_im |= mpfr_div (MPC_IM (res), MPC_IM (res), q, GMP_RNDU);
+ ok_im = mpfr_can_round (MPC_IM (res), prec - 4, GMP_RNDU,
+ MPC_RND_IM(rnd), MPFR_PREC(MPC_IM(a)));
+ }
+ else
+ {
+ inexact_im |= mpfr_div (MPC_IM (res), MPC_IM (res), q, GMP_RNDD);
+ ok_im = mpfr_can_round (MPC_IM (res), prec - 4, GMP_RNDD,
+ MPC_RND_IM(rnd), MPFR_PREC(MPC_IM(a)));
+ }
}
}
- else {
+ 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) {
+ /* We do not decide about the sign of the difference. */
+ if (mpfr_ui_div (q, 1, q, GMP_RNDU) || 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 */
- inexact_re |= ~mpfr_zero_p (mpc_realref (res));
- inexact_im |= ~mpfr_zero_p (mpc_imagref (res));
+ inexact_re = inexact_re || !mpfr_zero_p (MPC_RE (res));
+ inexact_im = inexact_im || !mpfr_zero_p (MPC_IM (res));
+ }
+ if (MPFR_SIGN (MPC_RE (res)) > 0)
+ {
+ inexact_re = mpfr_mul (MPC_RE (res), MPC_RE (res), q, GMP_RNDU)
+ || inexact_re;
+ ok_re = mpfr_can_round (MPC_RE (res), prec - 4, GMP_RNDU,
+ MPC_RND_RE(rnd), MPFR_PREC(MPC_RE(a)));
}
- mpfr_clear_underflow ();
- mpfr_clear_overflow ();
- inexact_re |= mpfr_mul (mpc_realref (res), mpc_realref (res), q, GMP_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));
-
- 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);
- 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));
+ else
+ {
+ inexact_re = mpfr_mul (MPC_RE (res), MPC_RE (res), q, GMP_RNDD)
+ || inexact_re;
+ ok_re = mpfr_can_round (MPC_RE (res), prec - 4, GMP_RNDD,
+ MPC_RND_RE(rnd), MPFR_PREC(MPC_RE(a)));
+ }
+
+ if (ok_re) /* compute imaginary part */
+ {
+ if (MPFR_SIGN (MPC_IM (res)) > 0)
+ {
+ inexact_im = mpfr_mul (MPC_IM (res), MPC_IM (res), q, GMP_RNDU)
+ || inexact_im;
+ ok_im = mpfr_can_round (MPC_IM (res), prec - 4, GMP_RNDU,
+ MPC_RND_IM(rnd), MPFR_PREC(MPC_IM(a)));
+ }
+ else
+ {
+ inexact_im = mpfr_mul (MPC_IM (res), MPC_IM (res), q, GMP_RNDD)
+ || inexact_im;
+ ok_im = mpfr_can_round (MPC_IM (res), prec - 4, GMP_RNDD,
+ MPC_RND_IM(rnd), MPFR_PREC(MPC_IM(a)));
+ }
}
}
- } while ((!ok_re || !ok_im) && !underflow_norm && !overflow_norm
- && !underflow_prod && !overflow_prod);
-
- inex = mpc_set (a, res, rnd);
- inexact_re = MPC_INEX_RE (inex);
- inexact_im = MPC_INEX_IM (inex);
-
- end:
- /* fix values and inexact flags in case of overflow/underflow */
- /* FIXME: heuristic, certainly does not cover all cases */
- if (overflow_re || (underflow_norm && !underflow_prod)) {
- mpfr_set_inf (mpc_realref (a), mpfr_sgn (mpc_realref (res)));
- inexact_re = mpfr_sgn (mpc_realref (res));
- }
- else if (underflow_re || (overflow_norm && !overflow_prod)) {
- 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)) {
- inexact_im = mpfr_signbit (mpc_imagref (res)) ? 1 : -1;
- mpfr_set_zero (mpc_imagref (a), -inexact_im);
}
+ while ((!ok_re && inexact_re) || (!ok_im && inexact_im));
+
+ mpc_set (a, res, rnd);
mpc_clear (res);
mpfr_clear (q);
- /* restore underflow and overflow flags from MPFR */
- if (saved_underflow)
- mpfr_set_underflow ();
- if (saved_overflow)
- mpfr_set_overflow ();
-
- return MPC_INEX (inexact_re, inexact_im);
+ return (MPC_INEX (inexact_re, inexact_im));
+ /* Only exactness vs. inexactness is tested, not the sign. */
}
projects / platform / upstream / mpc.git / blobdiff