1 /* mpc_atan -- arctangent of a complex number. 2 3 Copyright (C) 2009, 2010, 2011, 2012, 2013, 2017, 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 <stdio.h> 22 #include "mpc-impl.h" 23 24 /* set rop to 25 -pi/2 if s < 0 26 +pi/2 else 27 rounded in the direction rnd 28 */ 29 int 30 set_pi_over_2 (mpfr_ptr rop, int s, mpfr_rnd_t rnd) 31 { 32 int inex; 33 34 inex = mpfr_const_pi (rop, s < 0 ? INV_RND (rnd) : rnd); 35 mpfr_div_2ui (rop, rop, 1, MPFR_RNDN); 36 if (s < 0) 37 { 38 inex = -inex; 39 mpfr_neg (rop, rop, MPFR_RNDN); 40 } 41 42 return inex; 43 } 44 45 int 46 mpc_atan (mpc_ptr rop, mpc_srcptr op, mpc_rnd_t rnd) 47 { 48 int s_re, s_im; 49 int inex_re, inex_im, inex; 50 mpfr_exp_t saved_emin, saved_emax; 51 52 inex_re = 0; 53 inex_im = 0; 54 s_re = mpfr_signbit (mpc_realref (op)); 55 s_im = mpfr_signbit (mpc_imagref (op)); 56 57 /* special values */ 58 if (mpfr_nan_p (mpc_realref (op)) || mpfr_nan_p (mpc_imagref (op))) 59 { 60 if (mpfr_nan_p (mpc_realref (op))) 61 { 62 mpfr_set_nan (mpc_realref (rop)); 63 if (mpfr_zero_p (mpc_imagref (op)) || mpfr_inf_p (mpc_imagref (op))) 64 { 65 mpfr_set_ui (mpc_imagref (rop), 0, MPFR_RNDN); 66 if (s_im) 67 mpc_conj (rop, rop, MPC_RNDNN); 68 } 69 else 70 mpfr_set_nan (mpc_imagref (rop)); 71 } 72 else 73 { 74 if (mpfr_inf_p (mpc_realref (op))) 75 { 76 inex_re = set_pi_over_2 (mpc_realref (rop), -s_re, MPC_RND_RE (rnd)); 77 mpfr_set_ui (mpc_imagref (rop), 0, MPFR_RNDN); 78 } 79 else 80 { 81 mpfr_set_nan (mpc_realref (rop)); 82 mpfr_set_nan (mpc_imagref (rop)); 83 } 84 } 85 return MPC_INEX (inex_re, 0); 86 } 87 88 if (mpfr_inf_p (mpc_realref (op)) || mpfr_inf_p (mpc_imagref (op))) 89 { 90 inex_re = set_pi_over_2 (mpc_realref (rop), -s_re, MPC_RND_RE (rnd)); 91 92 mpfr_set_ui (mpc_imagref (rop), 0, MPFR_RNDN); 93 if (s_im) 94 mpc_conj (rop, rop, MPFR_RNDN); 95 96 return MPC_INEX (inex_re, 0); 97 } 98 99 /* pure real argument */ 100 if (mpfr_zero_p (mpc_imagref (op))) 101 { 102 inex_re = mpfr_atan (mpc_realref (rop), mpc_realref (op), MPC_RND_RE (rnd)); 103 104 mpfr_set_ui (mpc_imagref (rop), 0, MPFR_RNDN); 105 if (s_im) 106 mpc_conj (rop, rop, MPFR_RNDN); 107 108 return MPC_INEX (inex_re, 0); 109 } 110 111 /* pure imaginary argument */ 112 if (mpfr_zero_p (mpc_realref (op))) 113 { 114 int cmp_1; 115 116 if (s_im) 117 cmp_1 = -mpfr_cmp_si (mpc_imagref (op), -1); 118 else 119 cmp_1 = mpfr_cmp_ui (mpc_imagref (op), +1); 120 121 if (cmp_1 < 0) 122 { 123 /* atan(+0+iy) = +0 +i*atanh(y), if |y| < 1 124 atan(-0+iy) = -0 +i*atanh(y), if |y| < 1 */ 125 126 mpfr_set_ui (mpc_realref (rop), 0, MPFR_RNDN); 127 if (s_re) 128 mpfr_neg (mpc_realref (rop), mpc_realref (rop), MPFR_RNDN); 129 130 inex_im = mpfr_atanh (mpc_imagref (rop), mpc_imagref (op), MPC_RND_IM (rnd)); 131 } 132 else if (cmp_1 == 0) 133 { 134 /* atan(+/-0 +i) = +/-0 +i*inf 135 atan(+/-0 -i) = +/-0 -i*inf */ 136 mpfr_set_zero (mpc_realref (rop), s_re ? -1 : +1); 137 mpfr_set_inf (mpc_imagref (rop), s_im ? -1 : +1); 138 } 139 else 140 { 141 /* atan(+0+iy) = +pi/2 +i*atanh(1/y), if |y| > 1 142 atan(-0+iy) = -pi/2 +i*atanh(1/y), if |y| > 1 */ 143 mpfr_rnd_t rnd_im, rnd_away; 144 mpfr_t y, z; 145 mpfr_prec_t p, p_im; 146 int ok = 0; 147 148 rnd_im = MPC_RND_IM (rnd); 149 mpfr_init (y); 150 mpfr_init (z); 151 p_im = mpfr_get_prec (mpc_imagref (rop)); 152 p = p_im; 153 154 /* a = o(1/y) with error(a) < ulp(a), rounded away 155 b = o(atanh(a)) with error(b) < ulp(b) + 1/|a^2-1|*ulp(a), 156 since if a = 1/y + eps, then atanh(a) = atanh(1/y) + eps * atanh'(t) 157 with t in (1/y, a). Since a is rounded away, we have 1/y <= a <= 1 158 if y > 1, and -1 <= a <= 1/y if y < -1, thus |atanh'(t)| = 159 1/|t^2-1| <= 1/|a^2-1|. 160 161 We round atanh(1/y) away from 0. 162 */ 163 do 164 { 165 mpfr_exp_t err, exp_a; 166 167 p += mpc_ceil_log2 (p) + 2; 168 mpfr_set_prec (y, p); 169 mpfr_set_prec (z, p); 170 rnd_away = s_im == 0 ? MPFR_RNDU : MPFR_RNDD; 171 inex_im = mpfr_ui_div (y, 1, mpc_imagref (op), rnd_away); 172 exp_a = mpfr_get_exp (y); 173 /* FIXME: should we consider the case with unreasonably huge 174 precision prec(y)>3*exp_min, where atanh(1/Im(op)) could be 175 representable while 1/Im(op) underflows ? 176 This corresponds to |y| = 0.5*2^emin, in which case the 177 result may be wrong. */ 178 179 /* We would like to compute a rounded-up error bound 1/|a^2-1|, 180 so we need to round down |a^2-1|, which means rounding up 181 a^2 since |a|<1. */ 182 mpfr_sqr (z, y, MPFR_RNDU); 183 /* since |y| > 1, we should have |a| <= 1, thus a^2 <= 1 */ 184 MPC_ASSERT(mpfr_cmp_ui (z, 1) <= 0); 185 /* in case z=1, we should try again with more precision */ 186 if (mpfr_cmp_ui (z, 1) == 0) 187 continue; 188 /* now z < 1 */ 189 mpfr_ui_sub (z, 1, z, MPFR_RNDZ); 190 191 /* atanh cannot underflow: |atanh(x)| > |x| for |x| < 1 */ 192 inex_im |= mpfr_atanh (y, y, rnd_away); 193 194 /* the error is now bounded by ulp(b) + 1/z*ulp(a), thus 195 ulp(b) + 2^(exp(a) - exp(b) + 1 - exp(z)) * ulp(b) */ 196 err = exp_a - mpfr_get_exp (y) + 1 - mpfr_get_exp (z); 197 if (err >= 0) /* 1 + 2^err <= 2^(err+1) */ 198 err = err + 1; 199 else 200 err = 1; /* 1 + 2^err <= 2^1 */ 201 202 /* the error is bounded by 2^err ulps */ 203 204 ok = inex_im == 0 205 || mpfr_can_round (y, p - err, rnd_away, MPFR_RNDZ, 206 p_im + (rnd_im == MPFR_RNDN)); 207 } while (ok == 0); 208 209 inex_re = set_pi_over_2 (mpc_realref (rop), -s_re, MPC_RND_RE (rnd)); 210 inex_im = mpfr_set (mpc_imagref (rop), y, rnd_im); 211 mpfr_clear (y); 212 mpfr_clear (z); 213 } 214 return MPC_INEX (inex_re, inex_im); 215 } 216 217 saved_emin = mpfr_get_emin (); 218 saved_emax = mpfr_get_emax (); 219 mpfr_set_emin (mpfr_get_emin_min ()); 220 mpfr_set_emax (mpfr_get_emax_max ()); 221 222 /* regular number argument */ 223 { 224 mpfr_t a, b, x, y; 225 mpfr_prec_t prec, p; 226 mpfr_exp_t err, expo; 227 int ok = 0; 228 mpfr_t minus_op_re; 229 mpfr_exp_t op_re_exp, op_im_exp; 230 mpfr_rnd_t rnd1, rnd2; 231 232 mpfr_inits2 (MPFR_PREC_MIN, a, b, x, y, (mpfr_ptr) 0); 233 234 /* real part: Re(arctan(x+i*y)) = [arctan2(x,1-y) - arctan2(-x,1+y)]/2 */ 235 minus_op_re[0] = mpc_realref (op)[0]; 236 MPFR_CHANGE_SIGN (minus_op_re); 237 op_re_exp = mpfr_get_exp (mpc_realref (op)); 238 op_im_exp = mpfr_get_exp (mpc_imagref (op)); 239 240 prec = mpfr_get_prec (mpc_realref (rop)); /* result precision */ 241 242 /* a = o(1-y) error(a) < 1 ulp(a) 243 b = o(atan2(x,a)) error(b) < [1+2^{3+Exp(x)-Exp(a)-Exp(b)}] ulp(b) 244 = kb ulp(b) 245 c = o(1+y) error(c) < 1 ulp(c) 246 d = o(atan2(-x,c)) error(d) < [1+2^{3+Exp(x)-Exp(c)-Exp(d)}] ulp(d) 247 = kd ulp(d) 248 e = o(b - d) error(e) < [1 + kb*2^{Exp(b}-Exp(e)} 249 + kd*2^{Exp(d)-Exp(e)}] ulp(e) 250 error(e) < [1 + 2^{4+Exp(x)-Exp(a)-Exp(e)} 251 + 2^{4+Exp(x)-Exp(c)-Exp(e)}] ulp(e) 252 because |atan(u)| < |u| 253 < [1 + 2^{5+Exp(x)-min(Exp(a),Exp(c)) 254 -Exp(e)}] ulp(e) 255 f = e/2 exact 256 */ 257 258 /* p: working precision */ 259 p = (op_im_exp > 0 || prec > SAFE_ABS (mpfr_prec_t, op_im_exp)) ? prec 260 : (prec - op_im_exp); 261 rnd1 = mpfr_sgn (mpc_realref (op)) > 0 ? MPFR_RNDD : MPFR_RNDU; 262 rnd2 = mpfr_sgn (mpc_realref (op)) < 0 ? MPFR_RNDU : MPFR_RNDD; 263 264 do 265 { 266 p += mpc_ceil_log2 (p) + 2; 267 mpfr_set_prec (a, p); 268 mpfr_set_prec (b, p); 269 mpfr_set_prec (x, p); 270 271 /* x = upper bound for atan (x/(1-y)). Since atan is increasing, we 272 need an upper bound on x/(1-y), i.e., a lower bound on 1-y for 273 x positive, and an upper bound on 1-y for x negative */ 274 mpfr_ui_sub (a, 1, mpc_imagref (op), rnd1); 275 if (mpfr_sgn (a) == 0) /* y is near 1, thus 1+y is near 2, and 276 expo will be 1 or 2 below */ 277 { 278 MPC_ASSERT (mpfr_cmp_ui (mpc_imagref(op), 1) == 0); 279 /* check for intermediate underflow */ 280 err = 2; /* ensures err will be expo below */ 281 } 282 else 283 err = mpfr_get_exp (a); /* err = Exp(a) with the notations above */ 284 mpfr_atan2 (x, mpc_realref (op), a, MPFR_RNDU); 285 286 /* b = lower bound for atan (-x/(1+y)): for x negative, we need a 287 lower bound on -x/(1+y), i.e., an upper bound on 1+y */ 288 mpfr_add_ui (a, mpc_imagref(op), 1, rnd2); 289 /* if a is exactly zero, i.e., Im(op) = -1, then the error on a is 0, 290 and we can simply ignore the terms involving Exp(a) in the error */ 291 if (mpfr_sgn (a) == 0) 292 { 293 MPC_ASSERT (mpfr_cmp_si (mpc_imagref(op), -1) == 0); 294 /* check for intermediate underflow */ 295 expo = err; /* will leave err unchanged below */ 296 } 297 else 298 expo = mpfr_get_exp (a); /* expo = Exp(c) with the notations above */ 299 mpfr_atan2 (b, minus_op_re, a, MPFR_RNDD); 300 301 err = err < expo ? err : expo; /* err = min(Exp(a),Exp(c)) */ 302 mpfr_sub (x, x, b, MPFR_RNDU); 303 304 err = 5 + op_re_exp - err - mpfr_get_exp (x); 305 /* error is bounded by [1 + 2^err] ulp(e) */ 306 err = err < 0 ? 1 : err + 1; 307 308 mpfr_div_2ui (x, x, 1, MPFR_RNDU); 309 310 /* Note: using RND2=RNDD guarantees that if x is exactly representable 311 on prec + ... bits, mpfr_can_round will return 0 */ 312 ok = mpfr_can_round (x, p - err, MPFR_RNDU, MPFR_RNDD, 313 prec + (MPC_RND_RE (rnd) == MPFR_RNDN)); 314 } while (ok == 0); 315 316 /* Imaginary part 317 Im(atan(x+I*y)) = 1/4 * [log(x^2+(1+y)^2) - log (x^2 +(1-y)^2)] */ 318 prec = mpfr_get_prec (mpc_imagref (rop)); /* result precision */ 319 320 /* a = o(1+y) error(a) < 1 ulp(a) 321 b = o(a^2) error(b) < 5 ulp(b) 322 c = o(x^2) error(c) < 1 ulp(c) 323 d = o(b+c) error(d) < 7 ulp(d) 324 e = o(log(d)) error(e) < [1 + 7*2^{2-Exp(e)}] ulp(e) = ke ulp(e) 325 f = o(1-y) error(f) < 1 ulp(f) 326 g = o(f^2) error(g) < 5 ulp(g) 327 h = o(c+f) error(h) < 7 ulp(h) 328 i = o(log(h)) error(i) < [1 + 7*2^{2-Exp(i)}] ulp(i) = ki ulp(i) 329 j = o(e-i) error(j) < [1 + ke*2^{Exp(e)-Exp(j)} 330 + ki*2^{Exp(i)-Exp(j)}] ulp(j) 331 error(j) < [1 + 2^{Exp(e)-Exp(j)} + 2^{Exp(i)-Exp(j)} 332 + 7*2^{3-Exp(j)}] ulp(j) 333 < [1 + 2^{max(Exp(e),Exp(i))-Exp(j)+1} 334 + 7*2^{3-Exp(j)}] ulp(j) 335 k = j/4 exact 336 */ 337 err = 2; 338 p = prec; /* working precision */ 339 340 do 341 { 342 p += mpc_ceil_log2 (p) + err; 343 mpfr_set_prec (a, p); 344 mpfr_set_prec (b, p); 345 mpfr_set_prec (y, p); 346 347 /* a = upper bound for log(x^2 + (1+y)^2) */ 348 mpfr_add_ui (a, mpc_imagref (op), 1, MPFR_RNDA); 349 mpfr_sqr (a, a, MPFR_RNDU); 350 mpfr_sqr (y, mpc_realref (op), MPFR_RNDU); 351 mpfr_add (a, a, y, MPFR_RNDU); 352 mpfr_log (a, a, MPFR_RNDU); 353 354 /* b = lower bound for log(x^2 + (1-y)^2) */ 355 mpfr_ui_sub (b, 1, mpc_imagref (op), MPFR_RNDZ); /* round to zero */ 356 mpfr_sqr (b, b, MPFR_RNDZ); 357 /* we could write mpfr_sqr (y, mpc_realref (op), MPFR_RNDZ) but it is 358 more efficient to reuse the value of y (x^2) above and subtract 359 one ulp */ 360 mpfr_nextbelow (y); 361 mpfr_add (b, b, y, MPFR_RNDZ); 362 mpfr_log (b, b, MPFR_RNDZ); 363 364 mpfr_sub (y, a, b, MPFR_RNDU); 365 366 if (mpfr_zero_p (y)) 367 /* FIXME: happens when x and y have very different magnitudes; 368 could be handled more efficiently */ 369 ok = 0; 370 else 371 { 372 expo = MPC_MAX (mpfr_get_exp (a), mpfr_get_exp (b)); 373 expo = expo - mpfr_get_exp (y) + 1; 374 err = 3 - mpfr_get_exp (y); 375 /* error(j) <= [1 + 2^expo + 7*2^err] ulp(j) */ 376 if (expo <= err) /* error(j) <= [1 + 2^{err+1}] ulp(j) */ 377 err = (err < 0) ? 1 : err + 2; 378 else 379 err = (expo < 0) ? 1 : expo + 2; 380 381 mpfr_div_2ui (y, y, 2, MPFR_RNDN); 382 MPC_ASSERT (!mpfr_zero_p (y)); 383 /* FIXME: underflow. Since the main term of the Taylor series 384 in y=0 is 1/(x^2+1) * y, this means that y is very small 385 and/or x very large; but then the mpfr_zero_p (y) above 386 should be true. This needs a proof, or better yet, 387 special code. */ 388 389 ok = mpfr_can_round (y, p - err, MPFR_RNDU, MPFR_RNDD, 390 prec + (MPC_RND_IM (rnd) == MPFR_RNDN)); 391 } 392 } while (ok == 0); 393 394 inex = mpc_set_fr_fr (rop, x, y, rnd); 395 396 mpfr_clears (a, b, x, y, (mpfr_ptr) 0); 397 398 /* restore the exponent range, and check the range of results */ 399 mpfr_set_emin (saved_emin); 400 mpfr_set_emax (saved_emax); 401 inex_re = mpfr_check_range (mpc_realref (rop), MPC_INEX_RE (inex), 402 MPC_RND_RE (rnd)); 403 inex_im = mpfr_check_range (mpc_imagref (rop), MPC_INEX_IM (inex), 404 MPC_RND_IM (rnd)); 405 406 return MPC_INEX (inex_re, inex_im); 407 } 408 } 409