1 /* mpc_div -- Divide two complex numbers.
2
3 Copyright (C) 2002, 2003, 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 "mpc-impl.h"
22
23 /* this routine deals with the case where w is zero */
24 static int
mpc_div_zero(mpc_ptr a,mpc_srcptr z,mpc_srcptr w,mpc_rnd_t rnd)25 mpc_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 */
40 static int
mpc_div_inf_fin(mpc_ptr rop,mpc_srcptr z,mpc_srcptr w)41 mpc_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 */
112 static int
mpc_div_fin_inf(mpc_ptr rop,mpc_srcptr z,mpc_srcptr w)113 mpc_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
156 static int
mpc_div_real(mpc_ptr rop,mpc_srcptr z,mpc_srcptr w,mpc_rnd_t rnd)157 mpc_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
185 static int
mpc_div_imag(mpc_ptr rop,mpc_srcptr z,mpc_srcptr w,mpc_rnd_t rnd)186 mpc_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
233 int
mpc_div(mpc_ptr a,mpc_srcptr b,mpc_srcptr c,mpc_rnd_t rnd)234 mpc_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