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 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 prec_t = MPC_PREC_IM (a); 181 } 182 else { 183 prec_w = MPC_PREC_IM (a); 184 prec_t = MPC_PREC_RE (a); 185 if (im_cmp > 0) { 186 rnd_w = MPC_RND_IM(rnd); 187 rnd_t = MPC_RND_RE(rnd); 188 if (rnd_t == MPFR_RNDZ) 189 rnd_t = MPFR_RNDD; 190 } 191 else { 192 rnd_w = INV_RND(MPC_RND_IM (rnd)); 193 rnd_t = INV_RND(MPC_RND_RE (rnd)); 194 if (rnd_t == MPFR_RNDZ) 195 rnd_t = MPFR_RNDU; 196 } 197 } 198 199 saved_emin = mpfr_get_emin (); 200 saved_emax = mpfr_get_emax (); 201 mpfr_set_emin (mpfr_get_emin_min ()); 202 mpfr_set_emax (mpfr_get_emax_max ()); 203 204 do 205 { 206 loops ++; 207 prec += (loops <= 2) ? mpc_ceil_log2 (prec) + 4 : prec / 2; 208 mpfr_set_prec (w, prec); 209 mpfr_set_prec (t, prec); 210 /* let b = x + iy */ 211 /* w = sqrt ((|x| + sqrt (x^2 + y^2)) / 2), rounded down */ 212 /* final error on w bounded by 10 ulps, see algorithms.tex */ 213 inex_w = mpfr_sqr (w, mpc_realref (b), MPFR_RNDD); 214 inex_w |= mpfr_sqr (t, mpc_imagref (b), MPFR_RNDD); 215 inex_w |= mpfr_add (w, w, t, MPFR_RNDD); 216 inex_w |= mpfr_sqrt (w, w, MPFR_RNDD); 217 if (re_cmp < 0) 218 inex_w |= mpfr_sub (w, w, mpc_realref (b), MPFR_RNDD); 219 else 220 inex_w |= mpfr_add (w, w, mpc_realref (b), MPFR_RNDD); 221 inex_w |= mpfr_div_2ui (w, w, 1, MPFR_RNDD); 222 inex_w |= mpfr_sqrt (w, w, MPFR_RNDD); 223 224 repr_w = mpfr_min_prec (w) <= prec_w; 225 if (!repr_w) 226 /* use the usual trick for obtaining the ternary value */ 227 ok_w = mpfr_can_round (w, prec - 4, MPFR_RNDD, MPFR_RNDU, 228 prec_w + (rnd_w == MPFR_RNDN)); 229 else { 230 /* w is representable in the target precision and thus cannot be 231 rounded up */ 232 if (rnd_w == MPFR_RNDN) 233 /* If w can be rounded to nearest, then actually no rounding 234 occurs, and the ternary value is known from inex_w. */ 235 ok_w = mpfr_can_round (w, prec - 4, MPFR_RNDD, MPFR_RNDN, prec_w); 236 else 237 /* If w can be rounded down, then any direct rounding and the 238 ternary flag can be determined from inex_w. */ 239 ok_w = mpfr_can_round (w, prec - 4, MPFR_RNDD, MPFR_RNDD, prec_w); 240 } 241 242 if (!inex_w || ok_w) { 243 /* t = y / 2w, rounded away */ 244 /* total error bounded by 16 ulps, see algorithms.tex */ 245 inex_t = mpfr_div (t, mpc_imagref (b), w, r); 246 if (!inex_t && inex_w) 247 /* The division was exact, but w was not. */ 248 inex_t = im_sgn ? -1 : 1; 249 inex_t |= mpfr_div_2ui (t, t, 1, r); 250 repr_t = mpfr_min_prec (t) <= prec_t; 251 if (!repr_t) 252 /* As for w; since t was rounded away, we check whether rounding to 0 253 is possible. */ 254 ok_t = mpfr_can_round (t, prec - 4, r, MPFR_RNDZ, 255 prec_t + (rnd_t == MPFR_RNDN)); 256 else { 257 if (rnd_t == MPFR_RNDN) 258 ok_t = mpfr_can_round (t, prec - 4, r, MPFR_RNDN, prec_t); 259 else 260 ok_t = mpfr_can_round (t, prec - 4, r, r, prec_t); 261 } 262 } 263 } 264 while ((inex_w && !ok_w) || (inex_t && !ok_t)); 265 266 if (re_cmp > 0) { 267 inex_re = mpfr_set (mpc_realref (a), w, MPC_RND_RE(rnd)); 268 inex_im = mpfr_set (mpc_imagref (a), t, MPC_RND_IM(rnd)); 269 } 270 else if (im_cmp > 0) { 271 inex_re = mpfr_set (mpc_realref(a), t, MPC_RND_RE(rnd)); 272 inex_im = mpfr_set (mpc_imagref(a), w, MPC_RND_IM(rnd)); 273 } 274 else { 275 inex_re = mpfr_neg (mpc_realref (a), t, MPC_RND_RE(rnd)); 276 inex_im = mpfr_neg (mpc_imagref (a), w, MPC_RND_IM(rnd)); 277 } 278 279 if (repr_w && inex_w) { 280 if (rnd_w == MPFR_RNDN) { 281 /* w has not been rounded with mpfr_set/mpfr_neg, determine ternary 282 value from inex_w instead */ 283 if (re_cmp > 0) 284 inex_re = inex_w; 285 else if (im_cmp > 0) 286 inex_im = inex_w; 287 else 288 inex_im = -inex_w; 289 } 290 else { 291 /* determine ternary value, but also potentially add 1 ulp; can only 292 be done now when we are in the target precision */ 293 if (re_cmp > 0) { 294 if (rnd_w == MPFR_RNDU) { 295 MPFR_ADD_ONE_ULP (mpc_realref (a)); 296 inex_re = +1; 297 } 298 else 299 inex_re = -1; 300 } 301 else if (im_cmp > 0) { 302 if (rnd_w == MPFR_RNDU) { 303 MPFR_ADD_ONE_ULP (mpc_imagref (a)); 304 inex_im = +1; 305 } 306 else 307 inex_im = -1; 308 } 309 else { 310 if (rnd_w == MPFR_RNDU) { 311 MPFR_ADD_ONE_ULP (mpc_imagref (a)); 312 inex_im = -1; 313 } 314 else 315 inex_im = +1; 316 } 317 } 318 } 319 if (repr_t && inex_t) { 320 if (rnd_t == MPFR_RNDN) { 321 if (re_cmp > 0) 322 inex_im = inex_t; 323 else if (im_cmp > 0) 324 inex_re = inex_t; 325 else 326 inex_re = -inex_t; 327 } 328 else { 329 if (re_cmp > 0) { 330 if (rnd_t == r) 331 inex_im = inex_t; 332 else { 333 inex_im = -inex_t; 334 /* im_cmp > 0 implies that Im(b) > 0, thus im_sgn = 0 335 and r = MPFR_RNDU. 336 im_cmp < 0 implies that Im(b) < 0, thus im_sgn = -1 337 and r = MPFR_RNDD. */ 338 MPFR_SUB_ONE_ULP (mpc_imagref (a)); 339 } 340 } 341 else if (im_cmp > 0) { 342 if (rnd_t == r) 343 inex_re = inex_t; 344 else { 345 inex_re = -inex_t; 346 /* im_cmp > 0 implies r = MPFR_RNDU (see above) */ 347 MPFR_SUB_ONE_ULP (mpc_realref (a)); 348 } 349 } 350 else { /* im_cmp < 0 */ 351 if (rnd_t == r) 352 inex_re = -inex_t; 353 else { 354 inex_re = inex_t; 355 /* im_cmp < 0 implies r = MPFR_RNDD (see above) */ 356 MPFR_SUB_ONE_ULP (mpc_realref (a)); 357 } 358 } 359 } 360 } 361 362 mpfr_clear (w); 363 mpfr_clear (t); 364 365 /* restore the exponent range, and check the range of results */ 366 mpfr_set_emin (saved_emin); 367 mpfr_set_emax (saved_emax); 368 inex_re = mpfr_check_range (mpc_realref (a), inex_re, MPC_RND_RE (rnd)); 369 inex_im = mpfr_check_range (mpc_imagref (a), inex_im, MPC_RND_IM (rnd)); 370 371 return MPC_INEX (inex_re, inex_im); 372 } 373