1 /* mpc_sqrt -- Take the square root of a complex number.
3 Copyright (C) 2002, 2008, 2009, 2010, 2011, 2012 INRIA
5 This file is part of GNU MPC.
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.
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
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/ .
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))
31 mpc_sqrt (mpc_ptr a, mpc_srcptr b, mpc_rnd_t rnd)
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 */
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 */
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)))
58 mpfr_set_inf (mpc_realref (a), +1);
59 mpfr_set_inf (mpc_imagref (a), im_sgn);
60 return MPC_INEX (0, 0);
63 if (mpfr_inf_p (mpc_realref (b)))
65 if (mpfr_signbit (mpc_realref (b)))
67 if (mpfr_number_p (mpc_imagref (b)))
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);
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);
85 if (mpfr_number_p (mpc_imagref (b)))
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);
92 mpc_conj (a, a, MPC_RNDNN);
93 return MPC_INEX (0, 0);
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);
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)))
109 mpfr_set_nan (mpc_realref (a));
110 mpfr_set_nan (mpc_imagref (a));
111 return MPC_INEX (0, 0);
120 mpc_set_ui_ui (a, 0, 0, MPC_RNDNN);
122 mpc_conj (a, a, MPC_RNDNN);
123 return MPC_INEX (0, 0);
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);
130 mpc_conj (a, a, MPC_RNDNN);
131 return MPC_INEX (inex_w, 0);
135 mpfr_init2 (w, MPC_PREC_RE (b));
136 mpfr_neg (w, mpc_realref (b), GMP_RNDN);
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);
143 inex_w = mpfr_sqrt (mpc_imagref (a), w, MPC_RND_IM (rnd));
145 mpfr_set_ui (mpc_realref (a), 0, GMP_RNDN);
147 return MPC_INEX (0, inex_w);
151 /* purely imaginary */
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);
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));
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);
171 return MPC_INEX (inex_w, inex_t);
174 prec = MPC_MAX_PREC(a);
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);
189 prec_w = MPC_PREC_IM (a);
190 prec_t = MPC_PREC_RE (a);
192 rnd_w = MPC_RND_IM(rnd);
193 rnd_t = MPC_RND_RE(rnd);
194 if (rnd_t == GMP_RNDZ)
198 rnd_w = INV_RND(MPC_RND_IM (rnd));
199 rnd_t = INV_RND(MPC_RND_RE (rnd));
200 if (rnd_t == GMP_RNDZ)
208 prec += (loops <= 2) ? mpc_ceil_log2 (prec) + 4 : prec / 2;
209 mpfr_set_prec (w, prec);
210 mpfr_set_prec (t, prec);
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);
216 inex_w |= mpfr_sub (w, w, mpc_realref (b), GMP_RNDD);
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);
222 repr_w = mpfr_min_prec (w) <= prec_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));
228 /* w is representable in the target precision and thus cannot be
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);
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);
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;
250 /* As for w; since t was rounded away, we check whether rounding to 0
252 ok_t = mpfr_can_round (t, prec - 3, r, GMP_RNDZ,
253 prec_t + (rnd_t == GMP_RNDN));
255 if (rnd_t == GMP_RNDN)
256 ok_t = mpfr_can_round (t, prec - 3, r, GMP_RNDN, prec_t);
258 ok_t = mpfr_can_round (t, prec - 3, r, r, prec_t);
262 while ((inex_w && !ok_w) || (inex_t && !ok_t));
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));
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));
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));
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 */
289 /* determine ternary value, but also potentially add 1 ulp; can only
290 be done now when we are in the target precision */
292 if (rnd_w == GMP_RNDU) {
293 MPFR_ADD_ONE_ULP (mpc_realref (a));
299 else if (im_cmp > 0) {
300 if (rnd_w == GMP_RNDU) {
301 MPFR_ADD_ONE_ULP (mpc_imagref (a));
308 if (rnd_w == GMP_RNDU) {
309 MPFR_ADD_ONE_ULP (mpc_imagref (a));
317 if (repr_t && inex_t) {
318 if (rnd_t == GMP_RNDN) {
332 /* im_cmp > 0 implies that Im(b) > 0, thus im_sgn = 0
334 im_cmp < 0 implies that Im(b) < 0, thus im_sgn = -1
336 MPFR_SUB_ONE_ULP (mpc_imagref (a));
339 else if (im_cmp > 0) {
344 /* im_cmp > 0 implies r = GMP_RNDU (see above) */
345 MPFR_SUB_ONE_ULP (mpc_realref (a));
348 else { /* im_cmp < 0 */
353 /* im_cmp < 0 implies r = GMP_RNDD (see above) */
354 MPFR_SUB_ONE_ULP (mpc_realref (a));
363 return MPC_INEX (inex_re, inex_im);