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