|
Packit Service |
2e9770 |
/* mpc_div -- Divide two complex numbers.
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
Copyright (C) 2002, 2003, 2004, 2005, 2008, 2009, 2010, 2011, 2012 INRIA
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
This file is part of GNU MPC.
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
GNU MPC is free software; you can redistribute it and/or modify it under
|
|
Packit Service |
2e9770 |
the terms of the GNU Lesser General Public License as published by the
|
|
Packit Service |
2e9770 |
Free Software Foundation; either version 3 of the License, or (at your
|
|
Packit Service |
2e9770 |
option) any later version.
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY
|
|
Packit Service |
2e9770 |
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS
|
|
Packit Service |
2e9770 |
FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for
|
|
Packit Service |
2e9770 |
more details.
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
You should have received a copy of the GNU Lesser General Public License
|
|
Packit Service |
2e9770 |
along with this program. If not, see http://www.gnu.org/licenses/ .
|
|
Packit Service |
2e9770 |
*/
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
#include "mpc-impl.h"
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* this routine deals with the case where w is zero */
|
|
Packit Service |
2e9770 |
static int
|
|
Packit Service |
2e9770 |
mpc_div_zero (mpc_ptr a, mpc_srcptr z, mpc_srcptr w, mpc_rnd_t rnd)
|
|
Packit Service |
2e9770 |
/* Assumes w==0, implementation according to C99 G.5.1.8 */
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
int sign = MPFR_SIGNBIT (mpc_realref (w));
|
|
Packit Service |
2e9770 |
mpfr_t infty;
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
mpfr_init2 (infty, MPFR_PREC_MIN);
|
|
Packit Service |
2e9770 |
mpfr_set_inf (infty, sign);
|
|
Packit Service |
2e9770 |
mpfr_mul (mpc_realref (a), infty, mpc_realref (z), MPC_RND_RE (rnd));
|
|
Packit Service |
2e9770 |
mpfr_mul (mpc_imagref (a), infty, mpc_imagref (z), MPC_RND_IM (rnd));
|
|
Packit Service |
2e9770 |
mpfr_clear (infty);
|
|
Packit Service |
2e9770 |
return MPC_INEX (0, 0); /* exact */
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* this routine deals with the case where z is infinite and w finite */
|
|
Packit Service |
2e9770 |
static int
|
|
Packit Service |
2e9770 |
mpc_div_inf_fin (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w)
|
|
Packit Service |
2e9770 |
/* Assumes w finite and non-zero and z infinite; implementation
|
|
Packit Service |
2e9770 |
according to C99 G.5.1.8 */
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
int a, b, x, y;
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
a = (mpfr_inf_p (mpc_realref (z)) ? MPFR_SIGNBIT (mpc_realref (z)) : 0);
|
|
Packit Service |
2e9770 |
b = (mpfr_inf_p (mpc_imagref (z)) ? MPFR_SIGNBIT (mpc_imagref (z)) : 0);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* a is -1 if Re(z) = -Inf, 1 if Re(z) = +Inf, 0 if Re(z) is finite
|
|
Packit Service |
2e9770 |
b is -1 if Im(z) = -Inf, 1 if Im(z) = +Inf, 0 if Im(z) is finite */
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* x = MPC_MPFR_SIGN (a * mpc_realref (w) + b * mpc_imagref (w)) */
|
|
Packit Service |
2e9770 |
/* y = MPC_MPFR_SIGN (b * mpc_realref (w) - a * mpc_imagref (w)) */
|
|
Packit Service |
2e9770 |
if (a == 0 || b == 0) {
|
|
Packit Service |
2e9770 |
/* only one of a or b can be zero, since z is infinite */
|
|
Packit Service |
2e9770 |
x = a * MPC_MPFR_SIGN (mpc_realref (w)) + b * MPC_MPFR_SIGN (mpc_imagref (w));
|
|
Packit Service |
2e9770 |
y = b * MPC_MPFR_SIGN (mpc_realref (w)) - a * MPC_MPFR_SIGN (mpc_imagref (w));
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else {
|
|
Packit Service |
2e9770 |
/* Both parts of z are infinite; x could be determined by sign
|
|
Packit Service |
2e9770 |
considerations and comparisons. Since operations with non-finite
|
|
Packit Service |
2e9770 |
numbers are not considered time-critical, we let mpfr do the work. */
|
|
Packit Service |
2e9770 |
mpfr_t sign;
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
mpfr_init2 (sign, 2);
|
|
Packit Service |
2e9770 |
/* This is enough to determine the sign of sums and differences. */
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
if (a == 1)
|
|
Packit Service |
2e9770 |
if (b == 1) {
|
|
Packit Service |
2e9770 |
mpfr_add (sign, mpc_realref (w), mpc_imagref (w), MPFR_RNDN);
|
|
Packit Service |
2e9770 |
x = MPC_MPFR_SIGN (sign);
|
|
Packit Service |
2e9770 |
mpfr_sub (sign, mpc_realref (w), mpc_imagref (w), MPFR_RNDN);
|
|
Packit Service |
2e9770 |
y = MPC_MPFR_SIGN (sign);
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else { /* b == -1 */
|
|
Packit Service |
2e9770 |
mpfr_sub (sign, mpc_realref (w), mpc_imagref (w), MPFR_RNDN);
|
|
Packit Service |
2e9770 |
x = MPC_MPFR_SIGN (sign);
|
|
Packit Service |
2e9770 |
mpfr_add (sign, mpc_realref (w), mpc_imagref (w), MPFR_RNDN);
|
|
Packit Service |
2e9770 |
y = -MPC_MPFR_SIGN (sign);
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else /* a == -1 */
|
|
Packit Service |
2e9770 |
if (b == 1) {
|
|
Packit Service |
2e9770 |
mpfr_sub (sign, mpc_imagref (w), mpc_realref (w), MPFR_RNDN);
|
|
Packit Service |
2e9770 |
x = MPC_MPFR_SIGN (sign);
|
|
Packit Service |
2e9770 |
mpfr_add (sign, mpc_realref (w), mpc_imagref (w), MPFR_RNDN);
|
|
Packit Service |
2e9770 |
y = MPC_MPFR_SIGN (sign);
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else { /* b == -1 */
|
|
Packit Service |
2e9770 |
mpfr_add (sign, mpc_realref (w), mpc_imagref (w), MPFR_RNDN);
|
|
Packit Service |
2e9770 |
x = -MPC_MPFR_SIGN (sign);
|
|
Packit Service |
2e9770 |
mpfr_sub (sign, mpc_imagref (w), mpc_realref (w), MPFR_RNDN);
|
|
Packit Service |
2e9770 |
y = MPC_MPFR_SIGN (sign);
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
mpfr_clear (sign);
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
if (x == 0)
|
|
Packit Service |
2e9770 |
mpfr_set_nan (mpc_realref (rop));
|
|
Packit Service |
2e9770 |
else
|
|
Packit Service |
2e9770 |
mpfr_set_inf (mpc_realref (rop), x);
|
|
Packit Service |
2e9770 |
if (y == 0)
|
|
Packit Service |
2e9770 |
mpfr_set_nan (mpc_imagref (rop));
|
|
Packit Service |
2e9770 |
else
|
|
Packit Service |
2e9770 |
mpfr_set_inf (mpc_imagref (rop), y);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
return MPC_INEX (0, 0); /* exact */
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* this routine deals with the case where z if finite and w infinite */
|
|
Packit Service |
2e9770 |
static int
|
|
Packit Service |
2e9770 |
mpc_div_fin_inf (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w)
|
|
Packit Service |
2e9770 |
/* Assumes z finite and w infinite; implementation according to
|
|
Packit Service |
2e9770 |
C99 G.5.1.8 */
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpfr_t c, d, a, b, x, y, zero;
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
mpfr_init2 (c, 2); /* needed to hold a signed zero, +1 or -1 */
|
|
Packit Service |
2e9770 |
mpfr_init2 (d, 2);
|
|
Packit Service |
2e9770 |
mpfr_init2 (x, 2);
|
|
Packit Service |
2e9770 |
mpfr_init2 (y, 2);
|
|
Packit Service |
2e9770 |
mpfr_init2 (zero, 2);
|
|
Packit Service |
2e9770 |
mpfr_set_ui (zero, 0ul, MPFR_RNDN);
|
|
Packit Service |
2e9770 |
mpfr_init2 (a, mpfr_get_prec (mpc_realref (z)));
|
|
Packit Service |
2e9770 |
mpfr_init2 (b, mpfr_get_prec (mpc_imagref (z)));
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
mpfr_set_ui (c, (mpfr_inf_p (mpc_realref (w)) ? 1 : 0), MPFR_RNDN);
|
|
Packit Service |
2e9770 |
MPFR_COPYSIGN (c, c, mpc_realref (w), MPFR_RNDN);
|
|
Packit Service |
2e9770 |
mpfr_set_ui (d, (mpfr_inf_p (mpc_imagref (w)) ? 1 : 0), MPFR_RNDN);
|
|
Packit Service |
2e9770 |
MPFR_COPYSIGN (d, d, mpc_imagref (w), MPFR_RNDN);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
mpfr_mul (a, mpc_realref (z), c, MPFR_RNDN); /* exact */
|
|
Packit Service |
2e9770 |
mpfr_mul (b, mpc_imagref (z), d, MPFR_RNDN);
|
|
Packit Service |
2e9770 |
mpfr_add (x, a, b, MPFR_RNDN);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
mpfr_mul (b, mpc_imagref (z), c, MPFR_RNDN);
|
|
Packit Service |
2e9770 |
mpfr_mul (a, mpc_realref (z), d, MPFR_RNDN);
|
|
Packit Service |
2e9770 |
mpfr_sub (y, b, a, MPFR_RNDN);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
MPFR_COPYSIGN (mpc_realref (rop), zero, x, MPFR_RNDN);
|
|
Packit Service |
2e9770 |
MPFR_COPYSIGN (mpc_imagref (rop), zero, y, MPFR_RNDN);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
mpfr_clear (c);
|
|
Packit Service |
2e9770 |
mpfr_clear (d);
|
|
Packit Service |
2e9770 |
mpfr_clear (x);
|
|
Packit Service |
2e9770 |
mpfr_clear (y);
|
|
Packit Service |
2e9770 |
mpfr_clear (zero);
|
|
Packit Service |
2e9770 |
mpfr_clear (a);
|
|
Packit Service |
2e9770 |
mpfr_clear (b);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
return MPC_INEX (0, 0); /* exact */
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
static int
|
|
Packit Service |
2e9770 |
mpc_div_real (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w, mpc_rnd_t rnd)
|
|
Packit Service |
2e9770 |
/* Assumes z finite and w finite and non-zero, with imaginary part
|
|
Packit Service |
2e9770 |
of w a signed zero. */
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
int inex_re, inex_im;
|
|
Packit Service |
2e9770 |
/* save signs of operands in case there are overlaps */
|
|
Packit Service |
2e9770 |
int zrs = MPFR_SIGNBIT (mpc_realref (z));
|
|
Packit Service |
2e9770 |
int zis = MPFR_SIGNBIT (mpc_imagref (z));
|
|
Packit Service |
2e9770 |
int wrs = MPFR_SIGNBIT (mpc_realref (w));
|
|
Packit Service |
2e9770 |
int wis = MPFR_SIGNBIT (mpc_imagref (w));
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* warning: rop may overlap with z,w so treat the imaginary part first */
|
|
Packit Service |
2e9770 |
inex_im = mpfr_div (mpc_imagref(rop), mpc_imagref(z), mpc_realref(w), MPC_RND_IM(rnd));
|
|
Packit Service |
2e9770 |
inex_re = mpfr_div (mpc_realref(rop), mpc_realref(z), mpc_realref(w), MPC_RND_RE(rnd));
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* correct signs of zeroes if necessary, which does not affect the
|
|
Packit Service |
2e9770 |
inexact flags */
|
|
Packit Service |
2e9770 |
if (mpfr_zero_p (mpc_realref (rop)))
|
|
Packit Service |
2e9770 |
mpfr_setsign (mpc_realref (rop), mpc_realref (rop), (zrs != wrs && zis != wis),
|
|
Packit Service |
2e9770 |
MPFR_RNDN); /* exact */
|
|
Packit Service |
2e9770 |
if (mpfr_zero_p (mpc_imagref (rop)))
|
|
Packit Service |
2e9770 |
mpfr_setsign (mpc_imagref (rop), mpc_imagref (rop), (zis != wrs && zrs == wis),
|
|
Packit Service |
2e9770 |
MPFR_RNDN);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
return MPC_INEX(inex_re, inex_im);
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
static int
|
|
Packit Service |
2e9770 |
mpc_div_imag (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w, mpc_rnd_t rnd)
|
|
Packit Service |
2e9770 |
/* Assumes z finite and w finite and non-zero, with real part
|
|
Packit Service |
2e9770 |
of w a signed zero. */
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
int inex_re, inex_im;
|
|
Packit Service |
2e9770 |
int overlap = (rop == z) || (rop == w);
|
|
Packit Service |
2e9770 |
int imag_z = mpfr_zero_p (mpc_realref (z));
|
|
Packit Service |
2e9770 |
mpfr_t wloc;
|
|
Packit Service |
2e9770 |
mpc_t tmprop;
|
|
Packit Service |
2e9770 |
mpc_ptr dest = (overlap) ? tmprop : rop;
|
|
Packit Service |
2e9770 |
/* save signs of operands in case there are overlaps */
|
|
Packit Service |
2e9770 |
int zrs = MPFR_SIGNBIT (mpc_realref (z));
|
|
Packit Service |
2e9770 |
int zis = MPFR_SIGNBIT (mpc_imagref (z));
|
|
Packit Service |
2e9770 |
int wrs = MPFR_SIGNBIT (mpc_realref (w));
|
|
Packit Service |
2e9770 |
int wis = MPFR_SIGNBIT (mpc_imagref (w));
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
if (overlap)
|
|
Packit Service |
2e9770 |
mpc_init3 (tmprop, MPC_PREC_RE (rop), MPC_PREC_IM (rop));
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
wloc[0] = mpc_imagref(w)[0]; /* copies mpfr struct IM(w) into wloc */
|
|
Packit Service |
2e9770 |
inex_re = mpfr_div (mpc_realref(dest), mpc_imagref(z), wloc, MPC_RND_RE(rnd));
|
|
Packit Service |
2e9770 |
mpfr_neg (wloc, wloc, MPFR_RNDN);
|
|
Packit Service |
2e9770 |
/* changes the sign only in wloc, not in w; no need to correct later */
|
|
Packit Service |
2e9770 |
inex_im = mpfr_div (mpc_imagref(dest), mpc_realref(z), wloc, MPC_RND_IM(rnd));
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
if (overlap) {
|
|
Packit Service |
2e9770 |
/* Note: we could use mpc_swap here, but this might cause problems
|
|
Packit Service |
2e9770 |
if rop and tmprop have been allocated using different methods, since
|
|
Packit Service |
2e9770 |
it will swap the significands of rop and tmprop. See
|
|
Packit Service |
2e9770 |
http://lists.gforge.inria.fr/pipermail/mpc-discuss/2009-August/000504.html */
|
|
Packit Service |
2e9770 |
mpc_set (rop, tmprop, MPC_RNDNN); /* exact */
|
|
Packit Service |
2e9770 |
mpc_clear (tmprop);
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* correct signs of zeroes if necessary, which does not affect the
|
|
Packit Service |
2e9770 |
inexact flags */
|
|
Packit Service |
2e9770 |
if (mpfr_zero_p (mpc_realref (rop)))
|
|
Packit Service |
2e9770 |
mpfr_setsign (mpc_realref (rop), mpc_realref (rop), (zrs != wrs && zis != wis),
|
|
Packit Service |
2e9770 |
MPFR_RNDN); /* exact */
|
|
Packit Service |
2e9770 |
if (imag_z)
|
|
Packit Service |
2e9770 |
mpfr_setsign (mpc_imagref (rop), mpc_imagref (rop), (zis != wrs && zrs == wis),
|
|
Packit Service |
2e9770 |
MPFR_RNDN);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
return MPC_INEX(inex_re, inex_im);
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
int
|
|
Packit Service |
2e9770 |
mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
int ok_re = 0, ok_im = 0;
|
|
Packit Service |
2e9770 |
mpc_t res, c_conj;
|
|
Packit Service |
2e9770 |
mpfr_t q;
|
|
Packit Service |
2e9770 |
mpfr_prec_t prec;
|
|
Packit Service |
2e9770 |
int inex, inexact_prod, inexact_norm, inexact_re, inexact_im, loops = 0;
|
|
Packit Service |
2e9770 |
int underflow_norm, overflow_norm, underflow_prod, overflow_prod;
|
|
Packit Service |
2e9770 |
int underflow_re = 0, overflow_re = 0, underflow_im = 0, overflow_im = 0;
|
|
Packit Service |
2e9770 |
mpfr_rnd_t rnd_re = MPC_RND_RE (rnd), rnd_im = MPC_RND_IM (rnd);
|
|
Packit Service |
2e9770 |
int saved_underflow, saved_overflow;
|
|
Packit Service |
2e9770 |
int tmpsgn;
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* According to the C standard G.3, there are three types of numbers: */
|
|
Packit Service |
2e9770 |
/* finite (both parts are usual real numbers; contains 0), infinite */
|
|
Packit Service |
2e9770 |
/* (at least one part is a real infinity) and all others; the latter */
|
|
Packit Service |
2e9770 |
/* are numbers containing a nan, but no infinity, and could reasonably */
|
|
Packit Service |
2e9770 |
/* be called nan. */
|
|
Packit Service |
2e9770 |
/* By G.5.1.4, infinite/finite=infinite; finite/infinite=0; */
|
|
Packit Service |
2e9770 |
/* all other divisions that are not finite/finite return nan+i*nan. */
|
|
Packit Service |
2e9770 |
/* Division by 0 could be handled by the following case of division by */
|
|
Packit Service |
2e9770 |
/* a real; we handle it separately instead. */
|
|
Packit Service |
2e9770 |
if (mpc_zero_p (c))
|
|
Packit Service |
2e9770 |
return mpc_div_zero (a, b, c, rnd);
|
|
Packit Service |
2e9770 |
else if (mpc_inf_p (b) && mpc_fin_p (c))
|
|
Packit Service |
2e9770 |
return mpc_div_inf_fin (a, b, c);
|
|
Packit Service |
2e9770 |
else if (mpc_fin_p (b) && mpc_inf_p (c))
|
|
Packit Service |
2e9770 |
return mpc_div_fin_inf (a, b, c);
|
|
Packit Service |
2e9770 |
else if (!mpc_fin_p (b) || !mpc_fin_p (c)) {
|
|
Packit Service |
2e9770 |
mpc_set_nan (a);
|
|
Packit Service |
2e9770 |
return MPC_INEX (0, 0);
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else if (mpfr_zero_p(mpc_imagref(c)))
|
|
Packit Service |
2e9770 |
return mpc_div_real (a, b, c, rnd);
|
|
Packit Service |
2e9770 |
else if (mpfr_zero_p(mpc_realref(c)))
|
|
Packit Service |
2e9770 |
return mpc_div_imag (a, b, c, rnd);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
prec = MPC_MAX_PREC(a);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
mpc_init2 (res, 2);
|
|
Packit Service |
2e9770 |
mpfr_init (q);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* create the conjugate of c in c_conj without allocating new memory */
|
|
Packit Service |
2e9770 |
mpc_realref (c_conj)[0] = mpc_realref (c)[0];
|
|
Packit Service |
2e9770 |
mpc_imagref (c_conj)[0] = mpc_imagref (c)[0];
|
|
Packit Service |
2e9770 |
MPFR_CHANGE_SIGN (mpc_imagref (c_conj));
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* save the underflow or overflow flags from MPFR */
|
|
Packit Service |
2e9770 |
saved_underflow = mpfr_underflow_p ();
|
|
Packit Service |
2e9770 |
saved_overflow = mpfr_overflow_p ();
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
do {
|
|
Packit Service |
2e9770 |
loops ++;
|
|
Packit Service |
2e9770 |
prec += loops <= 2 ? mpc_ceil_log2 (prec) + 5 : prec / 2;
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
mpc_set_prec (res, prec);
|
|
Packit Service |
2e9770 |
mpfr_set_prec (q, prec);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* first compute norm(c) */
|
|
Packit Service |
2e9770 |
mpfr_clear_underflow ();
|
|
Packit Service |
2e9770 |
mpfr_clear_overflow ();
|
|
Packit Service |
2e9770 |
inexact_norm = mpc_norm (q, c, MPFR_RNDU);
|
|
Packit Service |
2e9770 |
underflow_norm = mpfr_underflow_p ();
|
|
Packit Service |
2e9770 |
overflow_norm = mpfr_overflow_p ();
|
|
Packit Service |
2e9770 |
if (underflow_norm)
|
|
Packit Service |
2e9770 |
mpfr_set_ui (q, 0ul, MPFR_RNDN);
|
|
Packit Service |
2e9770 |
/* to obtain divisions by 0 later on */
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* now compute b*conjugate(c) */
|
|
Packit Service |
2e9770 |
mpfr_clear_underflow ();
|
|
Packit Service |
2e9770 |
mpfr_clear_overflow ();
|
|
Packit Service |
2e9770 |
inexact_prod = mpc_mul (res, b, c_conj, MPC_RNDZZ);
|
|
Packit Service |
2e9770 |
inexact_re = MPC_INEX_RE (inexact_prod);
|
|
Packit Service |
2e9770 |
inexact_im = MPC_INEX_IM (inexact_prod);
|
|
Packit Service |
2e9770 |
underflow_prod = mpfr_underflow_p ();
|
|
Packit Service |
2e9770 |
overflow_prod = mpfr_overflow_p ();
|
|
Packit Service |
2e9770 |
/* unfortunately, does not distinguish between under-/overflow
|
|
Packit Service |
2e9770 |
in real or imaginary parts
|
|
Packit Service |
2e9770 |
hopefully, the side-effects of mpc_mul do indeed raise the
|
|
Packit Service |
2e9770 |
mpfr exceptions */
|
|
Packit Service |
2e9770 |
if (overflow_prod) {
|
|
Packit Service |
2e9770 |
int isinf = 0;
|
|
Packit Service |
2e9770 |
tmpsgn = mpfr_sgn (mpc_realref(res));
|
|
Packit Service |
2e9770 |
if (tmpsgn > 0)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpfr_nextabove (mpc_realref(res));
|
|
Packit Service |
2e9770 |
isinf = mpfr_inf_p (mpc_realref(res));
|
|
Packit Service |
2e9770 |
mpfr_nextbelow (mpc_realref(res));
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else if (tmpsgn < 0)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpfr_nextbelow (mpc_realref(res));
|
|
Packit Service |
2e9770 |
isinf = mpfr_inf_p (mpc_realref(res));
|
|
Packit Service |
2e9770 |
mpfr_nextabove (mpc_realref(res));
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
if (isinf)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpfr_set_inf (mpc_realref(res), tmpsgn);
|
|
Packit Service |
2e9770 |
overflow_re = 1;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
tmpsgn = mpfr_sgn (mpc_imagref(res));
|
|
Packit Service |
2e9770 |
isinf = 0;
|
|
Packit Service |
2e9770 |
if (tmpsgn > 0)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpfr_nextabove (mpc_imagref(res));
|
|
Packit Service |
2e9770 |
isinf = mpfr_inf_p (mpc_imagref(res));
|
|
Packit Service |
2e9770 |
mpfr_nextbelow (mpc_imagref(res));
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else if (tmpsgn < 0)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpfr_nextbelow (mpc_imagref(res));
|
|
Packit Service |
2e9770 |
isinf = mpfr_inf_p (mpc_imagref(res));
|
|
Packit Service |
2e9770 |
mpfr_nextabove (mpc_imagref(res));
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
if (isinf)
|
|
Packit Service |
2e9770 |
{
|
|
Packit Service |
2e9770 |
mpfr_set_inf (mpc_imagref(res), tmpsgn);
|
|
Packit Service |
2e9770 |
overflow_im = 1;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
mpc_set (a, res, rnd);
|
|
Packit Service |
2e9770 |
goto end;
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* divide the product by the norm */
|
|
Packit Service |
2e9770 |
if (inexact_norm == 0 && (inexact_re == 0 || inexact_im == 0)) {
|
|
Packit Service |
2e9770 |
/* The division has good chances to be exact in at least one part. */
|
|
Packit Service |
2e9770 |
/* Since this can cause problems when not rounding to the nearest, */
|
|
Packit Service |
2e9770 |
/* we use the division code of mpfr, which handles the situation. */
|
|
Packit Service |
2e9770 |
mpfr_clear_underflow ();
|
|
Packit Service |
2e9770 |
mpfr_clear_overflow ();
|
|
Packit Service |
2e9770 |
inexact_re |= mpfr_div (mpc_realref (res), mpc_realref (res), q, MPFR_RNDZ);
|
|
Packit Service |
2e9770 |
underflow_re = mpfr_underflow_p ();
|
|
Packit Service |
2e9770 |
overflow_re = mpfr_overflow_p ();
|
|
Packit Service |
2e9770 |
ok_re = !inexact_re || underflow_re || overflow_re
|
|
Packit Service |
2e9770 |
|| mpfr_can_round (mpc_realref (res), prec - 4, MPFR_RNDN,
|
|
Packit Service |
2e9770 |
MPFR_RNDZ, MPC_PREC_RE(a) + (rnd_re == MPFR_RNDN));
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
if (ok_re) /* compute imaginary part */ {
|
|
Packit Service |
2e9770 |
mpfr_clear_underflow ();
|
|
Packit Service |
2e9770 |
mpfr_clear_overflow ();
|
|
Packit Service |
2e9770 |
inexact_im |= mpfr_div (mpc_imagref (res), mpc_imagref (res), q, MPFR_RNDZ);
|
|
Packit Service |
2e9770 |
underflow_im = mpfr_underflow_p ();
|
|
Packit Service |
2e9770 |
overflow_im = mpfr_overflow_p ();
|
|
Packit Service |
2e9770 |
ok_im = !inexact_im || underflow_im || overflow_im
|
|
Packit Service |
2e9770 |
|| mpfr_can_round (mpc_imagref (res), prec - 4, MPFR_RNDN,
|
|
Packit Service |
2e9770 |
MPFR_RNDZ, MPC_PREC_IM(a) + (rnd_im == MPFR_RNDN));
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else {
|
|
Packit Service |
2e9770 |
/* The division is inexact, so for efficiency reasons we invert q */
|
|
Packit Service |
2e9770 |
/* only once and multiply by the inverse. */
|
|
Packit Service |
2e9770 |
if (mpfr_ui_div (q, 1ul, q, MPFR_RNDZ) || inexact_norm) {
|
|
Packit Service |
2e9770 |
/* if 1/q is inexact, the approximations of the real and
|
|
Packit Service |
2e9770 |
imaginary part below will be inexact, unless RE(res)
|
|
Packit Service |
2e9770 |
or IM(res) is zero */
|
|
Packit Service |
2e9770 |
inexact_re |= !mpfr_zero_p (mpc_realref (res));
|
|
Packit Service |
2e9770 |
inexact_im |= !mpfr_zero_p (mpc_imagref (res));
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
mpfr_clear_underflow ();
|
|
Packit Service |
2e9770 |
mpfr_clear_overflow ();
|
|
Packit Service |
2e9770 |
inexact_re |= mpfr_mul (mpc_realref (res), mpc_realref (res), q, MPFR_RNDZ);
|
|
Packit Service |
2e9770 |
underflow_re = mpfr_underflow_p ();
|
|
Packit Service |
2e9770 |
overflow_re = mpfr_overflow_p ();
|
|
Packit Service |
2e9770 |
ok_re = !inexact_re || underflow_re || overflow_re
|
|
Packit Service |
2e9770 |
|| mpfr_can_round (mpc_realref (res), prec - 4, MPFR_RNDN,
|
|
Packit Service |
2e9770 |
MPFR_RNDZ, MPC_PREC_RE(a) + (rnd_re == MPFR_RNDN));
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
if (ok_re) /* compute imaginary part */ {
|
|
Packit Service |
2e9770 |
mpfr_clear_underflow ();
|
|
Packit Service |
2e9770 |
mpfr_clear_overflow ();
|
|
Packit Service |
2e9770 |
inexact_im |= mpfr_mul (mpc_imagref (res), mpc_imagref (res), q, MPFR_RNDZ);
|
|
Packit Service |
2e9770 |
underflow_im = mpfr_underflow_p ();
|
|
Packit Service |
2e9770 |
overflow_im = mpfr_overflow_p ();
|
|
Packit Service |
2e9770 |
ok_im = !inexact_im || underflow_im || overflow_im
|
|
Packit Service |
2e9770 |
|| mpfr_can_round (mpc_imagref (res), prec - 4, MPFR_RNDN,
|
|
Packit Service |
2e9770 |
MPFR_RNDZ, MPC_PREC_IM(a) + (rnd_im == MPFR_RNDN));
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
} while ((!ok_re || !ok_im) && !underflow_norm && !overflow_norm
|
|
Packit Service |
2e9770 |
&& !underflow_prod && !overflow_prod);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
inex = mpc_set (a, res, rnd);
|
|
Packit Service |
2e9770 |
inexact_re = MPC_INEX_RE (inex);
|
|
Packit Service |
2e9770 |
inexact_im = MPC_INEX_IM (inex);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
end:
|
|
Packit Service |
2e9770 |
/* fix values and inexact flags in case of overflow/underflow */
|
|
Packit Service |
2e9770 |
/* FIXME: heuristic, certainly does not cover all cases */
|
|
Packit Service |
2e9770 |
if (overflow_re || (underflow_norm && !underflow_prod)) {
|
|
Packit Service |
2e9770 |
mpfr_set_inf (mpc_realref (a), mpfr_sgn (mpc_realref (res)));
|
|
Packit Service |
2e9770 |
inexact_re = mpfr_sgn (mpc_realref (res));
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else if (underflow_re || (overflow_norm && !overflow_prod)) {
|
|
Packit Service |
2e9770 |
inexact_re = mpfr_signbit (mpc_realref (res)) ? 1 : -1;
|
|
Packit Service |
2e9770 |
mpfr_set_zero (mpc_realref (a), -inexact_re);
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
if (overflow_im || (underflow_norm && !underflow_prod)) {
|
|
Packit Service |
2e9770 |
mpfr_set_inf (mpc_imagref (a), mpfr_sgn (mpc_imagref (res)));
|
|
Packit Service |
2e9770 |
inexact_im = mpfr_sgn (mpc_imagref (res));
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
else if (underflow_im || (overflow_norm && !overflow_prod)) {
|
|
Packit Service |
2e9770 |
inexact_im = mpfr_signbit (mpc_imagref (res)) ? 1 : -1;
|
|
Packit Service |
2e9770 |
mpfr_set_zero (mpc_imagref (a), -inexact_im);
|
|
Packit Service |
2e9770 |
}
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
mpc_clear (res);
|
|
Packit Service |
2e9770 |
mpfr_clear (q);
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
/* restore underflow and overflow flags from MPFR */
|
|
Packit Service |
2e9770 |
if (saved_underflow)
|
|
Packit Service |
2e9770 |
mpfr_set_underflow ();
|
|
Packit Service |
2e9770 |
if (saved_overflow)
|
|
Packit Service |
2e9770 |
mpfr_set_overflow ();
|
|
Packit Service |
2e9770 |
|
|
Packit Service |
2e9770 |
return MPC_INEX (inexact_re, inexact_im);
|
|
Packit Service |
2e9770 |
}
|