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