1/* mpc_mul -- Multiply two complex numbers 2 3Copyright (C) 2002, 2004, 2005, 2008, 2009, 2010, 2011, 2012, 2016, 2020, 2022 INRIA 4 5This file is part of GNU MPC. 6 7GNU MPC is free software; you can redistribute it and/or modify it under 8 he 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 <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 */ 32static int 33mul_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 negative, 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 */ 112static int 113mul_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) == MPFR_RNDD 129 || (xrs != yrs && xis == yis), MPFR_RNDN); 130 if (mpfr_zero_p (mpc_imagref (z))) 131 mpfr_setsign (mpc_imagref (z), mpc_imagref (z), MPC_RND_IM (rnd) == MPFR_RNDD 132 || (xrs != yis && xis != yrs), MPFR_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 */ 139static int 140mul_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), MPFR_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) == MPFR_RNDD 165 || sign, MPFR_RNDN); 166 167 if (overlap) 168 mpc_clear (rop); 169 170 return MPC_INEX (inex_re, inex_im); 171} 172 173 174int 175mpc_mul_naive (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) 176{ 177 /* computes z=x*y by the schoolbook method, where x and y are assumed 178 to be finite and without zero parts */ 179 int overlap, inex_re, inex_im; 180 mpc_t rop; 181 182 MPC_ASSERT ( mpfr_regular_p (mpc_realref (x)) && mpfr_regular_p (mpc_imagref (x)) 183 && mpfr_regular_p (mpc_realref (y)) && mpfr_regular_p (mpc_imagref (y))); 184 overlap = (z == x) || (z == y); 185 if (overlap) 186 mpc_init3 (rop, MPC_PREC_RE (z), MPC_PREC_IM (z)); 187 else 188 rop [0] = z [0]; 189 190 inex_re = mpfr_fmms (mpc_realref (rop), mpc_realref (x), mpc_realref (y), 191 mpc_imagref (x), mpc_imagref (y), MPC_RND_RE (rnd)); 192 inex_im = mpfr_fmma (mpc_imagref (rop), mpc_realref (x), mpc_imagref (y), 193 mpc_imagref (x), mpc_realref (y), MPC_RND_IM (rnd)); 194 195 mpc_set (z, rop, MPC_RNDNN); 196 if (overlap) 197 mpc_clear (rop); 198 199 return MPC_INEX (inex_re, inex_im); 200} 201 202 203int 204mpc_mul_karatsuba (mpc_ptr rop, mpc_srcptr op1, mpc_srcptr op2, mpc_rnd_t rnd) 205{ 206 /* computes rop=op1*op2 by a Karatsuba algorithm, where op1 and op2 207 are assumed to be finite and without zero parts */ 208 mpfr_srcptr a, b, c, d; 209 int mul_i, ok, inexact, mul_a, mul_c, inex_re = 0, inex_im = 0, sign_x, sign_u; 210 mpfr_t u, v, w, x; 211 mpfr_prec_t prec, prec_re, prec_u, prec_v, prec_w; 212 mpfr_rnd_t rnd_re, rnd_u; 213 int overlap; 214 /* true if rop == op1 or rop == op2 */ 215 mpc_t result; 216 /* overlap is quite difficult to handle, because we have to tentatively 217 round the variable u in the end to either the real or the imaginary 218 part of rop (it is not possible to tell now whether the real or 219 imaginary part is used). If this fails, we have to start again and 220 need the correct values of op1 and op2. 221 So we just create a new variable for the result in this case. */ 222 mpfr_ptr ref; 223 int loop; 224 const int MAX_MUL_LOOP = 1; 225 226 overlap = (rop == op1) || (rop == op2); 227 if (overlap) 228 mpc_init3 (result, MPC_PREC_RE (rop), MPC_PREC_IM (rop)); 229 else 230 result [0] = rop [0]; 231 232 a = mpc_realref(op1); 233 b = mpc_imagref(op1); 234 c = mpc_realref(op2); 235 d = mpc_imagref(op2); 236 237 /* (a + i*b) * (c + i*d) = [ac - bd] + i*[ad + bc] */ 238 239 mul_i = 0; /* number of multiplications by i */ 240 mul_a = 1; /* implicit factor for a */ 241 mul_c = 1; /* implicit factor for c */ 242 243 if (mpfr_cmp_abs (a, b) < 0) 244 { 245 MPFR_SWAP (a, b); 246 mul_i ++; 247 mul_a = -1; /* consider i * (a+i*b) = -b + i*a */ 248 } 249 250 if (mpfr_cmp_abs (c, d) < 0) 251 { 252 MPFR_SWAP (c, d); 253 mul_i ++; 254 mul_c = -1; /* consider -d + i*c instead of c + i*d */ 255 } 256 257 /* find the precision and rounding mode for the new real part */ 258 if (mul_i % 2) 259 { 260 prec_re = MPC_PREC_IM(rop); 261 rnd_re = MPC_RND_IM(rnd); 262 } 263 else /* mul_i = 0 or 2 */ 264 { 265 prec_re = MPC_PREC_RE(rop); 266 rnd_re = MPC_RND_RE(rnd); 267 } 268 269 if (mul_i) 270 rnd_re = INV_RND(rnd_re); 271 272 /* now |a| >= |b| and |c| >= |d| */ 273 prec = MPC_MAX_PREC(rop); 274 275 mpfr_init2 (v, prec_v = mpfr_get_prec (a) + mpfr_get_prec (d)); 276 mpfr_init2 (w, prec_w = mpfr_get_prec (b) + mpfr_get_prec (c)); 277 mpfr_init2 (u, 2); 278 mpfr_init2 (x, 2); 279 280 inexact = mpfr_mul (v, a, d, MPFR_RNDN); 281 if (inexact) { 282 /* over- or underflow */ 283 ok = 0; 284 goto clear; 285 } 286 if (mul_a == -1) 287 mpfr_neg (v, v, MPFR_RNDN); 288 289 inexact = mpfr_mul (w, b, c, MPFR_RNDN); 290 if (inexact) { 291 /* over- or underflow */ 292 ok = 0; 293 goto clear; 294 } 295 if (mul_c == -1) 296 mpfr_neg (w, w, MPFR_RNDN); 297 298 /* compute sign(v-w) */ 299 sign_x = mpfr_cmp_abs (v, w); 300 if (sign_x > 0) 301 sign_x = 2 * mpfr_sgn (v) - mpfr_sgn (w); 302 else if (sign_x == 0) 303 sign_x = mpfr_sgn (v) - mpfr_sgn (w); 304 else 305 sign_x = mpfr_sgn (v) - 2 * mpfr_sgn (w); 306 307 sign_u = mul_a * mpfr_sgn (a) * mul_c * mpfr_sgn (c); 308 309 if (sign_x * sign_u < 0) 310 { /* swap inputs */ 311 MPFR_SWAP (a, c); 312 MPFR_SWAP (b, d); 313 mpfr_swap (v, w); 314 { int tmp; tmp = mul_a; mul_a = mul_c; mul_c = tmp; } 315 sign_x = - sign_x; 316 } 317 318 /* now sign_x * sign_u >= 0 */ 319 loop = 0; 320 do 321 { 322 loop++; 323 /* the following should give failures with prob. <= 1/prec */ 324 prec += mpc_ceil_log2 (prec) + 3; 325 326 mpfr_set_prec (u, prec_u = prec); 327 mpfr_set_prec (x, prec); 328 329 /* first compute away(b +/- a) and store it in u */ 330 inexact = (mul_a == -1 ? 331 mpfr_sub (u, b, a, MPFR_RNDA) : 332 mpfr_add (u, b, a, MPFR_RNDA)); 333 334 /* then compute away(+/-c - d) and store it in x */ 335 inexact |= (mul_c == -1 ? 336 mpfr_add (x, c, d, MPFR_RNDA) : 337 mpfr_sub (x, c, d, MPFR_RNDA)); 338 if (mul_c == -1) 339 mpfr_neg (x, x, MPFR_RNDN); 340 341 if (inexact == 0) 342 mpfr_prec_round (u, prec_u = 2 * prec, MPFR_RNDN); 343 344 /* compute away(u*x) and store it in u */ 345 inexact |= mpfr_mul (u, u, x, MPFR_RNDA); 346 /* (a+b)*(c-d) */ 347 348 /* if all computations are exact up to here, it may be that 349 the real part is exact, thus we need if possible to 350 compute v - w exactly */ 351 if (inexact == 0) 352 { 353 mpfr_prec_t prec_x; 354 /* v and w are different from 0, so mpfr_get_exp is safe to use */ 355 prec_x = SAFE_ABS (mpfr_exp_t, mpfr_get_exp (v) - mpfr_get_exp (w)) 356 + MPC_MAX (prec_v, prec_w) + 1; 357 /* +1 is necessary for a potential carry */ 358 /* ensure we do not use a too large precision */ 359 if (prec_x > prec_u) 360 prec_x = prec_u; 361 if (prec_x > prec) 362 mpfr_prec_round (x, prec_x, MPFR_RNDN); 363 } 364 365 rnd_u = (sign_u > 0) ? MPFR_RNDU : MPFR_RNDD; 366 inexact |= mpfr_sub (x, v, w, rnd_u); /* ad - bc */ 367 368 /* in case u=0, ensure that rnd_u rounds x away from zero */ 369 if (mpfr_sgn (u) == 0) 370 rnd_u = (mpfr_sgn (x) > 0) ? MPFR_RNDU : MPFR_RNDD; 371 inexact |= mpfr_add (u, u, x, rnd_u); /* ac - bd */ 372 373 ok = inexact == 0 || 374 mpfr_can_round (u, prec_u - 3, rnd_u, MPFR_RNDZ, 375 prec_re + (rnd_re == MPFR_RNDN)); 376 /* this ensures both we can round correctly and determine the correct 377 inexact flag (for rounding to nearest) */ 378 } 379 while (!ok && loop <= MAX_MUL_LOOP); 380 /* after MAX_MUL_LOOP rounds, use mpc_naive instead */ 381 382 if (ok) { 383 /* if inexact is zero, then u is exactly ac-bd, otherwise fix the sign 384 of the inexact flag for u, which was rounded away from ac-bd */ 385 if (inexact != 0) 386 inexact = mpfr_sgn (u); 387 388 if (mul_i == 0) 389 { 390 inex_re = mpfr_set (mpc_realref(result), u, MPC_RND_RE(rnd)); 391 if (inex_re == 0) 392 { 393 inex_re = inexact; /* u is rounded away from 0 */ 394 inex_im = mpfr_add (mpc_imagref(result), v, w, MPC_RND_IM(rnd)); 395 } 396 else 397 inex_im = mpfr_add (mpc_imagref(result), v, w, MPC_RND_IM(rnd)); 398 } 399 else if (mul_i == 1) /* (x+i*y)/i = y - i*x */ 400 { 401 inex_im = mpfr_neg (mpc_imagref(result), u, MPC_RND_IM(rnd)); 402 if (inex_im == 0) 403 { 404 inex_im = -inexact; /* u is rounded away from 0 */ 405 inex_re = mpfr_add (mpc_realref(result), v, w, MPC_RND_RE(rnd)); 406 } 407 else 408 inex_re = mpfr_add (mpc_realref(result), v, w, MPC_RND_RE(rnd)); 409 } 410 else /* mul_i = 2, z/i^2 = -z */ 411 { 412 inex_re = mpfr_neg (mpc_realref(result), u, MPC_RND_RE(rnd)); 413 if (inex_re == 0) 414 { 415 inex_re = -inexact; /* u is rounded away from 0 */ 416 inex_im = -mpfr_add (mpc_imagref(result), v, w, 417 INV_RND(MPC_RND_IM(rnd))); 418 mpfr_neg (mpc_imagref(result), mpc_imagref(result), MPC_RND_IM(rnd)); 419 } 420 else 421 { 422 inex_im = -mpfr_add (mpc_imagref(result), v, w, 423 INV_RND(MPC_RND_IM(rnd))); 424 mpfr_neg (mpc_imagref(result), mpc_imagref(result), MPC_RND_IM(rnd)); 425 } 426 } 427 428 /* Clear potential signs of 0. */ 429 if (!inex_re) { 430 ref = mpc_realref (result); 431 if (mpfr_zero_p (ref) && mpfr_signbit (ref)) 432 MPFR_CHANGE_SIGN (ref); 433 } 434 if (!inex_im) { 435 ref = mpc_imagref (result); 436 if (mpfr_zero_p (ref) && mpfr_signbit (ref)) 437 MPFR_CHANGE_SIGN (ref); 438 } 439 440 mpc_set (rop, result, MPC_RNDNN); 441 } 442 443clear: 444 mpfr_clear (u); 445 mpfr_clear (v); 446 mpfr_clear (w); 447 mpfr_clear (x); 448 if (overlap) 449 mpc_clear (result); 450 451 if (ok) 452 return MPC_INEX(inex_re, inex_im); 453 else 454 return mpc_mul_naive (rop, op1, op2, rnd); 455} 456 457 458int 459mpc_mul (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd) 460{ 461 /* Conforming to ISO C99 standard (G.5.1 multiplicative operators), 462 infinities are treated specially if both parts are NaN when computed 463 naively. */ 464 if (mpc_inf_p (b)) 465 return mul_infinite (a, b, c); 466 if (mpc_inf_p (c)) 467 return mul_infinite (a, c, b); 468 469 /* NaN contamination of both parts in result */ 470 if (mpfr_nan_p (mpc_realref (b)) || mpfr_nan_p (mpc_imagref (b)) 471 || mpfr_nan_p (mpc_realref (c)) || mpfr_nan_p (mpc_imagref (c))) { 472 mpfr_set_nan (mpc_realref (a)); 473 mpfr_set_nan (mpc_imagref (a)); 474 return MPC_INEX (0, 0); 475 } 476 477 /* check for real multiplication */ 478 if (mpfr_zero_p (mpc_imagref (b))) 479 return mul_real (a, c, b, rnd); 480 if (mpfr_zero_p (mpc_imagref (c))) 481 return mul_real (a, b, c, rnd); 482 483 /* check for purely imaginary multiplication */ 484 if (mpfr_zero_p (mpc_realref (b))) 485 return mul_imag (a, c, b, rnd); 486 if (mpfr_zero_p (mpc_realref (c))) 487 return mul_imag (a, b, c, rnd); 488 489 /* If the real and imaginary part of one argument have a very different */ 490 /* exponent, it is not reasonable to use Karatsuba multiplication. */ 491 if ( SAFE_ABS (mpfr_exp_t, 492 mpfr_get_exp (mpc_realref (b)) - mpfr_get_exp (mpc_imagref (b))) 493 > (mpfr_exp_t) MPC_MAX_PREC (b) / 2 494 || SAFE_ABS (mpfr_exp_t, 495 mpfr_get_exp (mpc_realref (c)) - mpfr_get_exp (mpc_imagref (c))) 496 > (mpfr_exp_t) MPC_MAX_PREC (c) / 2) 497 return mpc_mul_naive (a, b, c, rnd); 498 else 499 return ((MPC_MAX_PREC(a) 500 <= (mpfr_prec_t) MUL_KARATSUBA_THRESHOLD * BITS_PER_MP_LIMB) 501 ? mpc_mul_naive : mpc_mul_karatsuba) (a, b, c, rnd); 502} 503