diff options
Diffstat (limited to 'mpc/src/mul.c')
| -rw-r--r-- | mpc/src/mul.c | 639 |
1 files changed, 639 insertions, 0 deletions
diff --git a/mpc/src/mul.c b/mpc/src/mul.c new file mode 100644 index 0000000000..2be9b8d646 --- /dev/null +++ b/mpc/src/mul.c @@ -0,0 +1,639 @@ +/* mpc_mul -- Multiply two complex numbers + +Copyright (C) 2002, 2004, 2005, 2008, 2009, 2010, 2011, 2012 INRIA + +This file is part of GNU MPC. + +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 +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. + +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/ . +*/ + +#include <stdio.h> /* for MPC_ASSERT */ +#include "mpc-impl.h" + +#define mpz_add_si(z,x,y) do { \ + if (y >= 0) \ + mpz_add_ui (z, x, (long int) y); \ + else \ + mpz_sub_ui (z, x, (long int) (-y)); \ + } while (0); + +/* compute z=x*y when x has an infinite part */ +static int +mul_infinite (mpc_ptr z, mpc_srcptr x, mpc_srcptr y) +{ + /* Let x=xr+i*xi and y=yr+i*yi; extract the signs of the operands */ + int xrs = mpfr_signbit (mpc_realref (x)) ? -1 : 1; + int xis = mpfr_signbit (mpc_imagref (x)) ? -1 : 1; + int yrs = mpfr_signbit (mpc_realref (y)) ? -1 : 1; + int yis = mpfr_signbit (mpc_imagref (y)) ? -1 : 1; + + int u, v; + + /* compute the sign of + u = xrs * yrs * xr * yr - xis * yis * xi * yi + v = xrs * yis * xr * yi + xis * yrs * xi * yr + +1 if positive, -1 if negatiye, 0 if NaN */ + if ( mpfr_nan_p (mpc_realref (x)) || mpfr_nan_p (mpc_imagref (x)) + || mpfr_nan_p (mpc_realref (y)) || mpfr_nan_p (mpc_imagref (y))) { + u = 0; + v = 0; + } + else if (mpfr_inf_p (mpc_realref (x))) { + /* x = (+/-inf) xr + i*xi */ + u = ( mpfr_zero_p (mpc_realref (y)) + || (mpfr_inf_p (mpc_imagref (x)) && mpfr_zero_p (mpc_imagref (y))) + || (mpfr_zero_p (mpc_imagref (x)) && mpfr_inf_p (mpc_imagref (y))) + || ( (mpfr_inf_p (mpc_imagref (x)) || mpfr_inf_p (mpc_imagref (y))) + && xrs*yrs == xis*yis) + ? 0 : xrs * yrs); + v = ( mpfr_zero_p (mpc_imagref (y)) + || (mpfr_inf_p (mpc_imagref (x)) && mpfr_zero_p (mpc_realref (y))) + || (mpfr_zero_p (mpc_imagref (x)) && mpfr_inf_p (mpc_realref (y))) + || ( (mpfr_inf_p (mpc_imagref (x)) || mpfr_inf_p (mpc_imagref (x))) + && xrs*yis != xis*yrs) + ? 0 : xrs * yis); + } + else { + /* x = xr + i*(+/-inf) with |xr| != inf */ + u = ( mpfr_zero_p (mpc_imagref (y)) + || (mpfr_zero_p (mpc_realref (x)) && mpfr_inf_p (mpc_realref (y))) + || (mpfr_inf_p (mpc_realref (y)) && xrs*yrs == xis*yis) + ? 0 : -xis * yis); + v = ( mpfr_zero_p (mpc_realref (y)) + || (mpfr_zero_p (mpc_realref (x)) && mpfr_inf_p (mpc_imagref (y))) + || (mpfr_inf_p (mpc_imagref (y)) && xrs*yis != xis*yrs) + ? 0 : xis * yrs); + } + + if (u == 0 && v == 0) { + /* Naive result is NaN+i*NaN. Obtain an infinity using the algorithm + given in Annex G.5.1 of the ISO C99 standard */ + int xr = (mpfr_zero_p (mpc_realref (x)) || mpfr_nan_p (mpc_realref (x)) ? 0 + : (mpfr_inf_p (mpc_realref (x)) ? 1 : 0)); + int xi = (mpfr_zero_p (mpc_imagref (x)) || mpfr_nan_p (mpc_imagref (x)) ? 0 + : (mpfr_inf_p (mpc_imagref (x)) ? 1 : 0)); + int yr = (mpfr_zero_p (mpc_realref (y)) || mpfr_nan_p (mpc_realref (y)) ? 0 : 1); + int yi = (mpfr_zero_p (mpc_imagref (y)) || mpfr_nan_p (mpc_imagref (y)) ? 0 : 1); + if (mpc_inf_p (y)) { + yr = mpfr_inf_p (mpc_realref (y)) ? 1 : 0; + yi = mpfr_inf_p (mpc_imagref (y)) ? 1 : 0; + } + + u = xrs * xr * yrs * yr - xis * xi * yis * yi; + v = xrs * xr * yis * yi + xis * xi * yrs * yr; + } + + if (u == 0) + mpfr_set_nan (mpc_realref (z)); + else + mpfr_set_inf (mpc_realref (z), u); + + if (v == 0) + mpfr_set_nan (mpc_imagref (z)); + else + mpfr_set_inf (mpc_imagref (z), v); + + return MPC_INEX (0, 0); /* exact */ +} + + +/* compute z = x*y for Im(y) == 0 */ +static int +mul_real (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) +{ + int xrs, xis, yrs, yis; + int inex; + + /* save signs of operands */ + xrs = MPFR_SIGNBIT (mpc_realref (x)); + xis = MPFR_SIGNBIT (mpc_imagref (x)); + yrs = MPFR_SIGNBIT (mpc_realref (y)); + yis = MPFR_SIGNBIT (mpc_imagref (y)); + + inex = mpc_mul_fr (z, x, mpc_realref (y), rnd); + /* Signs of zeroes may be wrong. Their correction does not change the + inexact flag. */ + if (mpfr_zero_p (mpc_realref (z))) + mpfr_setsign (mpc_realref (z), mpc_realref (z), MPC_RND_RE(rnd) == GMP_RNDD + || (xrs != yrs && xis == yis), GMP_RNDN); + if (mpfr_zero_p (mpc_imagref (z))) + mpfr_setsign (mpc_imagref (z), mpc_imagref (z), MPC_RND_IM (rnd) == GMP_RNDD + || (xrs != yis && xis != yrs), GMP_RNDN); + + return inex; +} + + +/* compute z = x*y for Re(y) == 0, and Im(x) != 0 and Im(y) != 0 */ +static int +mul_imag (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) +{ + int sign; + int inex_re, inex_im; + int overlap = z == x || z == y; + mpc_t rop; + + if (overlap) + mpc_init3 (rop, MPC_PREC_RE (z), MPC_PREC_IM (z)); + else + rop [0] = z[0]; + + sign = (MPFR_SIGNBIT (mpc_realref (y)) != MPFR_SIGNBIT (mpc_imagref (x))) + && (MPFR_SIGNBIT (mpc_imagref (y)) != MPFR_SIGNBIT (mpc_realref (x))); + + inex_re = -mpfr_mul (mpc_realref (rop), mpc_imagref (x), mpc_imagref (y), + INV_RND (MPC_RND_RE (rnd))); + mpfr_neg (mpc_realref (rop), mpc_realref (rop), GMP_RNDN); /* exact */ + inex_im = mpfr_mul (mpc_imagref (rop), mpc_realref (x), mpc_imagref (y), + MPC_RND_IM (rnd)); + mpc_set (z, rop, MPC_RNDNN); + + /* Sign of zeroes may be wrong (note that Re(z) cannot be zero) */ + if (mpfr_zero_p (mpc_imagref (z))) + mpfr_setsign (mpc_imagref (z), mpc_imagref (z), MPC_RND_IM (rnd) == GMP_RNDD + || sign, GMP_RNDN); + + if (overlap) + mpc_clear (rop); + + return MPC_INEX (inex_re, inex_im); +} + + +static int +mpfr_fmma (mpfr_ptr z, mpfr_srcptr a, mpfr_srcptr b, mpfr_srcptr c, + mpfr_srcptr d, int sign, mpfr_rnd_t rnd) +{ + /* Computes z = ab+cd if sign >= 0, or z = ab-cd if sign < 0. + Assumes that a, b, c, d are finite and non-zero; so any multiplication + of two of them yielding an infinity is an overflow, and a + multiplication yielding 0 is an underflow. + Assumes further that z is distinct from a, b, c, d. */ + + int inex; + mpfr_t u, v; + + /* u=a*b, v=sign*c*d exactly */ + mpfr_init2 (u, mpfr_get_prec (a) + mpfr_get_prec (b)); + mpfr_init2 (v, mpfr_get_prec (c) + mpfr_get_prec (d)); + mpfr_mul (u, a, b, GMP_RNDN); + mpfr_mul (v, c, d, GMP_RNDN); + if (sign < 0) + mpfr_neg (v, v, GMP_RNDN); + + /* tentatively compute z as u+v; here we need z to be distinct + from a, b, c, d to not lose the latter */ + inex = mpfr_add (z, u, v, rnd); + + if (mpfr_inf_p (z)) { + /* replace by "correctly rounded overflow" */ + mpfr_set_si (z, (mpfr_signbit (z) ? -1 : 1), GMP_RNDN); + inex = mpfr_mul_2ui (z, z, mpfr_get_emax (), rnd); + } + else if (mpfr_zero_p (u) && !mpfr_zero_p (v)) { + /* exactly u underflowed, determine inexact flag */ + inex = (mpfr_signbit (u) ? 1 : -1); + } + else if (mpfr_zero_p (v) && !mpfr_zero_p (u)) { + /* exactly v underflowed, determine inexact flag */ + inex = (mpfr_signbit (v) ? 1 : -1); + } + else if (mpfr_nan_p (z) || (mpfr_zero_p (u) && mpfr_zero_p (v))) { + /* In the first case, u and v are infinities with opposite signs. + In the second case, u and v are zeroes; their sum may be 0 or the + least representable number, with a sign to be determined. + Redo the computations with mpz_t exponents */ + mpfr_exp_t ea, eb, ec, ed; + mpz_t eu, ev; + /* cheat to work around the const qualifiers */ + + /* Normalise the input by shifting and keep track of the shifts in + the exponents of u and v */ + ea = mpfr_get_exp (a); + eb = mpfr_get_exp (b); + ec = mpfr_get_exp (c); + ed = mpfr_get_exp (d); + + mpfr_set_exp ((mpfr_ptr) a, (mpfr_prec_t) 0); + mpfr_set_exp ((mpfr_ptr) b, (mpfr_prec_t) 0); + mpfr_set_exp ((mpfr_ptr) c, (mpfr_prec_t) 0); + mpfr_set_exp ((mpfr_ptr) d, (mpfr_prec_t) 0); + + mpz_init (eu); + mpz_init (ev); + mpz_set_si (eu, (long int) ea); + mpz_add_si (eu, eu, (long int) eb); + mpz_set_si (ev, (long int) ec); + mpz_add_si (ev, ev, (long int) ed); + + /* recompute u and v and move exponents to eu and ev */ + mpfr_mul (u, a, b, GMP_RNDN); + /* exponent of u is non-positive */ + mpz_sub_ui (eu, eu, (unsigned long int) (-mpfr_get_exp (u))); + mpfr_set_exp (u, (mpfr_prec_t) 0); + mpfr_mul (v, c, d, GMP_RNDN); + if (sign < 0) + mpfr_neg (v, v, GMP_RNDN); + mpz_sub_ui (ev, ev, (unsigned long int) (-mpfr_get_exp (v))); + mpfr_set_exp (v, (mpfr_prec_t) 0); + + if (mpfr_nan_p (z)) { + mpfr_exp_t emax = mpfr_get_emax (); + int overflow; + /* We have a = ma * 2^ea with 1/2 <= |ma| < 1 and ea <= emax, and + analogously for b. So eu <= 2*emax, and eu > emax since we have + an overflow. The same holds for ev. Shift u and v by as much as + possible so that one of them has exponent emax and the + remaining exponents in eu and ev are the same. Then carry out + the addition. Shifting u and v prevents an underflow. */ + if (mpz_cmp (eu, ev) >= 0) { + mpfr_set_exp (u, emax); + mpz_sub_ui (eu, eu, (long int) emax); + mpz_sub (ev, ev, eu); + mpfr_set_exp (v, (mpfr_exp_t) mpz_get_ui (ev)); + /* remaining common exponent is now in eu */ + } + else { + mpfr_set_exp (v, emax); + mpz_sub_ui (ev, ev, (long int) emax); + mpz_sub (eu, eu, ev); + mpfr_set_exp (u, (mpfr_exp_t) mpz_get_ui (eu)); + mpz_set (eu, ev); + /* remaining common exponent is now also in eu */ + } + inex = mpfr_add (z, u, v, rnd); + /* Result is finite since u and v have different signs. */ + overflow = mpfr_mul_2ui (z, z, mpz_get_ui (eu), rnd); + if (overflow) + inex = overflow; + } + else { + int underflow; + /* Addition of two zeroes with same sign. We have a = ma * 2^ea + with 1/2 <= |ma| < 1 and ea >= emin and similarly for b. + So 2*emin < 2*emin+1 <= eu < emin < 0, and analogously for v. */ + mpfr_exp_t emin = mpfr_get_emin (); + if (mpz_cmp (eu, ev) <= 0) { + mpfr_set_exp (u, emin); + mpz_add_ui (eu, eu, (unsigned long int) (-emin)); + mpz_sub (ev, ev, eu); + mpfr_set_exp (v, (mpfr_exp_t) mpz_get_si (ev)); + } + else { + mpfr_set_exp (v, emin); + mpz_add_ui (ev, ev, (unsigned long int) (-emin)); + mpz_sub (eu, eu, ev); + mpfr_set_exp (u, (mpfr_exp_t) mpz_get_si (eu)); + mpz_set (eu, ev); + } + inex = mpfr_add (z, u, v, rnd); + mpz_neg (eu, eu); + underflow = mpfr_div_2ui (z, z, mpz_get_ui (eu), rnd); + if (underflow) + inex = underflow; + } + + mpz_clear (eu); + mpz_clear (ev); + + mpfr_set_exp ((mpfr_ptr) a, ea); + mpfr_set_exp ((mpfr_ptr) b, eb); + mpfr_set_exp ((mpfr_ptr) c, ec); + mpfr_set_exp ((mpfr_ptr) d, ed); + /* works also when some of a, b, c, d are not all distinct */ + } + + mpfr_clear (u); + mpfr_clear (v); + + return inex; +} + + +int +mpc_mul_naive (mpc_ptr z, mpc_srcptr x, mpc_srcptr y, mpc_rnd_t rnd) +{ + /* computes z=x*y by the schoolbook method, where x and y are assumed + to be finite and without zero parts */ + int overlap, inex; + mpc_t rop; + + MPC_ASSERT ( mpfr_regular_p (mpc_realref (x)) && mpfr_regular_p (mpc_imagref (x)) + && mpfr_regular_p (mpc_realref (y)) && mpfr_regular_p (mpc_imagref (y))); + overlap = (z == x) || (z == y); + if (overlap) + mpc_init3 (rop, MPC_PREC_RE (z), MPC_PREC_IM (z)); + else + rop [0] = z [0]; + + inex = MPC_INEX (mpfr_fmma (mpc_realref (rop), mpc_realref (x), mpc_realref (y), mpc_imagref (x), + mpc_imagref (y), -1, MPC_RND_RE (rnd)), + mpfr_fmma (mpc_imagref (rop), mpc_realref (x), mpc_imagref (y), mpc_imagref (x), + mpc_realref (y), +1, MPC_RND_IM (rnd))); + + mpc_set (z, rop, MPC_RNDNN); + if (overlap) + mpc_clear (rop); + + return inex; +} + + +int +mpc_mul_karatsuba (mpc_ptr rop, mpc_srcptr op1, mpc_srcptr op2, mpc_rnd_t rnd) +{ + /* computes rop=op1*op2 by a Karatsuba algorithm, where op1 and op2 + are assumed to be finite and without zero parts */ + mpfr_srcptr a, b, c, d; + int mul_i, ok, inexact, mul_a, mul_c, inex_re = 0, inex_im = 0, sign_x, sign_u; + mpfr_t u, v, w, x; + mpfr_prec_t prec, prec_re, prec_u, prec_v, prec_w; + mpfr_rnd_t rnd_re, rnd_u; + int overlap; + /* true if rop == op1 or rop == op2 */ + mpc_t result; + /* overlap is quite difficult to handle, because we have to tentatively + round the variable u in the end to either the real or the imaginary + part of rop (it is not possible to tell now whether the real or + imaginary part is used). If this fails, we have to start again and + need the correct values of op1 and op2. + So we just create a new variable for the result in this case. */ + int loop; + const int MAX_MUL_LOOP = 1; + + overlap = (rop == op1) || (rop == op2); + if (overlap) + mpc_init3 (result, MPC_PREC_RE (rop), MPC_PREC_IM (rop)); + else + result [0] = rop [0]; + + a = mpc_realref(op1); + b = mpc_imagref(op1); + c = mpc_realref(op2); + d = mpc_imagref(op2); + + /* (a + i*b) * (c + i*d) = [ac - bd] + i*[ad + bc] */ + + mul_i = 0; /* number of multiplications by i */ + mul_a = 1; /* implicit factor for a */ + mul_c = 1; /* implicit factor for c */ + + if (mpfr_cmp_abs (a, b) < 0) + { + MPFR_SWAP (a, b); + mul_i ++; + mul_a = -1; /* consider i * (a+i*b) = -b + i*a */ + } + + if (mpfr_cmp_abs (c, d) < 0) + { + MPFR_SWAP (c, d); + mul_i ++; + mul_c = -1; /* consider -d + i*c instead of c + i*d */ + } + + /* find the precision and rounding mode for the new real part */ + if (mul_i % 2) + { + prec_re = MPC_PREC_IM(rop); + rnd_re = MPC_RND_IM(rnd); + } + else /* mul_i = 0 or 2 */ + { + prec_re = MPC_PREC_RE(rop); + rnd_re = MPC_RND_RE(rnd); + } + + if (mul_i) + rnd_re = INV_RND(rnd_re); + + /* now |a| >= |b| and |c| >= |d| */ + prec = MPC_MAX_PREC(rop); + + mpfr_init2 (v, prec_v = mpfr_get_prec (a) + mpfr_get_prec (d)); + mpfr_init2 (w, prec_w = mpfr_get_prec (b) + mpfr_get_prec (c)); + mpfr_init2 (u, 2); + mpfr_init2 (x, 2); + + inexact = mpfr_mul (v, a, d, GMP_RNDN); + if (inexact) { + /* over- or underflow */ + ok = 0; + goto clear; + } + if (mul_a == -1) + mpfr_neg (v, v, GMP_RNDN); + + inexact = mpfr_mul (w, b, c, GMP_RNDN); + if (inexact) { + /* over- or underflow */ + ok = 0; + goto clear; + } + if (mul_c == -1) + mpfr_neg (w, w, GMP_RNDN); + + /* compute sign(v-w) */ + sign_x = mpfr_cmp_abs (v, w); + if (sign_x > 0) + sign_x = 2 * mpfr_sgn (v) - mpfr_sgn (w); + else if (sign_x == 0) + sign_x = mpfr_sgn (v) - mpfr_sgn (w); + else + sign_x = mpfr_sgn (v) - 2 * mpfr_sgn (w); + + sign_u = mul_a * mpfr_sgn (a) * mul_c * mpfr_sgn (c); + + if (sign_x * sign_u < 0) + { /* swap inputs */ + MPFR_SWAP (a, c); + MPFR_SWAP (b, d); + mpfr_swap (v, w); + { int tmp; tmp = mul_a; mul_a = mul_c; mul_c = tmp; } + sign_x = - sign_x; + } + + /* now sign_x * sign_u >= 0 */ + loop = 0; + do + { + loop++; + /* the following should give failures with prob. <= 1/prec */ + prec += mpc_ceil_log2 (prec) + 3; + + mpfr_set_prec (u, prec_u = prec); + mpfr_set_prec (x, prec); + + /* first compute away(b +/- a) and store it in u */ + inexact = (mul_a == -1 ? + ROUND_AWAY (mpfr_sub (u, b, a, MPFR_RNDA), u) : + ROUND_AWAY (mpfr_add (u, b, a, MPFR_RNDA), u)); + + /* then compute away(+/-c - d) and store it in x */ + inexact |= (mul_c == -1 ? + ROUND_AWAY (mpfr_add (x, c, d, MPFR_RNDA), x) : + ROUND_AWAY (mpfr_sub (x, c, d, MPFR_RNDA), x)); + if (mul_c == -1) + mpfr_neg (x, x, GMP_RNDN); + + if (inexact == 0) + mpfr_prec_round (u, prec_u = 2 * prec, GMP_RNDN); + + /* compute away(u*x) and store it in u */ + inexact |= ROUND_AWAY (mpfr_mul (u, u, x, MPFR_RNDA), u); + /* (a+b)*(c-d) */ + + /* if all computations are exact up to here, it may be that + the real part is exact, thus we need if possible to + compute v - w exactly */ + if (inexact == 0) + { + mpfr_prec_t prec_x; + /* v and w are different from 0, so mpfr_get_exp is safe to use */ + prec_x = SAFE_ABS (mpfr_exp_t, mpfr_get_exp (v) - mpfr_get_exp (w)) + + MPC_MAX (prec_v, prec_w) + 1; + /* +1 is necessary for a potential carry */ + /* ensure we do not use a too large precision */ + if (prec_x > prec_u) + prec_x = prec_u; + if (prec_x > prec) + mpfr_prec_round (x, prec_x, GMP_RNDN); + } + + rnd_u = (sign_u > 0) ? GMP_RNDU : GMP_RNDD; + inexact |= mpfr_sub (x, v, w, rnd_u); /* ad - bc */ + + /* in case u=0, ensure that rnd_u rounds x away from zero */ + if (mpfr_sgn (u) == 0) + rnd_u = (mpfr_sgn (x) > 0) ? GMP_RNDU : GMP_RNDD; + inexact |= mpfr_add (u, u, x, rnd_u); /* ac - bd */ + + ok = inexact == 0 || + mpfr_can_round (u, prec_u - 3, rnd_u, GMP_RNDZ, + prec_re + (rnd_re == GMP_RNDN)); + /* this ensures both we can round correctly and determine the correct + inexact flag (for rounding to nearest) */ + } + while (!ok && loop <= MAX_MUL_LOOP); + /* after MAX_MUL_LOOP rounds, use mpc_naive instead */ + + if (ok) { + /* if inexact is zero, then u is exactly ac-bd, otherwise fix the sign + of the inexact flag for u, which was rounded away from ac-bd */ + if (inexact != 0) + inexact = mpfr_sgn (u); + + if (mul_i == 0) + { + inex_re = mpfr_set (mpc_realref(result), u, MPC_RND_RE(rnd)); + if (inex_re == 0) + { + inex_re = inexact; /* u is rounded away from 0 */ + inex_im = mpfr_add (mpc_imagref(result), v, w, MPC_RND_IM(rnd)); + } + else + inex_im = mpfr_add (mpc_imagref(result), v, w, MPC_RND_IM(rnd)); + } + else if (mul_i == 1) /* (x+i*y)/i = y - i*x */ + { + inex_im = mpfr_neg (mpc_imagref(result), u, MPC_RND_IM(rnd)); + if (inex_im == 0) + { + inex_im = -inexact; /* u is rounded away from 0 */ + inex_re = mpfr_add (mpc_realref(result), v, w, MPC_RND_RE(rnd)); + } + else + inex_re = mpfr_add (mpc_realref(result), v, w, MPC_RND_RE(rnd)); + } + else /* mul_i = 2, z/i^2 = -z */ + { + inex_re = mpfr_neg (mpc_realref(result), u, MPC_RND_RE(rnd)); + if (inex_re == 0) + { + inex_re = -inexact; /* u is rounded away from 0 */ + inex_im = -mpfr_add (mpc_imagref(result), v, w, + INV_RND(MPC_RND_IM(rnd))); + mpfr_neg (mpc_imagref(result), mpc_imagref(result), MPC_RND_IM(rnd)); + } + else + { + inex_im = -mpfr_add (mpc_imagref(result), v, w, + INV_RND(MPC_RND_IM(rnd))); + mpfr_neg (mpc_imagref(result), mpc_imagref(result), MPC_RND_IM(rnd)); + } + } + + mpc_set (rop, result, MPC_RNDNN); + } + +clear: + mpfr_clear (u); + mpfr_clear (v); + mpfr_clear (w); + mpfr_clear (x); + if (overlap) + mpc_clear (result); + + if (ok) + return MPC_INEX(inex_re, inex_im); + else + return mpc_mul_naive (rop, op1, op2, rnd); +} + + +int +mpc_mul (mpc_ptr a, mpc_srcptr b, mpc_srcptr c, mpc_rnd_t rnd) +{ + /* Conforming to ISO C99 standard (G.5.1 multiplicative operators), + infinities are treated specially if both parts are NaN when computed + naively. */ + if (mpc_inf_p (b)) + return mul_infinite (a, b, c); + if (mpc_inf_p (c)) + return mul_infinite (a, c, b); + + /* NaN contamination of both parts in result */ + if (mpfr_nan_p (mpc_realref (b)) || mpfr_nan_p (mpc_imagref (b)) + || mpfr_nan_p (mpc_realref (c)) || mpfr_nan_p (mpc_imagref (c))) { + mpfr_set_nan (mpc_realref (a)); + mpfr_set_nan (mpc_imagref (a)); + return MPC_INEX (0, 0); + } + + /* check for real multiplication */ + if (mpfr_zero_p (mpc_imagref (b))) + return mul_real (a, c, b, rnd); + if (mpfr_zero_p (mpc_imagref (c))) + return mul_real (a, b, c, rnd); + + /* check for purely imaginary multiplication */ + if (mpfr_zero_p (mpc_realref (b))) + return mul_imag (a, c, b, rnd); + if (mpfr_zero_p (mpc_realref (c))) + return mul_imag (a, b, c, rnd); + + /* If the real and imaginary part of one argument have a very different */ + /* exponent, it is not reasonable to use Karatsuba multiplication. */ + if ( SAFE_ABS (mpfr_exp_t, + mpfr_get_exp (mpc_realref (b)) - mpfr_get_exp (mpc_imagref (b))) + > (mpfr_exp_t) MPC_MAX_PREC (b) / 2 + || SAFE_ABS (mpfr_exp_t, + mpfr_get_exp (mpc_realref (c)) - mpfr_get_exp (mpc_imagref (c))) + > (mpfr_exp_t) MPC_MAX_PREC (c) / 2) + return mpc_mul_naive (a, b, c, rnd); + else + return ((MPC_MAX_PREC(a) + <= (mpfr_prec_t) MUL_KARATSUBA_THRESHOLD * BITS_PER_MP_LIMB) + ? mpc_mul_naive : mpc_mul_karatsuba) (a, b, c, rnd); +} |