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
);