1 /* mpc_atan -- arctangent of a complex number.
3 Copyright (C) 2009, 2010, 2011, 2012, 2013 INRIA
5 This file is part of GNU MPC.
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.
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
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/ .
27 rounded in the direction rnd
30 set_pi_over_2 (mpfr_ptr rop
, int s
, mpfr_rnd_t rnd
)
34 inex
= mpfr_const_pi (rop
, s
< 0 ? INV_RND (rnd
) : rnd
);
35 mpfr_div_2ui (rop
, rop
, 1, GMP_RNDN
);
39 mpfr_neg (rop
, rop
, GMP_RNDN
);
46 mpc_atan (mpc_ptr rop
, mpc_srcptr op
, mpc_rnd_t rnd
)
56 s_re
= mpfr_signbit (mpc_realref (op
));
57 s_im
= mpfr_signbit (mpc_imagref (op
));
60 if (mpfr_nan_p (mpc_realref (op
)) || mpfr_nan_p (mpc_imagref (op
)))
62 if (mpfr_nan_p (mpc_realref (op
)))
64 mpfr_set_nan (mpc_realref (rop
));
65 if (mpfr_zero_p (mpc_imagref (op
)) || mpfr_inf_p (mpc_imagref (op
)))
67 mpfr_set_ui (mpc_imagref (rop
), 0, GMP_RNDN
);
69 mpc_conj (rop
, rop
, MPC_RNDNN
);
72 mpfr_set_nan (mpc_imagref (rop
));
76 if (mpfr_inf_p (mpc_realref (op
)))
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
);
83 mpfr_set_nan (mpc_realref (rop
));
84 mpfr_set_nan (mpc_imagref (rop
));
87 return MPC_INEX (inex_re
, 0);
90 if (mpfr_inf_p (mpc_realref (op
)) || mpfr_inf_p (mpc_imagref (op
)))
92 inex_re
= set_pi_over_2 (mpc_realref (rop
), -s_re
, MPC_RND_RE (rnd
));
94 mpfr_set_ui (mpc_imagref (rop
), 0, GMP_RNDN
);
96 mpc_conj (rop
, rop
, GMP_RNDN
);
98 return MPC_INEX (inex_re
, 0);
101 /* pure real argument */
102 if (mpfr_zero_p (mpc_imagref (op
)))
104 inex_re
= mpfr_atan (mpc_realref (rop
), mpc_realref (op
), MPC_RND_RE (rnd
));
106 mpfr_set_ui (mpc_imagref (rop
), 0, GMP_RNDN
);
108 mpc_conj (rop
, rop
, GMP_RNDN
);
110 return MPC_INEX (inex_re
, 0);
113 /* pure imaginary argument */
114 if (mpfr_zero_p (mpc_realref (op
)))
119 cmp_1
= -mpfr_cmp_si (mpc_imagref (op
), -1);
121 cmp_1
= mpfr_cmp_ui (mpc_imagref (op
), +1);
125 /* atan(+0+iy) = +0 +i*atanh(y), if |y| < 1
126 atan(-0+iy) = -0 +i*atanh(y), if |y| < 1 */
128 mpfr_set_ui (mpc_realref (rop
), 0, GMP_RNDN
);
130 mpfr_neg (mpc_realref (rop
), mpc_realref (rop
), GMP_RNDN
);
132 inex_im
= mpfr_atanh (mpc_imagref (rop
), mpc_imagref (op
), MPC_RND_IM (rnd
));
136 /* atan(+/-0 +i) = +/-0 +i*inf
137 atan(+/-0 -i) = +/-0 -i*inf */
138 mpfr_set_zero (mpc_realref (rop
), s_re
? -1 : +1);
139 mpfr_set_inf (mpc_imagref (rop
), s_im
? -1 : +1);
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
;
150 rnd_im
= MPC_RND_IM (rnd
);
152 p_im
= mpfr_get_prec (mpc_imagref (rop
));
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)
158 As |atanh (1/y)| > |1/y| we have Exp(a)-Exp(b) <=0 so, at most,
159 2 bits of precision are lost.
161 We round atanh(1/y) away from 0.
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. */
175 /* atanh cannot underflow: |atanh(x)| > |x| for |x| < 1 */
176 inex_im
|= mpfr_atanh (y
, y
, rnd_away
);
179 || mpfr_can_round (y
, p
- 2, rnd_away
, GMP_RNDZ
,
180 p_im
+ (rnd_im
== GMP_RNDN
));
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
);
187 return MPC_INEX (inex_re
, inex_im
);
190 /* regular number argument */
194 mpfr_exp_t err
, expo
;
197 mpfr_exp_t op_re_exp
, op_im_exp
;
198 mpfr_rnd_t rnd1
, rnd2
;
200 mpfr_inits2 (MPFR_PREC_MIN
, a
, b
, x
, y
, (mpfr_ptr
) 0);
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
));
208 prec
= mpfr_get_prec (mpc_realref (rop
)); /* result precision */
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)
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)
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))
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
;
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
);
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 */
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 */
251 err
= mpfr_get_exp (a
); /* err = Exp(a) with the notations above */
252 mpfr_atan2 (x
, mpc_realref (op
), a
, GMP_RNDU
);
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)
261 MPC_ASSERT (mpfr_cmp_si (mpc_imagref(op
), -1) == 0);
262 /* check for intermediate underflow */
263 expo
= err
; /* will leave err unchanged below */
266 expo
= mpfr_get_exp (a
); /* expo = Exp(c) with the notations above */
267 mpfr_atan2 (b
, minus_op_re
, a
, GMP_RNDD
);
269 err
= err
< expo
? err
: expo
; /* err = min(Exp(a),Exp(c)) */
270 mpfr_sub (x
, x
, b
, GMP_RNDU
);
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;
276 mpfr_div_2ui (x
, x
, 1, GMP_RNDU
);
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
));
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 */
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)
306 p
= prec
; /* working precision */
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
);
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
);
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
329 mpfr_add (b
, b
, y
, GMP_RNDZ
);
330 mpfr_log (b
, b
, GMP_RNDZ
);
332 mpfr_sub (y
, a
, b
, GMP_RNDU
);
335 /* FIXME: happens when x and y have very different magnitudes;
336 could be handled more efficiently */
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;
347 err
= (expo
< 0) ? 1 : expo
+ 2;
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,
357 ok
= mpfr_can_round (y
, p
- err
, GMP_RNDU
, GMP_RNDD
,
358 prec
+ (MPC_RND_IM (rnd
) == GMP_RNDN
));
362 inex
= mpc_set_fr_fr (rop
, x
, y
, rnd
);
364 mpfr_clears (a
, b
, x
, y
, (mpfr_ptr
) 0);