1 /* mpfr_atan2 -- arc-tan 2 of a floating-point number
3 Copyright 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramel projects, INRIA.
6 This file is part of the GNU MPFR Library.
8 The GNU MPFR Library is free software; you can redistribute it and/or modify
9 it under the terms of the GNU Lesser General Public License as published by
10 the Free Software Foundation; either version 3 of the License, or (at your
11 option) any later version.
13 The GNU MPFR Library is distributed in the hope that it will be useful, but
14 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
15 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
16 License for more details.
18 You should have received a copy of the GNU Lesser General Public License
19 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
20 http://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
21 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
23 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
27 pi_div_2ui (mpfr_ptr dest
, int i
, int neg
, mpfr_rnd_t rnd_mode
)
30 MPFR_SAVE_EXPO_DECL (expo
);
32 MPFR_SAVE_EXPO_MARK (expo
);
33 if (neg
) /* -PI/2^i */
35 inexact
= - mpfr_const_pi (dest
, MPFR_INVERT_RND (rnd_mode
));
36 MPFR_CHANGE_SIGN (dest
);
40 inexact
= mpfr_const_pi (dest
, rnd_mode
);
42 mpfr_div_2ui (dest
, dest
, i
, rnd_mode
); /* exact */
43 MPFR_SAVE_EXPO_FREE (expo
);
44 return mpfr_check_range (dest
, inexact
, rnd_mode
);
48 mpfr_atan2 (mpfr_ptr dest
, mpfr_srcptr y
, mpfr_srcptr x
, mpfr_rnd_t rnd_mode
)
54 MPFR_SAVE_EXPO_DECL (expo
);
58 (("y[%Pu]=%.*Rg x[%Pu]=%.*Rg rnd=%d",
59 mpfr_get_prec (y
), mpfr_log_prec
, y
,
60 mpfr_get_prec (x
), mpfr_log_prec
, x
, rnd_mode
),
61 ("atan[%Pu]=%.*Rg inexact=%d",
62 mpfr_get_prec (dest
), mpfr_log_prec
, dest
, inexact
));
65 if (MPFR_ARE_SINGULAR (x
, y
))
67 /* atan2(0, 0) does not raise the "invalid" floating-point
68 exception, nor does atan2(y, 0) raise the "divide-by-zero"
69 floating-point exception.
70 -- atan2(±0, -0) returns ±pi.313)
71 -- atan2(±0, +0) returns ±0.
72 -- atan2(±0, x) returns ±pi, for x < 0.
73 -- atan2(±0, x) returns ±0, for x > 0.
74 -- atan2(y, ±0) returns -pi/2 for y < 0.
75 -- atan2(y, ±0) returns pi/2 for y > 0.
76 -- atan2(±oo, -oo) returns ±3pi/4.
77 -- atan2(±oo, +oo) returns ±pi/4.
78 -- atan2(±oo, x) returns ±pi/2, for finite x.
79 -- atan2(±y, -oo) returns ±pi, for finite y > 0.
80 -- atan2(±y, +oo) returns ±0, for finite y > 0.
82 if (MPFR_IS_NAN (x
) || MPFR_IS_NAN (y
))
89 if (MPFR_IS_NEG (x
)) /* +/- PI */
94 inexact
= mpfr_const_pi (dest
, MPFR_INVERT_RND (rnd_mode
));
95 MPFR_CHANGE_SIGN (dest
);
99 return mpfr_const_pi (dest
, rnd_mode
);
104 MPFR_SET_ZERO (dest
);
105 MPFR_SET_SAME_SIGN (dest
, y
);
109 if (MPFR_IS_ZERO (x
))
111 return pi_div_2ui (dest
, 1, MPFR_IS_NEG (y
), rnd_mode
);
115 if (!MPFR_IS_INF (x
)) /* +/- PI/2 */
116 return pi_div_2ui (dest
, 1, MPFR_IS_NEG (y
), rnd_mode
);
117 else if (MPFR_IS_POS (x
)) /* +/- PI/4 */
118 return pi_div_2ui (dest
, 2, MPFR_IS_NEG (y
), rnd_mode
);
119 else /* +/- 3*PI/4: Ugly since we have to round properly */
122 MPFR_ZIV_DECL (loop2
);
123 mpfr_prec_t prec2
= MPFR_PREC (dest
) + 10;
125 MPFR_SAVE_EXPO_MARK (expo
);
126 mpfr_init2 (tmp2
, prec2
);
127 MPFR_ZIV_INIT (loop2
, prec2
);
130 mpfr_const_pi (tmp2
, MPFR_RNDN
);
131 mpfr_mul_ui (tmp2
, tmp2
, 3, MPFR_RNDN
); /* Error <= 2 */
132 mpfr_div_2ui (tmp2
, tmp2
, 2, MPFR_RNDN
);
133 if (mpfr_round_p (MPFR_MANT (tmp2
), MPFR_LIMB_SIZE (tmp2
),
134 MPFR_PREC (tmp2
) - 2,
135 MPFR_PREC (dest
) + (rnd_mode
== MPFR_RNDN
)))
137 MPFR_ZIV_NEXT (loop2
, prec2
);
138 mpfr_set_prec (tmp2
, prec2
);
140 MPFR_ZIV_FREE (loop2
);
142 MPFR_CHANGE_SIGN (tmp2
);
143 inexact
= mpfr_set (dest
, tmp2
, rnd_mode
);
145 MPFR_SAVE_EXPO_FREE (expo
);
146 return mpfr_check_range (dest
, inexact
, rnd_mode
);
149 MPFR_ASSERTD (MPFR_IS_INF (x
));
156 /* When x is a power of two, we call directly atan(y/x) since y/x is
158 if (MPFR_UNLIKELY (MPFR_IS_POWER_OF_2 (x
)))
162 unsigned int saved_flags
= __gmpfr_flags
;
164 mpfr_init2 (yoverx
, MPFR_PREC (y
));
165 if (MPFR_LIKELY (mpfr_div_2si (yoverx
, y
, MPFR_GET_EXP (x
) - 1,
168 /* Here the flags have not changed due to mpfr_div_2si. */
169 r
= mpfr_atan (dest
, yoverx
, rnd_mode
);
175 /* Division is inexact because of a small exponent range */
177 __gmpfr_flags
= saved_flags
;
181 MPFR_SAVE_EXPO_MARK (expo
);
183 /* Set up initial prec */
184 prec
= MPFR_PREC (dest
) + 3 + MPFR_INT_CEIL_LOG2 (MPFR_PREC (dest
));
185 mpfr_init2 (tmp
, prec
);
187 MPFR_ZIV_INIT (loop
, prec
);
189 /* use atan2(y,x) = atan(y/x) */
193 MPFR_BLOCK_DECL (flags
);
195 MPFR_BLOCK (flags
, div_inex
= mpfr_div (tmp
, y
, x
, MPFR_RNDN
));
198 /* Result is exact. */
199 inexact
= mpfr_atan (dest
, tmp
, rnd_mode
);
203 /* Error <= ulp (tmp) except in case of underflow or overflow. */
205 /* If the division underflowed, since |atan(z)/z| < 1, we have
207 if (MPFR_UNDERFLOW (flags
))
211 /* In the case MPFR_RNDN with 2^(emin-2) < |y/x| < 2^(emin-1):
212 The smallest significand value S > 1 of |y/x| is:
213 * 1 / (1 - 2^(-px)) if py <= px,
214 * (1 - 2^(-px) + 2^(-py)) / (1 - 2^(-px)) if py >= px.
215 Therefore S - 1 > 2^(-pz), where pz = max(px,py). We have:
216 atan(|y/x|) > atan(z), where z = 2^(emin-2) * (1 + 2^(-pz)).
218 > 2^(emin-2) * (1 + 2^(-pz) - 2^(2 emin - 5))
219 Assuming pz <= -2 emin + 5, we can round away from zero
220 (this is what mpfr_underflow always does on MPFR_RNDN).
221 In the case MPFR_RNDN with |y/x| <= 2^(emin-2), we round
222 toward zero, as |atan(z)/z| < 1. */
223 MPFR_ASSERTN (MPFR_PREC_MAX
<=
224 2 * (mpfr_uexp_t
) - MPFR_EMIN_MIN
+ 5);
225 if (rnd_mode
== MPFR_RNDN
&& MPFR_IS_ZERO (tmp
))
226 rnd_mode
= MPFR_RNDZ
;
227 sign
= MPFR_SIGN (tmp
);
229 MPFR_SAVE_EXPO_FREE (expo
);
230 return mpfr_underflow (dest
, rnd_mode
, sign
);
233 mpfr_atan (tmp
, tmp
, MPFR_RNDN
); /* Error <= 2*ulp (tmp) since
234 abs(D(arctan)) <= 1 */
235 /* TODO: check that the error bound is correct in case of overflow. */
236 /* FIXME: Error <= ulp(tmp) ? */
237 if (MPFR_LIKELY (MPFR_CAN_ROUND (tmp
, prec
- 2, MPFR_PREC (dest
),
240 MPFR_ZIV_NEXT (loop
, prec
);
241 mpfr_set_prec (tmp
, prec
);
244 /* Use sign(y)*(PI - atan (|y/x|)) */
246 mpfr_init2 (pi
, prec
);
249 mpfr_div (tmp
, y
, x
, MPFR_RNDN
); /* Error <= ulp (tmp) */
250 /* If tmp is 0, we have |y/x| <= 2^(-emin-2), thus
251 atan|y/x| < 2^(-emin-2). */
252 MPFR_SET_POS (tmp
); /* no error */
253 mpfr_atan (tmp
, tmp
, MPFR_RNDN
); /* Error <= 2*ulp (tmp) since
254 abs(D(arctan)) <= 1 */
255 mpfr_const_pi (pi
, MPFR_RNDN
); /* Error <= ulp(pi) /2 */
256 e
= MPFR_NOTZERO(tmp
) ? MPFR_GET_EXP (tmp
) : __gmpfr_emin
- 1;
257 mpfr_sub (tmp
, pi
, tmp
, MPFR_RNDN
); /* see above */
259 MPFR_CHANGE_SIGN (tmp
);
260 /* Error(tmp) <= (1/2+2^(EXP(pi)-EXP(tmp)-1)+2^(e-EXP(tmp)+1))*ulp
261 <= 2^(MAX (MAX (EXP(PI)-EXP(tmp)-1, e-EXP(tmp)+1),
263 e
= MAX (MAX (MPFR_GET_EXP (pi
)-MPFR_GET_EXP (tmp
) - 1,
264 e
- MPFR_GET_EXP (tmp
) + 1), -1) + 2;
265 if (MPFR_LIKELY (MPFR_CAN_ROUND (tmp
, prec
- e
, MPFR_PREC (dest
),
268 MPFR_ZIV_NEXT (loop
, prec
);
269 mpfr_set_prec (tmp
, prec
);
270 mpfr_set_prec (pi
, prec
);
274 inexact
= mpfr_set (dest
, tmp
, rnd_mode
);
277 MPFR_ZIV_FREE (loop
);
279 MPFR_SAVE_EXPO_FREE (expo
);
280 return mpfr_check_range (dest
, inexact
, rnd_mode
);