1 /* mpfr_atan2 -- arc-tan 2 of a floating-point number
3 Copyright 2005, 2006, 2007, 2008, 2009 Free Software Foundation, Inc.
4 Contributed by the Arenaire and Cacao projects, INRIA.
6 This file is part of the GNU MPFR Library, and was contributed by Mathieu Dutour.
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 2.1 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.LIB. If not, write to
20 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston,
21 MA 02110-1301, USA. */
23 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
27 mpfr_atan2 (mpfr_ptr dest
, mpfr_srcptr y
, mpfr_srcptr x
, mp_rnd_t rnd_mode
)
33 MPFR_SAVE_EXPO_DECL (expo
);
36 MPFR_LOG_FUNC (("y[%#R]=%R x[%#R]=%R rnd=%d", y
, y
, x
, x
, rnd_mode
),
37 ("atan[%#R]=%R inexact=%d", dest
, dest
, inexact
));
40 if (MPFR_ARE_SINGULAR (x
, y
))
42 /* atan2(0, 0) does not raise the "invalid" floating-point
43 exception, nor does atan2(y, 0) raise the "divide-by-zero"
44 floating-point exception.
45 -- atan2(±0, -0) returns ±pi.313)
46 -- atan2(±0, +0) returns ±0.
47 -- atan2(±0, x) returns ±pi, for x < 0.
48 -- atan2(±0, x) returns ±0, for x > 0.
49 -- atan2(y, ±0) returns -pi/2 for y < 0.
50 -- atan2(y, ±0) returns pi/2 for y > 0.
51 -- atan2(±oo, -oo) returns ±3pi/4.
52 -- atan2(±oo, +oo) returns ±pi/4.
53 -- atan2(±oo, x) returns ±pi/2, for finite x.
54 -- atan2(±y, -oo) returns ±pi, for finite y > 0.
55 -- atan2(±y, +oo) returns ±0, for finite y > 0.
57 if (MPFR_IS_NAN (x
) || MPFR_IS_NAN (y
))
64 if (MPFR_IS_NEG (x
)) /* +/- PI */
69 inexact
= mpfr_const_pi (dest
, MPFR_INVERT_RND (rnd_mode
));
70 MPFR_CHANGE_SIGN (dest
);
74 return mpfr_const_pi (dest
, rnd_mode
);
80 MPFR_SET_SAME_SIGN (dest
, y
);
87 if (MPFR_IS_NEG (y
)) /* -PI/2 */
89 inexact
= mpfr_const_pi (dest
, MPFR_INVERT_RND(rnd_mode
));
90 MPFR_CHANGE_SIGN (dest
);
91 mpfr_div_2ui (dest
, dest
, 1, rnd_mode
);
96 inexact
= mpfr_const_pi (dest
, rnd_mode
);
97 mpfr_div_2ui (dest
, dest
, 1, rnd_mode
);
103 if (!MPFR_IS_INF (x
)) /* +/- PI/2 */
105 else if (MPFR_IS_POS (x
)) /* +/- PI/4 */
109 rnd_mode
= MPFR_INVERT_RND (rnd_mode
);
110 inexact
= mpfr_const_pi (dest
, rnd_mode
);
111 MPFR_CHANGE_SIGN (dest
);
112 mpfr_div_2ui (dest
, dest
, 2, rnd_mode
);
117 inexact
= mpfr_const_pi (dest
, rnd_mode
);
118 mpfr_div_2ui (dest
, dest
, 2, rnd_mode
);
122 else /* +/- 3*PI/4: Ugly since we have to round properly */
125 MPFR_ZIV_DECL (loop2
);
126 mp_prec_t prec2
= MPFR_PREC (dest
) + BITS_PER_MP_LIMB
;
128 mpfr_init2 (tmp2
, prec2
);
129 MPFR_ZIV_INIT (loop2
, prec2
);
132 mpfr_const_pi (tmp2
, GMP_RNDN
);
133 mpfr_mul_ui (tmp2
, tmp2
, 3, GMP_RNDN
); /* Error <= 2 */
134 mpfr_div_2ui (tmp2
, tmp2
, 2, GMP_RNDN
);
135 if (mpfr_round_p (MPFR_MANT (tmp2
), MPFR_LIMB_SIZE (tmp2
),
136 MPFR_PREC (tmp2
) - 2,
137 MPFR_PREC (dest
) + (rnd_mode
== GMP_RNDN
)))
139 MPFR_ZIV_NEXT (loop2
, prec2
);
140 mpfr_set_prec (tmp2
, prec2
);
142 MPFR_ZIV_FREE (loop2
);
144 MPFR_CHANGE_SIGN (tmp2
);
145 inexact
= mpfr_set (dest
, tmp2
, rnd_mode
);
150 MPFR_ASSERTD (MPFR_IS_INF (x
));
157 /* When x=1, atan2(y,x) = atan(y). FIXME: more generally, if x is a power
158 of two, we could call directly atan(y/x) since y/x is exact. */
159 if (mpfr_cmp_ui (x
, 1) == 0)
160 return mpfr_atan (dest
, y
, rnd_mode
);
162 MPFR_SAVE_EXPO_MARK (expo
);
164 /* Set up initial prec */
165 prec
= MPFR_PREC (dest
) + 3 + MPFR_INT_CEIL_LOG2 (MPFR_PREC (dest
));
166 mpfr_init2 (tmp
, prec
);
168 MPFR_ZIV_INIT (loop
, prec
);
170 /* use atan2(y,x) = atan(y/x) */
174 MPFR_BLOCK_DECL (flags
);
176 MPFR_BLOCK (flags
, div_inex
= mpfr_div (tmp
, y
, x
, GMP_RNDN
));
179 /* Result is exact. */
180 inexact
= mpfr_atan (dest
, tmp
, rnd_mode
);
184 /* Error <= ulp (tmp) except in case of underflow or overflow. */
186 /* If the division underflowed, since |atan(z)/z| < 1, we have
188 if (MPFR_UNDERFLOW (flags
))
192 /* In the case GMP_RNDN with 2^(emin-2) < |y/x| < 2^(emin-1):
193 The smallest significand value S > 1 of |y/x| is:
194 * 1 / (1 - 2^(-px)) if py <= px,
195 * (1 - 2^(-px) + 2^(-py)) / (1 - 2^(-px)) if py >= px.
196 Therefore S - 1 > 2^(-pz), where pz = max(px,py). We have:
197 atan(|y/x|) > atan(z), where z = 2^(emin-2) * (1 + 2^(-pz)).
199 > 2^(emin-2) * (1 + 2^(-pz) - 2^(2 emin - 5))
200 Assuming pz <= -2 emin + 5, we can round away from zero
201 (this is what mpfr_underflow always does on GMP_RNDN).
202 In the case GMP_RNDN with |y/x| <= 2^(emin-2), we round
203 towards zero, as |atan(z)/z| < 1. */
204 MPFR_ASSERTN (MPFR_PREC_MAX
<=
205 2 * (mpfr_uexp_t
) - MPFR_EMIN_MIN
+ 5);
206 if (rnd_mode
== GMP_RNDN
&& MPFR_IS_ZERO (tmp
))
208 sign
= MPFR_SIGN (tmp
);
210 MPFR_SAVE_EXPO_FREE (expo
);
211 return mpfr_underflow (dest
, rnd_mode
, sign
);
214 mpfr_atan (tmp
, tmp
, GMP_RNDN
); /* Error <= 2*ulp (tmp) since
215 abs(D(arctan)) <= 1 */
216 /* TODO: check that the error bound is correct in case of overflow. */
217 /* FIXME: Error <= ulp(tmp) ? */
218 if (MPFR_LIKELY (MPFR_CAN_ROUND (tmp
, prec
- 2, MPFR_PREC (dest
),
221 MPFR_ZIV_NEXT (loop
, prec
);
222 mpfr_set_prec (tmp
, prec
);
225 /* Use sign(y)*(PI - atan (|y/x|)) */
227 mpfr_init2 (pi
, prec
);
230 mpfr_div (tmp
, y
, x
, GMP_RNDN
); /* Error <= ulp (tmp) */
231 /* If tmp is 0, we have |y/x| <= 2^(-emin-2), thus
232 atan|y/x| < 2^(-emin-2). */
233 MPFR_SET_POS (tmp
); /* no error */
234 mpfr_atan (tmp
, tmp
, GMP_RNDN
); /* Error <= 2*ulp (tmp) since
235 abs(D(arctan)) <= 1 */
236 mpfr_const_pi (pi
, GMP_RNDN
); /* Error <= ulp(pi) /2 */
237 e
= MPFR_NOTZERO(tmp
) ? MPFR_GET_EXP (tmp
) : __gmpfr_emin
- 1;
238 mpfr_sub (tmp
, pi
, tmp
, GMP_RNDN
); /* see above */
240 MPFR_CHANGE_SIGN (tmp
);
241 /* Error(tmp) <= (1/2+2^(EXP(pi)-EXP(tmp)-1)+2^(e-EXP(tmp)+1))*ulp
242 <= 2^(MAX (MAX (EXP(PI)-EXP(tmp)-1, e-EXP(tmp)+1),
244 e
= MAX (MAX (MPFR_GET_EXP (pi
)-MPFR_GET_EXP (tmp
) - 1,
245 e
- MPFR_GET_EXP (tmp
) + 1), -1) + 2;
246 if (MPFR_LIKELY (MPFR_CAN_ROUND (tmp
, prec
- e
, MPFR_PREC (dest
),
249 MPFR_ZIV_NEXT (loop
, prec
);
250 mpfr_set_prec (tmp
, prec
);
251 mpfr_set_prec (pi
, prec
);
255 inexact
= mpfr_set (dest
, tmp
, rnd_mode
);
258 MPFR_ZIV_FREE (loop
);
260 MPFR_SAVE_EXPO_FREE (expo
);
261 return mpfr_check_range (dest
, inexact
, rnd_mode
);