1 /* mpc_sqrt -- Take the square root of a complex number. 2 3 Copyright (C) 2002, 2008, 2009, 2010, 2011, 2012 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 #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 30 int 31 mpc_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