1 /* mpc_mul -- Multiply two complex numbers 2 3 Copyright (C) 2002, 2004, 2005, 2008, 2009, 2010, 2011, 2012 INRIA 4 5 This file is part of GNU MPC. 6 7 GNU MPC is free software; you can redistribute it and/or modify it under 8 the terms of the GNU Lesser General Public License as published by the 9 Free Software Foundation; either version 3 of the License, or (at your 10 option) any later version. 11 12 GNU MPC is distributed in the hope that it will be useful, but WITHOUT ANY 13 WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS 14 FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public License for 15 more details. 16 17 You should have received a copy of the GNU Lesser General Public License 18 along with this program. If not, see http://www.gnu.org/licenses/ . 19 */ 20 21 #include <stdio.h> /* for MPC_ASSERT */ 22 #include "mpc-impl.h" 23 24 #define mpz_add_si(z,x,y) do { \ 25 if (y >= 0) \ 26 mpz_add_ui (z, x, (long int) y); \ 27 else \ 28 mpz_sub_ui (z, x, (long int) (-y)); \ 29 } while (0); 30 31 /* compute z=x*y when x has an infinite part */ 32 static int 33 mul_infinite (mpc_ptr z, mpc_srcptr x, mpc_srcptr y) 34 { 35 /* Let x=xr+i*xi and y=yr+i*yi; extract the signs of the operands */ 36 int xrs = mpfr_signbit (mpc_realref (x)) ? -1 : 1; 37 int xis = mpfr_signbit (mpc_imagref (x)) ? -1 : 1; 38 int yrs = mpfr_signbit (mpc_realref (y)) ? -1 : 1; 39 int yis = mpfr_signbit (mpc_imagref (y)) ? -1 : 1; 40 41 int u, v; 42 43 /* compute the sign of 44 u = xrs * yrs * xr * yr - xis * yis * xi * yi 45 v = xrs * yis * xr * yi + xis * yrs * xi * yr 46 +1 if positive, -1 if negatiye, 0 if NaN */ 47 if ( mpfr_nan_p (mpc_realref (x)) || mpfr_nan_p (mpc_imagref (x)) 48 || mpfr_nan_p (mpc_realref (y)) || mpfr_nan_p (mpc_imagref (y))) { 49 u = 0; 50 v = 0; 51 } 52 else if (mpfr_inf_p (mpc_realref (x))) { 53 /* x = (+/-inf) xr + i*xi */ 54 u = ( mpfr_zero_p (mpc_realref (y)) 55 || (mpfr_inf_p (mpc_imagref (x)) && mpfr_zero_p (mpc_imagref (y))) 56 || (mpfr_zero_p (mpc_imagref (x)) && mpfr_inf_p (mpc_imagref (y))) 57 || ( (mpfr_inf_p (mpc_imagref (x)) || mpfr_inf_p (mpc_imagref (y))) 58 && xrs*yrs == xis*yis) 59 ? 0 : xrs * yrs); 60 v = ( mpfr_zero_p (mpc_imagref (y)) 61 || (mpfr_inf_p (mpc_imagref (x)) && mpfr_zero_p (mpc_realref (y))) 62 || (mpfr_zero_p (mpc_imagref (x)) && mpfr_inf_p (mpc_realref (y))) 63 || ( (mpfr_inf_p (mpc_imagref (x)) || mpfr_inf_p (mpc_imagref (x))) 64 && xrs*yis != xis*yrs) 65 ? 0 : xrs * yis); 66 } 67 else { 68 /* x = xr + i*(+/-inf) with |xr| != inf */ 69 u = ( mpfr_zero_p (mpc_imagref (y)) 70 || (mpfr_zero_p (mpc_realref (x)) && mpfr_inf_p (mpc_realref (y))) 71 || (mpfr_inf_p (mpc_realref (y)) && xrs*yrs == xis*yis) 72 ? 0 : -xis * yis); 73 v = ( mpfr_zero_p (mpc_realref (y)) 74 || (mpfr_zero_p (mpc_realref (x)) && mpfr_inf_p (mpc_imagref (y))) 75 || (mpfr_inf_p (mpc_imagref (y)) && xrs*yis != xis*yrs) 76 ? 0 : xis * yrs); 77 } 78 79 if (u == 0 && v == 0) { 80 /* Naive result is NaN+i*NaN. Obtain an infinity using the algorithm 81 given in Annex G.5.1 of the ISO C99 standard */ 82 int xr = (mpfr_zero_p (mpc_realref (x)) || mpfr_nan_p (mpc_realref (x)) ? 0 83 : (mpfr_inf_p (mpc_realref (x)) ? 1 : 0)); 84 int xi = (mpfr_zero_p (mpc_imagref (x)) || mpfr_nan_p (mpc_imagref (x)) ? 0 85 : (mpfr_inf_p (mpc_imagref (x)) ? 1 : 0)); 86 int yr = (mpfr_zero_p (mpc_realref (y)) || mpfr_nan_p (mpc_realref (y)) ? 0 : 1); 87 int yi = (mpfr_zero_p (mpc_imagref (y)) || mpfr_nan_p (mpc_imagref (y)) ? 0 : 1); 88 if (mpc_inf_p (y)) { 89 yr = mpfr_inf_p (mpc_realref (y)) ? 1 : 0; 90 yi = mpfr_inf_p (mpc_imagref (y)) ? 1 : 0; 91 } 92 93 u = xrs * xr * yrs * yr - xis * xi * yis * yi; 94 v = xrs * xr * yis * yi + xis * xi * yrs * yr; 95 } 96 97 if (u == 0) 98 mpfr_set_nan (mpc_realref (z)); 99 else 100 mpfr_set_inf (mpc_realref (z), u); 101 102 if (v == 0) 103 mpfr_set_nan (mpc_imagref (z)); 104 else 105 mpfr_set_inf (mpc_imagref (z), v); 106 107 return MPC_INEX (0, 0); /* exact */ 108 } 109 110 111 /* compute z = x*y for Im(y) == 0 */ 112 static int 113 mul_real (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) 114 { 115 int xrs, xis, yrs, yis; 116 int inex; 117 118 /* save signs of operands */ 119 xrs = MPFR_SIGNBIT (mpc_realref (x)); 120 xis = MPFR_SIGNBIT (mpc_imagref (x)); 121 yrs = MPFR_SIGNBIT (mpc_realref (y)); 122 yis = MPFR_SIGNBIT (mpc_imagref (y)); 123 124 inex = mpc_mul_fr (z, x, mpc_realref (y), rnd); 125 /* Signs of zeroes may be wrong. Their correction does not change the 126 inexact flag. */ 127 if (mpfr_zero_p (mpc_realref (z))) 128 mpfr_setsign (mpc_realref (z), mpc_realref (z), MPC_RND_RE(rnd) == GMP_RNDD 129 || (xrs != yrs && xis == yis), GMP_RNDN); 130 if (mpfr_zero_p (mpc_imagref (z))) 131 mpfr_setsign (mpc_imagref (z), mpc_imagref (z), MPC_RND_IM (rnd) == GMP_RNDD 132 || (xrs != yis && xis != yrs), GMP_RNDN); 133 134 return inex; 135 } 136 137 138 /* compute z = x*y for Re(y) == 0, and Im(x) != 0 and Im(y) != 0 */ 139 static int 140 mul_imag (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) 141 { 142 int sign; 143 int inex_re, inex_im; 144 int overlap = z == x || z == y; 145 mpc_t rop; 146 147 if (overlap) 148 mpc_init3 (rop, MPC_PREC_RE (z), MPC_PREC_IM (z)); 149 else 150 rop [0] = z[0]; 151 152 sign = (MPFR_SIGNBIT (mpc_realref (y)) != MPFR_SIGNBIT (mpc_imagref (x))) 153 && (MPFR_SIGNBIT (mpc_imagref (y)) != MPFR_SIGNBIT (mpc_realref (x))); 154 155 inex_re = -mpfr_mul (mpc_realref (rop), mpc_imagref (x), mpc_imagref (y), 156 INV_RND (MPC_RND_RE (rnd))); 157 mpfr_neg (mpc_realref (rop), mpc_realref (rop), GMP_RNDN); /* exact */ 158 inex_im = mpfr_mul (mpc_imagref (rop), mpc_realref (x), mpc_imagref (y), 159 MPC_RND_IM (rnd)); 160 mpc_set (z, rop, MPC_RNDNN); 161 162 /* Sign of zeroes may be wrong (note that Re(z) cannot be zero) */ 163 if (mpfr_zero_p (mpc_imagref (z))) 164 mpfr_setsign (mpc_imagref (z), mpc_imagref (z), MPC_RND_IM (rnd) == GMP_RNDD 165 || sign, GMP_RNDN); 166 167 if (overlap) 168 mpc_clear (rop); 169 170 return MPC_INEX (inex_re, inex_im); 171 } 172 173 174 static int 175 mpfr_fmma (mpfr_ptr z, mpfr_srcptr a, mpfr_srcptr b, mpfr_srcptr c, 176 mpfr_srcptr d, int sign, mpfr_rnd_t rnd) 177 { 178 /* Computes z = ab+cd if sign >= 0, or z = ab-cd if sign < 0. 179 Assumes that a, b, c, d are finite and non-zero; so any multiplication 180 of two of them yielding an infinity is an overflow, and a 181 multiplication yielding 0 is an underflow. 182 Assumes further that z is distinct from a, b, c, d. */ 183 184 int inex; 185 mpfr_t u, v; 186 187 /* u=a*b, v=sign*c*d exactly */ 188 mpfr_init2 (u, mpfr_get_prec (a) + mpfr_get_prec (b)); 189 mpfr_init2 (v, mpfr_get_prec (c) + mpfr_get_prec (d)); 190 mpfr_mul (u, a, b, GMP_RNDN); 191 mpfr_mul (v, c, d, GMP_RNDN); 192 if (sign < 0) 193 mpfr_neg (v, v, GMP_RNDN); 194 195 /* tentatively compute z as u+v; here we need z to be distinct 196 from a, b, c, d to not lose the latter */ 197 inex = mpfr_add (z, u, v, rnd); 198 199 if (mpfr_inf_p (z)) { 200 /* replace by "correctly rounded overflow" */ 201 mpfr_set_si (z, (mpfr_signbit (z) ? -1 : 1), GMP_RNDN); 202 inex = mpfr_mul_2ui (z, z, mpfr_get_emax (), rnd); 203 } 204 else if (mpfr_zero_p (u) && !mpfr_zero_p (v)) { 205 /* exactly u underflowed, determine inexact flag */ 206 inex = (mpfr_signbit (u) ? 1 : -1); 207 } 208 else if (mpfr_zero_p (v) && !mpfr_zero_p (u)) { 209 /* exactly v underflowed, determine inexact flag */ 210 inex = (mpfr_signbit (v) ? 1 : -1); 211 } 212 else if (mpfr_nan_p (z) || (mpfr_zero_p (u) && mpfr_zero_p (v))) { 213 /* In the first case, u and v are infinities with opposite signs. 214 In the second case, u and v are zeroes; their sum may be 0 or the 215 least representable number, with a sign to be determined. 216 Redo the computations with mpz_t exponents */ 217 mpfr_exp_t ea, eb, ec, ed; 218 mpz_t eu, ev; 219 /* cheat to work around the const qualifiers */ 220 221 /* Normalise the input by shifting and keep track of the shifts in 222 the exponents of u and v */ 223 ea = mpfr_get_exp (a); 224 eb = mpfr_get_exp (b); 225 ec = mpfr_get_exp (c); 226 ed = mpfr_get_exp (d); 227 228 mpfr_set_exp ((mpfr_ptr) a, (mpfr_prec_t) 0); 229 mpfr_set_exp ((mpfr_ptr) b, (mpfr_prec_t) 0); 230 mpfr_set_exp ((mpfr_ptr) c, (mpfr_prec_t) 0); 231 mpfr_set_exp ((mpfr_ptr) d, (mpfr_prec_t) 0); 232 233 mpz_init (eu); 234 mpz_init (ev); 235 mpz_set_si (eu, (long int) ea); 236 mpz_add_si (eu, eu, (long int) eb); 237 mpz_set_si (ev, (long int) ec); 238 mpz_add_si (ev, ev, (long int) ed); 239 240 /* recompute u and v and move exponents to eu and ev */ 241 mpfr_mul (u, a, b, GMP_RNDN); 242 /* exponent of u is non-positive */ 243 mpz_sub_ui (eu, eu, (unsigned long int) (-mpfr_get_exp (u))); 244 mpfr_set_exp (u, (mpfr_prec_t) 0); 245 mpfr_mul (v, c, d, GMP_RNDN); 246 if (sign < 0) 247 mpfr_neg (v, v, GMP_RNDN); 248 mpz_sub_ui (ev, ev, (unsigned long int) (-mpfr_get_exp (v))); 249 mpfr_set_exp (v, (mpfr_prec_t) 0); 250 251 if (mpfr_nan_p (z)) { 252 mpfr_exp_t emax = mpfr_get_emax (); 253 int overflow; 254 /* We have a = ma * 2^ea with 1/2 <= |ma| < 1 and ea <= emax, and 255 analogously for b. So eu <= 2*emax, and eu > emax since we have 256 an overflow. The same holds for ev. Shift u and v by as much as 257 possible so that one of them has exponent emax and the 258 remaining exponents in eu and ev are the same. Then carry out 259 the addition. Shifting u and v prevents an underflow. */ 260 if (mpz_cmp (eu, ev) >= 0) { 261 mpfr_set_exp (u, emax); 262 mpz_sub_ui (eu, eu, (long int) emax); 263 mpz_sub (ev, ev, eu); 264 mpfr_set_exp (v, (mpfr_exp_t) mpz_get_ui (ev)); 265 /* remaining common exponent is now in eu */ 266 } 267 else { 268 mpfr_set_exp (v, emax); 269 mpz_sub_ui (ev, ev, (long int) emax); 270 mpz_sub (eu, eu, ev); 271 mpfr_set_exp (u, (mpfr_exp_t) mpz_get_ui (eu)); 272 mpz_set (eu, ev); 273 /* remaining common exponent is now also in eu */ 274 } 275 inex = mpfr_add (z, u, v, rnd); 276 /* Result is finite since u and v have different signs. */ 277 overflow = mpfr_mul_2ui (z, z, mpz_get_ui (eu), rnd); 278 if (overflow) 279 inex = overflow; 280 } 281 else { 282 int underflow; 283 /* Addition of two zeroes with same sign. We have a = ma * 2^ea 284 with 1/2 <= |ma| < 1 and ea >= emin and similarly for b. 285 So 2*emin < 2*emin+1 <= eu < emin < 0, and analogously for v. */ 286 mpfr_exp_t emin = mpfr_get_emin (); 287 if (mpz_cmp (eu, ev) <= 0) { 288 mpfr_set_exp (u, emin); 289 mpz_add_ui (eu, eu, (unsigned long int) (-emin)); 290 mpz_sub (ev, ev, eu); 291 mpfr_set_exp (v, (mpfr_exp_t) mpz_get_si (ev)); 292 } 293 else { 294 mpfr_set_exp (v, emin); 295 mpz_add_ui (ev, ev, (unsigned long int) (-emin)); 296 mpz_sub (eu, eu, ev); 297 mpfr_set_exp (u, (mpfr_exp_t) mpz_get_si (eu)); 298 mpz_set (eu, ev); 299 } 300 inex = mpfr_add (z, u, v, rnd); 301 mpz_neg (eu, eu); 302 underflow = mpfr_div_2ui (z, z, mpz_get_ui (eu), rnd); 303 if (underflow) 304 inex = underflow; 305 } 306 307 mpz_clear (eu); 308 mpz_clear (ev); 309 310 mpfr_set_exp ((mpfr_ptr) a, ea); 311 mpfr_set_exp ((mpfr_ptr) b, eb); 312 mpfr_set_exp ((mpfr_ptr) c, ec); 313 mpfr_set_exp ((mpfr_ptr) d, ed); 314 /* works also when some of a, b, c, d are not all distinct */ 315 } 316 317 mpfr_clear (u); 318 mpfr_clear (v); 319 320 return inex; 321 } 322 323 324 int 325 mpc_mul_naive (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) 326 { 327 /* computes z=x*y by the schoolbook method, where x and y are assumed 328 to be finite and without zero parts */ 329 int overlap, inex; 330 mpc_t rop; 331 332 MPC_ASSERT ( mpfr_regular_p (mpc_realref (x)) && mpfr_regular_p (mpc_imagref (x)) 333 && mpfr_regular_p (mpc_realref (y)) && mpfr_regular_p (mpc_imagref (y))); 334 overlap = (z == x) || (z == y); 335 if (overlap) 336 mpc_init3 (rop, MPC_PREC_RE (z), MPC_PREC_IM (z)); 337 else 338 rop [0] = z [0]; 339 340 inex = MPC_INEX (mpfr_fmma (mpc_realref (rop), mpc_realref (x), mpc_realref (y), mpc_imagref (x), 341 mpc_imagref (y), -1, MPC_RND_RE (rnd)), 342 mpfr_fmma (mpc_imagref (rop), mpc_realref (x), mpc_imagref (y), mpc_imagref (x), 343 mpc_realref (y), +1, MPC_RND_IM (rnd))); 344 345 mpc_set (z, rop, MPC_RNDNN); 346 if (overlap) 347 mpc_clear (rop); 348 349 return inex; 350 } 351 352 353 int 354 mpc_mul_karatsuba (mpc_ptr rop, mpc_srcptr op1, mpc_srcptr op2, mpc_rnd_t rnd) 355 { 356 /* computes rop=op1*op2 by a Karatsuba algorithm, where op1 and op2 357 are assumed to be finite and without zero parts */ 358 mpfr_srcptr a, b, c, d; 359 int mul_i, ok, inexact, mul_a, mul_c, inex_re = 0, inex_im = 0, sign_x, sign_u; 360 mpfr_t u, v, w, x; 361 mpfr_prec_t prec, prec_re, prec_u, prec_v, prec_w; 362 mpfr_rnd_t rnd_re, rnd_u; 363 int overlap; 364 /* true if rop == op1 or rop == op2 */ 365 mpc_t result; 366 /* overlap is quite difficult to handle, because we have to tentatively 367 round the variable u in the end to either the real or the imaginary 368 part of rop (it is not possible to tell now whether the real or 369 imaginary part is used). If this fails, we have to start again and 370 need the correct values of op1 and op2. 371 So we just create a new variable for the result in this case. */ 372 int loop; 373 const int MAX_MUL_LOOP = 1; 374 375 overlap = (rop == op1) || (rop == op2); 376 if (overlap) 377 mpc_init3 (result, MPC_PREC_RE (rop), MPC_PREC_IM (rop)); 378 else 379 result [0] = rop [0]; 380 381 a = mpc_realref(op1); 382 b = mpc_imagref(op1); 383 c = mpc_realref(op2); 384 d = mpc_imagref(op2); 385 386 /* (a + i*b) * (c + i*d) = [ac - bd] + i*[ad + bc] */ 387 388 mul_i = 0; /* number of multiplications by i */ 389 mul_a = 1; /* implicit factor for a */ 390 mul_c = 1; /* implicit factor for c */ 391 392 if (mpfr_cmp_abs (a, b) < 0) 393 { 394 MPFR_SWAP (a, b); 395 mul_i ++; 396 mul_a = -1; /* consider i * (a+i*b) = -b + i*a */ 397 } 398 399 if (mpfr_cmp_abs (c, d) < 0) 400 { 401 MPFR_SWAP (c, d); 402 mul_i ++; 403 mul_c = -1; /* consider -d + i*c instead of c + i*d */ 404 } 405 406 /* find the precision and rounding mode for the new real part */ 407 if (mul_i % 2) 408 { 409 prec_re = MPC_PREC_IM(rop); 410 rnd_re = MPC_RND_IM(rnd); 411 } 412 else /* mul_i = 0 or 2 */ 413 { 414 prec_re = MPC_PREC_RE(rop); 415 rnd_re = MPC_RND_RE(rnd); 416 } 417 418 if (mul_i) 419 rnd_re = INV_RND(rnd_re); 420 421 /* now |a| >= |b| and |c| >= |d| */ 422 prec = MPC_MAX_PREC(rop); 423 424 mpfr_init2 (v, prec_v = mpfr_get_prec (a) + mpfr_get_prec (d)); 425 mpfr_init2 (w, prec_w = mpfr_get_prec (b) + mpfr_get_prec (c)); 426 mpfr_init2 (u, 2); 427 mpfr_init2 (x, 2); 428 429 inexact = mpfr_mul (v, a, d, GMP_RNDN); 430 if (inexact) { 431 /* over- or underflow */ 432 ok = 0; 433 goto clear; 434 } 435 if (mul_a == -1) 436 mpfr_neg (v, v, GMP_RNDN); 437 438 inexact = mpfr_mul (w, b, c, GMP_RNDN); 439 if (inexact) { 440 /* over- or underflow */ 441 ok = 0; 442 goto clear; 443 } 444 if (mul_c == -1) 445 mpfr_neg (w, w, GMP_RNDN); 446 447 /* compute sign(v-w) */ 448 sign_x = mpfr_cmp_abs (v, w); 449 if (sign_x > 0) 450 sign_x = 2 * mpfr_sgn (v) - mpfr_sgn (w); 451 else if (sign_x == 0) 452 sign_x = mpfr_sgn (v) - mpfr_sgn (w); 453 else 454 sign_x = mpfr_sgn (v) - 2 * mpfr_sgn (w); 455 456 sign_u = mul_a * mpfr_sgn (a) * mul_c * mpfr_sgn (c); 457 458 if (sign_x * sign_u < 0) 459 { /* swap inputs */ 460 MPFR_SWAP (a, c); 461 MPFR_SWAP (b, d); 462 mpfr_swap (v, w); 463 { int tmp; tmp = mul_a; mul_a = mul_c; mul_c = tmp; } 464 sign_x = - sign_x; 465 } 466 467 /* now sign_x * sign_u >= 0 */ 468 loop = 0; 469 do 470 { 471 loop++; 472 /* the following should give failures with prob. <= 1/prec */ 473 prec += mpc_ceil_log2 (prec) + 3; 474 475 mpfr_set_prec (u, prec_u = prec); 476 mpfr_set_prec (x, prec); 477 478 /* first compute away(b +/- a) and store it in u */ 479 inexact = (mul_a == -1 ? 480 ROUND_AWAY (mpfr_sub (u, b, a, MPFR_RNDA), u) : 481 ROUND_AWAY (mpfr_add (u, b, a, MPFR_RNDA), u)); 482 483 /* then compute away(+/-c - d) and store it in x */ 484 inexact |= (mul_c == -1 ? 485 ROUND_AWAY (mpfr_add (x, c, d, MPFR_RNDA), x) : 486 ROUND_AWAY (mpfr_sub (x, c, d, MPFR_RNDA), x)); 487 if (mul_c == -1) 488 mpfr_neg (x, x, GMP_RNDN); 489 490 if (inexact == 0) 491 mpfr_prec_round (u, prec_u = 2 * prec, GMP_RNDN); 492 493 /* compute away(u*x) and store it in u */ 494 inexact |= ROUND_AWAY (mpfr_mul (u, u, x, MPFR_RNDA), u); 495 /* (a+b)*(c-d) */ 496 497 /* if all computations are exact up to here, it may be that 498 the real part is exact, thus we need if possible to 499 compute v - w exactly */ 500 if (inexact == 0) 501 { 502 mpfr_prec_t prec_x; 503 /* v and w are different from 0, so mpfr_get_exp is safe to use */ 504 prec_x = SAFE_ABS (mpfr_exp_t, mpfr_get_exp (v) - mpfr_get_exp (w)) 505 + MPC_MAX (prec_v, prec_w) + 1; 506 /* +1 is necessary for a potential carry */ 507 /* ensure we do not use a too large precision */ 508 if (prec_x > prec_u) 509 prec_x = prec_u; 510 if (prec_x > prec) 511 mpfr_prec_round (x, prec_x, GMP_RNDN); 512 } 513 514 rnd_u = (sign_u > 0) ? GMP_RNDU : GMP_RNDD; 515 inexact |= mpfr_sub (x, v, w, rnd_u); /* ad - bc */ 516 517 /* in case u=0, ensure that rnd_u rounds x away from zero */ 518 if (mpfr_sgn (u) == 0) 519 rnd_u = (mpfr_sgn (x) > 0) ? GMP_RNDU : GMP_RNDD; 520 inexact |= mpfr_add (u, u, x, rnd_u); /* ac - bd */ 521 522 ok = inexact == 0 || 523 mpfr_can_round (u, prec_u - 3, rnd_u, GMP_RNDZ, 524 prec_re + (rnd_re == GMP_RNDN)); 525 /* this ensures both we can round correctly and determine the correct 526 inexact flag (for rounding to nearest) */ 527 } 528 while (!ok && loop <= MAX_MUL_LOOP); 529 /* after MAX_MUL_LOOP rounds, use mpc_naive instead */ 530 531 if (ok) { 532 /* if inexact is zero, then u is exactly ac-bd, otherwise fix the sign 533 of the inexact flag for u, which was rounded away from ac-bd */ 534 if (inexact != 0) 535 inexact = mpfr_sgn (u); 536 537 if (mul_i == 0) 538 { 539 inex_re = mpfr_set (mpc_realref(result), u, MPC_RND_RE(rnd)); 540 if (inex_re == 0) 541 { 542 inex_re = inexact; /* u is rounded away from 0 */ 543 inex_im = mpfr_add (mpc_imagref(result), v, w, MPC_RND_IM(rnd)); 544 } 545 else 546 inex_im = mpfr_add (mpc_imagref(result), v, w, MPC_RND_IM(rnd)); 547 } 548 else if (mul_i == 1) /* (x+i*y)/i = y - i*x */ 549 { 550 inex_im = mpfr_neg (mpc_imagref(result), u, MPC_RND_IM(rnd)); 551 if (inex_im == 0) 552 { 553 inex_im = -inexact; /* u is rounded away from 0 */ 554 inex_re = mpfr_add (mpc_realref(result), v, w, MPC_RND_RE(rnd)); 555 } 556 else 557 inex_re = mpfr_add (mpc_realref(result), v, w, MPC_RND_RE(rnd)); 558 } 559 else /* mul_i = 2, z/i^2 = -z */ 560 { 561 inex_re = mpfr_neg (mpc_realref(result), u, MPC_RND_RE(rnd)); 562 if (inex_re == 0) 563 { 564 inex_re = -inexact; /* u is rounded away from 0 */ 565 inex_im = -mpfr_add (mpc_imagref(result), v, w, 566 INV_RND(MPC_RND_IM(rnd))); 567 mpfr_neg (mpc_imagref(result), mpc_imagref(result), MPC_RND_IM(rnd)); 568 } 569 else 570 { 571 inex_im = -mpfr_add (mpc_imagref(result), v, w, 572 INV_RND(MPC_RND_IM(rnd))); 573 mpfr_neg (mpc_imagref(result), mpc_imagref(result), MPC_RND_IM(rnd)); 574 } 575 } 576 577 mpc_set (rop, result, MPC_RNDNN); 578 } 579 580 clear: 581 mpfr_clear (u); 582 mpfr_clear (v); 583 mpfr_clear (w); 584 mpfr_clear (x); 585 if (overlap) 586 mpc_clear (result); 587 588 if (ok) 589 return MPC_INEX(inex_re, inex_im); 590 else 591 return mpc_mul_naive (rop, op1, op2, rnd); 592 } 593 594 595 int 596 mpc_mul (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd) 597 { 598 /* Conforming to ISO C99 standard (G.5.1 multiplicative operators), 599 infinities are treated specially if both parts are NaN when computed 600 naively. */ 601 if (mpc_inf_p (b)) 602 return mul_infinite (a, b, c); 603 if (mpc_inf_p (c)) 604 return mul_infinite (a, c, b); 605 606 /* NaN contamination of both parts in result */ 607 if (mpfr_nan_p (mpc_realref (b)) || mpfr_nan_p (mpc_imagref (b)) 608 || mpfr_nan_p (mpc_realref (c)) || mpfr_nan_p (mpc_imagref (c))) { 609 mpfr_set_nan (mpc_realref (a)); 610 mpfr_set_nan (mpc_imagref (a)); 611 return MPC_INEX (0, 0); 612 } 613 614 /* check for real multiplication */ 615 if (mpfr_zero_p (mpc_imagref (b))) 616 return mul_real (a, c, b, rnd); 617 if (mpfr_zero_p (mpc_imagref (c))) 618 return mul_real (a, b, c, rnd); 619 620 /* check for purely imaginary multiplication */ 621 if (mpfr_zero_p (mpc_realref (b))) 622 return mul_imag (a, c, b, rnd); 623 if (mpfr_zero_p (mpc_realref (c))) 624 return mul_imag (a, b, c, rnd); 625 626 /* If the real and imaginary part of one argument have a very different */ 627 /* exponent, it is not reasonable to use Karatsuba multiplication. */ 628 if ( SAFE_ABS (mpfr_exp_t, 629 mpfr_get_exp (mpc_realref (b)) - mpfr_get_exp (mpc_imagref (b))) 630 > (mpfr_exp_t) MPC_MAX_PREC (b) / 2 631 || SAFE_ABS (mpfr_exp_t, 632 mpfr_get_exp (mpc_realref (c)) - mpfr_get_exp (mpc_imagref (c))) 633 > (mpfr_exp_t) MPC_MAX_PREC (c) / 2) 634 return mpc_mul_naive (a, b, c, rnd); 635 else 636 return ((MPC_MAX_PREC(a) 637 <= (mpfr_prec_t) MUL_KARATSUBA_THRESHOLD * BITS_PER_MP_LIMB) 638 ? mpc_mul_naive : mpc_mul_karatsuba) (a, b, c, rnd); 639 } 640