1/* mpc_div -- Divide two complex numbers. 2 3Copyright (C) 2002, 2003, 2004, 2005, 2008, 2009, 2010, 2011, 2012 INRIA 4 5This file is part of GNU MPC. 6 7GNU MPC is free software; you can redistribute it and/or modify it under 8the terms of the GNU Lesser General Public License as published by the 9Free Software Foundation; either version 3 of the License, or (at your 10option) any later version. 11 12GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY 13WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 14FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for 15more details. 16 17You should have received a copy of the GNU Lesser General Public License 18along with this program. If not, see http://www.gnu.org/licenses/ . 19*/ 20 21#include "mpc-impl.h" 22 23/* this routine deals with the case where w is zero */ 24static int 25mpc_div_zero (mpc_ptr a, mpc_srcptr z, mpc_srcptr w, mpc_rnd_t rnd) 26/* Assumes w==0, implementation according to C99 G.5.1.8 */ 27{ 28 int sign = MPFR_SIGNBIT (mpc_realref (w)); 29 mpfr_t infty; 30 31 mpfr_init2 (infty, MPFR_PREC_MIN); 32 mpfr_set_inf (infty, sign); 33 mpfr_mul (mpc_realref (a), infty, mpc_realref (z), MPC_RND_RE (rnd)); 34 mpfr_mul (mpc_imagref (a), infty, mpc_imagref (z), MPC_RND_IM (rnd)); 35 mpfr_clear (infty); 36 return MPC_INEX (0, 0); /* exact */ 37} 38 39/* this routine deals with the case where z is infinite and w finite */ 40static int 41mpc_div_inf_fin (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w) 42/* Assumes w finite and non-zero and z infinite; implementation 43 according to C99 G.5.1.8 */ 44{ 45 int a, b, x, y; 46 47 a = (mpfr_inf_p (mpc_realref (z)) ? MPFR_SIGNBIT (mpc_realref (z)) : 0); 48 b = (mpfr_inf_p (mpc_imagref (z)) ? MPFR_SIGNBIT (mpc_imagref (z)) : 0); 49 50 /* a is -1 if Re(z) = -Inf, 1 if Re(z) = +Inf, 0 if Re(z) is finite 51 b is -1 if Im(z) = -Inf, 1 if Im(z) = +Inf, 0 if Im(z) is finite */ 52 53 /* x = MPC_MPFR_SIGN (a * mpc_realref (w) + b * mpc_imagref (w)) */ 54 /* y = MPC_MPFR_SIGN (b * mpc_realref (w) - a * mpc_imagref (w)) */ 55 if (a == 0 || b == 0) { 56 /* only one of a or b can be zero, since z is infinite */ 57 x = a * MPC_MPFR_SIGN (mpc_realref (w)) + b * MPC_MPFR_SIGN (mpc_imagref (w)); 58 y = b * MPC_MPFR_SIGN (mpc_realref (w)) - a * MPC_MPFR_SIGN (mpc_imagref (w)); 59 } 60 else { 61 /* Both parts of z are infinite; x could be determined by sign 62 considerations and comparisons. Since operations with non-finite 63 numbers are not considered time-critical, we let mpfr do the work. */ 64 mpfr_t sign; 65 66 mpfr_init2 (sign, 2); 67 /* This is enough to determine the sign of sums and differences. */ 68 69 if (a == 1) 70 if (b == 1) { 71 mpfr_add (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN); 72 x = MPC_MPFR_SIGN (sign); 73 mpfr_sub (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN); 74 y = MPC_MPFR_SIGN (sign); 75 } 76 else { /* b == -1 */ 77 mpfr_sub (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN); 78 x = MPC_MPFR_SIGN (sign); 79 mpfr_add (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN); 80 y = -MPC_MPFR_SIGN (sign); 81 } 82 else /* a == -1 */ 83 if (b == 1) { 84 mpfr_sub (sign, mpc_imagref (w), mpc_realref (w), GMP_RNDN); 85 x = MPC_MPFR_SIGN (sign); 86 mpfr_add (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN); 87 y = MPC_MPFR_SIGN (sign); 88 } 89 else { /* b == -1 */ 90 mpfr_add (sign, mpc_realref (w), mpc_imagref (w), GMP_RNDN); 91 x = -MPC_MPFR_SIGN (sign); 92 mpfr_sub (sign, mpc_imagref (w), mpc_realref (w), GMP_RNDN); 93 y = MPC_MPFR_SIGN (sign); 94 } 95 mpfr_clear (sign); 96 } 97 98 if (x == 0) 99 mpfr_set_nan (mpc_realref (rop)); 100 else 101 mpfr_set_inf (mpc_realref (rop), x); 102 if (y == 0) 103 mpfr_set_nan (mpc_imagref (rop)); 104 else 105 mpfr_set_inf (mpc_imagref (rop), y); 106 107 return MPC_INEX (0, 0); /* exact */ 108} 109 110 111/* this routine deals with the case where z if finite and w infinite */ 112static int 113mpc_div_fin_inf (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w) 114/* Assumes z finite and w infinite; implementation according to 115 C99 G.5.1.8 */ 116{ 117 mpfr_t c, d, a, b, x, y, zero; 118 119 mpfr_init2 (c, 2); /* needed to hold a signed zero, +1 or -1 */ 120 mpfr_init2 (d, 2); 121 mpfr_init2 (x, 2); 122 mpfr_init2 (y, 2); 123 mpfr_init2 (zero, 2); 124 mpfr_set_ui (zero, 0ul, GMP_RNDN); 125 mpfr_init2 (a, mpfr_get_prec (mpc_realref (z))); 126 mpfr_init2 (b, mpfr_get_prec (mpc_imagref (z))); 127 128 mpfr_set_ui (c, (mpfr_inf_p (mpc_realref (w)) ? 1 : 0), GMP_RNDN); 129 MPFR_COPYSIGN (c, c, mpc_realref (w), GMP_RNDN); 130 mpfr_set_ui (d, (mpfr_inf_p (mpc_imagref (w)) ? 1 : 0), GMP_RNDN); 131 MPFR_COPYSIGN (d, d, mpc_imagref (w), GMP_RNDN); 132 133 mpfr_mul (a, mpc_realref (z), c, GMP_RNDN); /* exact */ 134 mpfr_mul (b, mpc_imagref (z), d, GMP_RNDN); 135 mpfr_add (x, a, b, GMP_RNDN); 136 137 mpfr_mul (b, mpc_imagref (z), c, GMP_RNDN); 138 mpfr_mul (a, mpc_realref (z), d, GMP_RNDN); 139 mpfr_sub (y, b, a, GMP_RNDN); 140 141 MPFR_COPYSIGN (mpc_realref (rop), zero, x, GMP_RNDN); 142 MPFR_COPYSIGN (mpc_imagref (rop), zero, y, GMP_RNDN); 143 144 mpfr_clear (c); 145 mpfr_clear (d); 146 mpfr_clear (x); 147 mpfr_clear (y); 148 mpfr_clear (zero); 149 mpfr_clear (a); 150 mpfr_clear (b); 151 152 return MPC_INEX (0, 0); /* exact */ 153} 154 155 156static int 157mpc_div_real (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w, mpc_rnd_t rnd) 158/* Assumes z finite and w finite and non-zero, with imaginary part 159 of w a signed zero. */ 160{ 161 int inex_re, inex_im; 162 /* save signs of operands in case there are overlaps */ 163 int zrs = MPFR_SIGNBIT (mpc_realref (z)); 164 int zis = MPFR_SIGNBIT (mpc_imagref (z)); 165 int wrs = MPFR_SIGNBIT (mpc_realref (w)); 166 int wis = MPFR_SIGNBIT (mpc_imagref (w)); 167 168 /* warning: rop may overlap with z,w so treat the imaginary part first */ 169 inex_im = mpfr_div (mpc_imagref(rop), mpc_imagref(z), mpc_realref(w), MPC_RND_IM(rnd)); 170 inex_re = mpfr_div (mpc_realref(rop), mpc_realref(z), mpc_realref(w), MPC_RND_RE(rnd)); 171 172 /* correct signs of zeroes if necessary, which does not affect the 173 inexact flags */ 174 if (mpfr_zero_p (mpc_realref (rop))) 175 mpfr_setsign (mpc_realref (rop), mpc_realref (rop), (zrs != wrs && zis != wis), 176 GMP_RNDN); /* exact */ 177 if (mpfr_zero_p (mpc_imagref (rop))) 178 mpfr_setsign (mpc_imagref (rop), mpc_imagref (rop), (zis != wrs && zrs == wis), 179 GMP_RNDN); 180 181 return MPC_INEX(inex_re, inex_im); 182} 183 184 185static int 186mpc_div_imag (mpc_ptr rop, mpc_srcptr z, mpc_srcptr w, mpc_rnd_t rnd) 187/* Assumes z finite and w finite and non-zero, with real part 188 of w a signed zero. */ 189{ 190 int inex_re, inex_im; 191 int overlap = (rop == z) || (rop == w); 192 int imag_z = mpfr_zero_p (mpc_realref (z)); 193 mpfr_t wloc; 194 mpc_t tmprop; 195 mpc_ptr dest = (overlap) ? tmprop : rop; 196 /* save signs of operands in case there are overlaps */ 197 int zrs = MPFR_SIGNBIT (mpc_realref (z)); 198 int zis = MPFR_SIGNBIT (mpc_imagref (z)); 199 int wrs = MPFR_SIGNBIT (mpc_realref (w)); 200 int wis = MPFR_SIGNBIT (mpc_imagref (w)); 201 202 if (overlap) 203 mpc_init3 (tmprop, MPC_PREC_RE (rop), MPC_PREC_IM (rop)); 204 205 wloc[0] = mpc_imagref(w)[0]; /* copies mpfr struct IM(w) into wloc */ 206 inex_re = mpfr_div (mpc_realref(dest), mpc_imagref(z), wloc, MPC_RND_RE(rnd)); 207 mpfr_neg (wloc, wloc, GMP_RNDN); 208 /* changes the sign only in wloc, not in w; no need to correct later */ 209 inex_im = mpfr_div (mpc_imagref(dest), mpc_realref(z), wloc, MPC_RND_IM(rnd)); 210 211 if (overlap) { 212 /* Note: we could use mpc_swap here, but this might cause problems 213 if rop and tmprop have been allocated using different methods, since 214 it will swap the significands of rop and tmprop. See 215 http://lists.gforge.inria.fr/pipermail/mpc-discuss/2009-August/000504.html */ 216 mpc_set (rop, tmprop, MPC_RNDNN); /* exact */ 217 mpc_clear (tmprop); 218 } 219 220 /* correct signs of zeroes if necessary, which does not affect the 221 inexact flags */ 222 if (mpfr_zero_p (mpc_realref (rop))) 223 mpfr_setsign (mpc_realref (rop), mpc_realref (rop), (zrs != wrs && zis != wis), 224 GMP_RNDN); /* exact */ 225 if (imag_z) 226 mpfr_setsign (mpc_imagref (rop), mpc_imagref (rop), (zis != wrs && zrs == wis), 227 GMP_RNDN); 228 229 return MPC_INEX(inex_re, inex_im); 230} 231 232 233int 234mpc_div (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd) 235{ 236 int ok_re = 0, ok_im = 0; 237 mpc_t res, c_conj; 238 mpfr_t q; 239 mpfr_prec_t prec; 240 int inex, inexact_prod, inexact_norm, inexact_re, inexact_im, loops = 0; 241 int underflow_norm, overflow_norm, underflow_prod, overflow_prod; 242 int underflow_re = 0, overflow_re = 0, underflow_im = 0, overflow_im = 0; 243 mpfr_rnd_t rnd_re = MPC_RND_RE (rnd), rnd_im = MPC_RND_IM (rnd); 244 int saved_underflow, saved_overflow; 245 int tmpsgn; 246 247 /* According to the C standard G.3, there are three types of numbers: */ 248 /* finite (both parts are usual real numbers; contains 0), infinite */ 249 /* (at least one part is a real infinity) and all others; the latter */ 250 /* are numbers containing a nan, but no infinity, and could reasonably */ 251 /* be called nan. */ 252 /* By G.5.1.4, infinite/finite=infinite; finite/infinite=0; */ 253 /* all other divisions that are not finite/finite return nan+i*nan. */ 254 /* Division by 0 could be handled by the following case of division by */ 255 /* a real; we handle it separately instead. */ 256 if (mpc_zero_p (c)) 257 return mpc_div_zero (a, b, c, rnd); 258 else if (mpc_inf_p (b) && mpc_fin_p (c)) 259 return mpc_div_inf_fin (a, b, c); 260 else if (mpc_fin_p (b) && mpc_inf_p (c)) 261 return mpc_div_fin_inf (a, b, c); 262 else if (!mpc_fin_p (b) || !mpc_fin_p (c)) { 263 mpc_set_nan (a); 264 return MPC_INEX (0, 0); 265 } 266 else if (mpfr_zero_p(mpc_imagref(c))) 267 return mpc_div_real (a, b, c, rnd); 268 else if (mpfr_zero_p(mpc_realref(c))) 269 return mpc_div_imag (a, b, c, rnd); 270 271 prec = MPC_MAX_PREC(a); 272 273 mpc_init2 (res, 2); 274 mpfr_init (q); 275 276 /* create the conjugate of c in c_conj without allocating new memory */ 277 mpc_realref (c_conj)[0] = mpc_realref (c)[0]; 278 mpc_imagref (c_conj)[0] = mpc_imagref (c)[0]; 279 MPFR_CHANGE_SIGN (mpc_imagref (c_conj)); 280 281 /* save the underflow or overflow flags from MPFR */ 282 saved_underflow = mpfr_underflow_p (); 283 saved_overflow = mpfr_overflow_p (); 284 285 do { 286 loops ++; 287 prec += loops <= 2 ? mpc_ceil_log2 (prec) + 5 : prec / 2; 288 289 mpc_set_prec (res, prec); 290 mpfr_set_prec (q, prec); 291 292 /* first compute norm(c) */ 293 mpfr_clear_underflow (); 294 mpfr_clear_overflow (); 295 inexact_norm = mpc_norm (q, c, GMP_RNDU); 296 underflow_norm = mpfr_underflow_p (); 297 overflow_norm = mpfr_overflow_p (); 298 if (underflow_norm) 299 mpfr_set_ui (q, 0ul, GMP_RNDN); 300 /* to obtain divisions by 0 later on */ 301 302 /* now compute b*conjugate(c) */ 303 mpfr_clear_underflow (); 304 mpfr_clear_overflow (); 305 inexact_prod = mpc_mul (res, b, c_conj, MPC_RNDZZ); 306 inexact_re = MPC_INEX_RE (inexact_prod); 307 inexact_im = MPC_INEX_IM (inexact_prod); 308 underflow_prod = mpfr_underflow_p (); 309 overflow_prod = mpfr_overflow_p (); 310 /* unfortunately, does not distinguish between under-/overflow 311 in real or imaginary parts 312 hopefully, the side-effects of mpc_mul do indeed raise the 313 mpfr exceptions */ 314 if (overflow_prod) { 315 int isinf = 0; 316 tmpsgn = mpfr_sgn (mpc_realref(res)); 317 if (tmpsgn > 0) 318 { 319 mpfr_nextabove (mpc_realref(res)); 320 isinf = mpfr_inf_p (mpc_realref(res)); 321 mpfr_nextbelow (mpc_realref(res)); 322 } 323 else if (tmpsgn < 0) 324 { 325 mpfr_nextbelow (mpc_realref(res)); 326 isinf = mpfr_inf_p (mpc_realref(res)); 327 mpfr_nextabove (mpc_realref(res)); 328 } 329 if (isinf) 330 { 331 mpfr_set_inf (mpc_realref(res), tmpsgn); 332 overflow_re = 1; 333 } 334 tmpsgn = mpfr_sgn (mpc_imagref(res)); 335 isinf = 0; 336 if (tmpsgn > 0) 337 { 338 mpfr_nextabove (mpc_imagref(res)); 339 isinf = mpfr_inf_p (mpc_imagref(res)); 340 mpfr_nextbelow (mpc_imagref(res)); 341 } 342 else if (tmpsgn < 0) 343 { 344 mpfr_nextbelow (mpc_imagref(res)); 345 isinf = mpfr_inf_p (mpc_imagref(res)); 346 mpfr_nextabove (mpc_imagref(res)); 347 } 348 if (isinf) 349 { 350 mpfr_set_inf (mpc_imagref(res), tmpsgn); 351 overflow_im = 1; 352 } 353 mpc_set (a, res, rnd); 354 goto end; 355 } 356 357 /* divide the product by the norm */ 358 if (inexact_norm == 0 && (inexact_re == 0 || inexact_im == 0)) { 359 /* The division has good chances to be exact in at least one part. */ 360 /* Since this can cause problems when not rounding to the nearest, */ 361 /* we use the division code of mpfr, which handles the situation. */ 362 mpfr_clear_underflow (); 363 mpfr_clear_overflow (); 364 inexact_re |= mpfr_div (mpc_realref (res), mpc_realref (res), q, GMP_RNDZ); 365 underflow_re = mpfr_underflow_p (); 366 overflow_re = mpfr_overflow_p (); 367 ok_re = !inexact_re || underflow_re || overflow_re 368 || mpfr_can_round (mpc_realref (res), prec - 4, GMP_RNDN, 369 GMP_RNDZ, MPC_PREC_RE(a) + (rnd_re == GMP_RNDN)); 370 371 if (ok_re) /* compute imaginary part */ { 372 mpfr_clear_underflow (); 373 mpfr_clear_overflow (); 374 inexact_im |= mpfr_div (mpc_imagref (res), mpc_imagref (res), q, GMP_RNDZ); 375 underflow_im = mpfr_underflow_p (); 376 overflow_im = mpfr_overflow_p (); 377 ok_im = !inexact_im || underflow_im || overflow_im 378 || mpfr_can_round (mpc_imagref (res), prec - 4, GMP_RNDN, 379 GMP_RNDZ, MPC_PREC_IM(a) + (rnd_im == GMP_RNDN)); 380 } 381 } 382 else { 383 /* The division is inexact, so for efficiency reasons we invert q */ 384 /* only once and multiply by the inverse. */ 385 if (mpfr_ui_div (q, 1ul, q, GMP_RNDZ) || inexact_norm) { 386 /* if 1/q is inexact, the approximations of the real and 387 imaginary part below will be inexact, unless RE(res) 388 or IM(res) is zero */ 389 inexact_re |= ~mpfr_zero_p (mpc_realref (res)); 390 inexact_im |= ~mpfr_zero_p (mpc_imagref (res)); 391 } 392 mpfr_clear_underflow (); 393 mpfr_clear_overflow (); 394 inexact_re |= mpfr_mul (mpc_realref (res), mpc_realref (res), q, GMP_RNDZ); 395 underflow_re = mpfr_underflow_p (); 396 overflow_re = mpfr_overflow_p (); 397 ok_re = !inexact_re || underflow_re || overflow_re 398 || mpfr_can_round (mpc_realref (res), prec - 4, GMP_RNDN, 399 GMP_RNDZ, MPC_PREC_RE(a) + (rnd_re == GMP_RNDN)); 400 401 if (ok_re) /* compute imaginary part */ { 402 mpfr_clear_underflow (); 403 mpfr_clear_overflow (); 404 inexact_im |= mpfr_mul (mpc_imagref (res), mpc_imagref (res), q, GMP_RNDZ); 405 underflow_im = mpfr_underflow_p (); 406 overflow_im = mpfr_overflow_p (); 407 ok_im = !inexact_im || underflow_im || overflow_im 408 || mpfr_can_round (mpc_imagref (res), prec - 4, GMP_RNDN, 409 GMP_RNDZ, MPC_PREC_IM(a) + (rnd_im == GMP_RNDN)); 410 } 411 } 412 } while ((!ok_re || !ok_im) && !underflow_norm && !overflow_norm 413 && !underflow_prod && !overflow_prod); 414 415 inex = mpc_set (a, res, rnd); 416 inexact_re = MPC_INEX_RE (inex); 417 inexact_im = MPC_INEX_IM (inex); 418 419 end: 420 /* fix values and inexact flags in case of overflow/underflow */ 421 /* FIXME: heuristic, certainly does not cover all cases */ 422 if (overflow_re || (underflow_norm && !underflow_prod)) { 423 mpfr_set_inf (mpc_realref (a), mpfr_sgn (mpc_realref (res))); 424 inexact_re = mpfr_sgn (mpc_realref (res)); 425 } 426 else if (underflow_re || (overflow_norm && !overflow_prod)) { 427 inexact_re = mpfr_signbit (mpc_realref (res)) ? 1 : -1; 428 mpfr_set_zero (mpc_realref (a), -inexact_re); 429 } 430 if (overflow_im || (underflow_norm && !underflow_prod)) { 431 mpfr_set_inf (mpc_imagref (a), mpfr_sgn (mpc_imagref (res))); 432 inexact_im = mpfr_sgn (mpc_imagref (res)); 433 } 434 else if (underflow_im || (overflow_norm && !overflow_prod)) { 435 inexact_im = mpfr_signbit (mpc_imagref (res)) ? 1 : -1; 436 mpfr_set_zero (mpc_imagref (a), -inexact_im); 437 } 438 439 mpc_clear (res); 440 mpfr_clear (q); 441 442 /* restore underflow and overflow flags from MPFR */ 443 if (saved_underflow) 444 mpfr_set_underflow (); 445 if (saved_overflow) 446 mpfr_set_overflow (); 447 448 return MPC_INEX (inexact_re, inexact_im); 449} 450