1/* mpc_sqrt -- Take the square root of a complex number.
2
3Copyright (C) 2002, 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#if MPFR_VERSION_MAJOR < 3
24#define mpfr_min_prec(x) \
25   ( ((prec + BITS_PER_MP_LIMB - 1) / BITS_PER_MP_LIMB) * BITS_PER_MP_LIMB \
26     - mpn_scan1 (x->_mpfr_d, 0))
27#endif
28
29
30int
31mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd)
32{
33  int ok_w, ok_t = 0;
34  mpfr_t    w, t;
35  mpfr_rnd_t  rnd_w, rnd_t;
36  mpfr_prec_t prec_w, prec_t;
37  /* the rounding mode and the precision required for w and t, which can */
38  /* be either the real or the imaginary part of a */
39  mpfr_prec_t prec;
40  int inex_w, inex_t = 1, inex_re, inex_im, loops = 0;
41  const int re_cmp = mpfr_cmp_ui (mpc_realref (b), 0),
42            im_cmp = mpfr_cmp_ui (mpc_imagref (b), 0);
43     /* comparison of the real/imaginary part of b with 0 */
44  int repr_w, repr_t = 0 /* to avoid gcc warning */ ;
45     /* flag indicating whether the computed value is already representable
46        at the target precision */
47  const int im_sgn = mpfr_signbit (mpc_imagref (b)) == 0 ? 0 : -1;
48     /* we need to know the sign of Im(b) when it is +/-0 */
49  const mpfr_rnd_t r = im_sgn ? GMP_RNDD : GMP_RNDU;
50     /* rounding mode used when computing t */
51
52  /* special values */
53  if (!mpc_fin_p (b)) {
54   /* sqrt(x +i*Inf) = +Inf +I*Inf, even if x = NaN */
55   /* sqrt(x -i*Inf) = +Inf -I*Inf, even if x = NaN */
56   if (mpfr_inf_p (mpc_imagref (b)))
57      {
58         mpfr_set_inf (mpc_realref (a), +1);
59         mpfr_set_inf (mpc_imagref (a), im_sgn);
60         return MPC_INEX (0, 0);
61      }
62
63   if (mpfr_inf_p (mpc_realref (b)))
64      {
65         if (mpfr_signbit (mpc_realref (b)))
66         {
67            if (mpfr_number_p (mpc_imagref (b)))
68               {
69               /* sqrt(-Inf +i*y) = +0 +i*Inf, when y positive */
70               /* sqrt(-Inf +i*y) = +0 -i*Inf, when y positive */
71               mpfr_set_ui (mpc_realref (a), 0, GMP_RNDN);
72               mpfr_set_inf (mpc_imagref (a), im_sgn);
73               return MPC_INEX (0, 0);
74               }
75            else
76               {
77               /* sqrt(-Inf +i*NaN) = NaN +/-i*Inf */
78               mpfr_set_nan (mpc_realref (a));
79               mpfr_set_inf (mpc_imagref (a), im_sgn);
80               return MPC_INEX (0, 0);
81               }
82         }
83         else
84         {
85            if (mpfr_number_p (mpc_imagref (b)))
86               {
87               /* sqrt(+Inf +i*y) = +Inf +i*0, when y positive */
88               /* sqrt(+Inf +i*y) = +Inf -i*0, when y positive */
89               mpfr_set_inf (mpc_realref (a), +1);
90               mpfr_set_ui (mpc_imagref (a), 0, GMP_RNDN);
91               if (im_sgn)
92                  mpc_conj (a, a, MPC_RNDNN);
93               return MPC_INEX (0, 0);
94               }
95            else
96               {
97               /* sqrt(+Inf -i*Inf) = +Inf -i*Inf */
98               /* sqrt(+Inf +i*Inf) = +Inf +i*Inf */
99               /* sqrt(+Inf +i*NaN) = +Inf +i*NaN */
100               return mpc_set (a, b, rnd);
101               }
102         }
103      }
104
105   /* sqrt(x +i*NaN) = NaN +i*NaN, if x is not infinite */
106   /* sqrt(NaN +i*y) = NaN +i*NaN, if y is not infinite */
107   if (mpfr_nan_p (mpc_realref (b)) || mpfr_nan_p (mpc_imagref (b)))
108      {
109         mpfr_set_nan (mpc_realref (a));
110         mpfr_set_nan (mpc_imagref (a));
111         return MPC_INEX (0, 0);
112      }
113  }
114
115  /* purely real */
116  if (im_cmp == 0)
117    {
118      if (re_cmp == 0)
119        {
120          mpc_set_ui_ui (a, 0, 0, MPC_RNDNN);
121          if (im_sgn)
122            mpc_conj (a, a, MPC_RNDNN);
123          return MPC_INEX (0, 0);
124        }
125      else if (re_cmp > 0)
126        {
127          inex_w = mpfr_sqrt (mpc_realref (a), mpc_realref (b), MPC_RND_RE (rnd));
128          mpfr_set_ui (mpc_imagref (a), 0, GMP_RNDN);
129          if (im_sgn)
130            mpc_conj (a, a, MPC_RNDNN);
131          return MPC_INEX (inex_w, 0);
132        }
133      else
134        {
135          mpfr_init2 (w, MPC_PREC_RE (b));
136          mpfr_neg (w, mpc_realref (b), GMP_RNDN);
137          if (im_sgn)
138            {
139              inex_w = -mpfr_sqrt (mpc_imagref (a), w, INV_RND (MPC_RND_IM (rnd)));
140              mpfr_neg (mpc_imagref (a), mpc_imagref (a), GMP_RNDN);
141            }
142          else
143            inex_w = mpfr_sqrt (mpc_imagref (a), w, MPC_RND_IM (rnd));
144
145          mpfr_set_ui (mpc_realref (a), 0, GMP_RNDN);
146          mpfr_clear (w);
147          return MPC_INEX (0, inex_w);
148        }
149    }
150
151  /* purely imaginary */
152  if (re_cmp == 0)
153    {
154      mpfr_t y;
155
156      y[0] = mpc_imagref (b)[0];
157      /* If y/2 underflows, so does sqrt(y/2) */
158      mpfr_div_2ui (y, y, 1, GMP_RNDN);
159      if (im_cmp > 0)
160        {
161          inex_w = mpfr_sqrt (mpc_realref (a), y, MPC_RND_RE (rnd));
162          inex_t = mpfr_sqrt (mpc_imagref (a), y, MPC_RND_IM (rnd));
163        }
164      else
165        {
166          mpfr_neg (y, y, GMP_RNDN);
167          inex_w = mpfr_sqrt (mpc_realref (a), y, MPC_RND_RE (rnd));
168          inex_t = -mpfr_sqrt (mpc_imagref (a), y, INV_RND (MPC_RND_IM (rnd)));
169          mpfr_neg (mpc_imagref (a), mpc_imagref (a), GMP_RNDN);
170        }
171      return MPC_INEX (inex_w, inex_t);
172    }
173
174  prec = MPC_MAX_PREC(a);
175
176  mpfr_init (w);
177  mpfr_init (t);
178
179   if (re_cmp > 0) {
180      rnd_w = MPC_RND_RE (rnd);
181      prec_w = MPC_PREC_RE (a);
182      rnd_t = MPC_RND_IM(rnd);
183      if (rnd_t == GMP_RNDZ)
184         /* force GMP_RNDD or GMP_RNDUP, using sign(t) = sign(y) */
185         rnd_t = (im_cmp > 0 ? GMP_RNDD : GMP_RNDU);
186      prec_t = MPC_PREC_IM (a);
187   }
188   else {
189      prec_w = MPC_PREC_IM (a);
190      prec_t = MPC_PREC_RE (a);
191      if (im_cmp > 0) {
192         rnd_w = MPC_RND_IM(rnd);
193         rnd_t = MPC_RND_RE(rnd);
194         if (rnd_t == GMP_RNDZ)
195            rnd_t = GMP_RNDD;
196      }
197      else {
198         rnd_w = INV_RND(MPC_RND_IM (rnd));
199         rnd_t = INV_RND(MPC_RND_RE (rnd));
200         if (rnd_t == GMP_RNDZ)
201            rnd_t = GMP_RNDU;
202      }
203   }
204
205  do
206    {
207      loops ++;
208      prec += (loops <= 2) ? mpc_ceil_log2 (prec) + 4 : prec / 2;
209      mpfr_set_prec (w, prec);
210      mpfr_set_prec (t, prec);
211      /* let b = x + iy */
212      /* w = sqrt ((|x| + sqrt (x^2 + y^2)) / 2), rounded down */
213      /* total error bounded by 3 ulps */
214      inex_w = mpc_abs (w, b, GMP_RNDD);
215      if (re_cmp < 0)
216        inex_w |= mpfr_sub (w, w, mpc_realref (b), GMP_RNDD);
217      else
218        inex_w |= mpfr_add (w, w, mpc_realref (b), GMP_RNDD);
219      inex_w |= mpfr_div_2ui (w, w, 1, GMP_RNDD);
220      inex_w |= mpfr_sqrt (w, w, GMP_RNDD);
221
222      repr_w = mpfr_min_prec (w) <= prec_w;
223      if (!repr_w)
224         /* use the usual trick for obtaining the ternary value */
225         ok_w = mpfr_can_round (w, prec - 2, GMP_RNDD, GMP_RNDU,
226                                prec_w + (rnd_w == GMP_RNDN));
227      else {
228            /* w is representable in the target precision and thus cannot be
229               rounded up */
230         if (rnd_w == GMP_RNDN)
231            /* If w can be rounded to nearest, then actually no rounding
232               occurs, and the ternary value is known from inex_w. */
233            ok_w = mpfr_can_round (w, prec - 2, GMP_RNDD, GMP_RNDN, prec_w);
234         else
235            /* If w can be rounded down, then any direct rounding and the
236               ternary flag can be determined from inex_w. */
237            ok_w = mpfr_can_round (w, prec - 2, GMP_RNDD, GMP_RNDD, prec_w);
238      }
239
240      if (!inex_w || ok_w) {
241         /* t = y / 2w, rounded away */
242         /* total error bounded by 7 ulps */
243         inex_t = mpfr_div (t, mpc_imagref (b), w, r);
244         if (!inex_t && inex_w)
245            /* The division was exact, but w was not. */
246            inex_t = im_sgn ? -1 : 1;
247         inex_t |= mpfr_div_2ui (t, t, 1, r);
248         repr_t = mpfr_min_prec (t) <= prec_t;
249         if (!repr_t)
250             /* As for w; since t was rounded away, we check whether rounding to 0
251                is possible. */
252            ok_t = mpfr_can_round (t, prec - 3, r, GMP_RNDZ,
253                                   prec_t + (rnd_t == GMP_RNDN));
254         else {
255            if (rnd_t == GMP_RNDN)
256               ok_t = mpfr_can_round (t, prec - 3, r, GMP_RNDN, prec_t);
257            else
258               ok_t = mpfr_can_round (t, prec - 3, r, r, prec_t);
259         }
260      }
261    }
262    while ((inex_w && !ok_w) || (inex_t && !ok_t));
263
264   if (re_cmp > 0) {
265         inex_re = mpfr_set (mpc_realref (a), w, MPC_RND_RE(rnd));
266         inex_im = mpfr_set (mpc_imagref (a), t, MPC_RND_IM(rnd));
267   }
268   else if (im_cmp > 0) {
269      inex_re = mpfr_set (mpc_realref(a), t, MPC_RND_RE(rnd));
270      inex_im = mpfr_set (mpc_imagref(a), w, MPC_RND_IM(rnd));
271   }
272   else {
273      inex_re = mpfr_neg (mpc_realref (a), t, MPC_RND_RE(rnd));
274      inex_im = mpfr_neg (mpc_imagref (a), w, MPC_RND_IM(rnd));
275   }
276
277   if (repr_w && inex_w) {
278      if (rnd_w == GMP_RNDN) {
279         /* w has not been rounded with mpfr_set/mpfr_neg, determine ternary
280            value from inex_w instead */
281         if (re_cmp > 0)
282            inex_re = inex_w;
283         else if (im_cmp > 0)
284            inex_im = inex_w;
285         else
286            inex_im = -inex_w;
287      }
288      else {
289         /* determine ternary value, but also potentially add 1 ulp; can only
290            be done now when we are in the target precision */
291         if (re_cmp > 0) {
292            if (rnd_w == GMP_RNDU) {
293               MPFR_ADD_ONE_ULP (mpc_realref (a));
294               inex_re = +1;
295            }
296            else
297               inex_re = -1;
298         }
299         else if (im_cmp > 0) {
300            if (rnd_w == GMP_RNDU) {
301               MPFR_ADD_ONE_ULP (mpc_imagref (a));
302               inex_im = +1;
303            }
304            else
305               inex_im = -1;
306         }
307         else {
308            if (rnd_w == GMP_RNDU) {
309               MPFR_ADD_ONE_ULP (mpc_imagref (a));
310               inex_im = -1;
311            }
312            else
313               inex_im = +1;
314         }
315      }
316   }
317   if (repr_t && inex_t) {
318      if (rnd_t == GMP_RNDN) {
319         if (re_cmp > 0)
320            inex_im = inex_t;
321         else if (im_cmp > 0)
322            inex_re = inex_t;
323         else
324            inex_re = -inex_t;
325      }
326      else {
327         if (re_cmp > 0) {
328            if (rnd_t == r)
329               inex_im = inex_t;
330            else {
331               inex_im = -inex_t;
332               /* im_cmp > 0 implies that Im(b) > 0, thus im_sgn = 0
333                  and r = GMP_RNDU.
334                  im_cmp < 0 implies that Im(b) < 0, thus im_sgn = -1
335                  and r = GMP_RNDD. */
336               MPFR_SUB_ONE_ULP (mpc_imagref (a));
337            }
338         }
339         else if (im_cmp > 0) {
340            if (rnd_t == r)
341               inex_re = inex_t;
342            else {
343               inex_re = -inex_t;
344               /* im_cmp > 0 implies r = GMP_RNDU (see above) */
345               MPFR_SUB_ONE_ULP (mpc_realref (a));
346            }
347         }
348         else { /* im_cmp < 0 */
349            if (rnd_t == r)
350               inex_re = -inex_t;
351            else {
352               inex_re = inex_t;
353               /* im_cmp < 0 implies r = GMP_RNDD (see above) */
354               MPFR_SUB_ONE_ULP (mpc_realref (a));
355            }
356         }
357      }
358   }
359
360  mpfr_clear (w);
361  mpfr_clear (t);
362
363  return MPC_INEX (inex_re, inex_im);
364}
365