1 /* mpfr_atan2 -- arc-tan 2 of a floating-point number
2
3 Copyright 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramel projects, INRIA.
5
6 This file is part of the GNU MPFR Library.
7
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.
12
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.
17
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. */
22
23 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
25
26 static int
pi_div_2ui(mpfr_ptr dest,int i,int neg,mpfr_rnd_t rnd_mode)27 pi_div_2ui (mpfr_ptr dest, int i, int neg, mpfr_rnd_t rnd_mode)
28 {
29 int inexact;
30 MPFR_SAVE_EXPO_DECL (expo);
31
32 MPFR_SAVE_EXPO_MARK (expo);
33 if (neg) /* -PI/2^i */
34 {
35 inexact = - mpfr_const_pi (dest, MPFR_INVERT_RND (rnd_mode));
36 MPFR_CHANGE_SIGN (dest);
37 }
38 else /* PI/2^i */
39 {
40 inexact = mpfr_const_pi (dest, rnd_mode);
41 }
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);
45 }
46
47 int
mpfr_atan2(mpfr_ptr dest,mpfr_srcptr y,mpfr_srcptr x,mpfr_rnd_t rnd_mode)48 mpfr_atan2 (mpfr_ptr dest, mpfr_srcptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
49 {
50 mpfr_t tmp, pi;
51 int inexact;
52 mpfr_prec_t prec;
53 mpfr_exp_t e;
54 MPFR_SAVE_EXPO_DECL (expo);
55 MPFR_ZIV_DECL (loop);
56
57 MPFR_LOG_FUNC
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));
63
64 /* Special cases */
65 if (MPFR_ARE_SINGULAR (x, y))
66 {
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.
81 */
82 if (MPFR_IS_NAN (x) || MPFR_IS_NAN (y))
83 {
84 MPFR_SET_NAN (dest);
85 MPFR_RET_NAN;
86 }
87 if (MPFR_IS_ZERO (y))
88 {
89 if (MPFR_IS_NEG (x)) /* +/- PI */
90 {
91 set_pi:
92 if (MPFR_IS_NEG (y))
93 {
94 inexact = mpfr_const_pi (dest, MPFR_INVERT_RND (rnd_mode));
95 MPFR_CHANGE_SIGN (dest);
96 return -inexact;
97 }
98 else
99 return mpfr_const_pi (dest, rnd_mode);
100 }
101 else /* +/- 0 */
102 {
103 set_zero:
104 MPFR_SET_ZERO (dest);
105 MPFR_SET_SAME_SIGN (dest, y);
106 return 0;
107 }
108 }
109 if (MPFR_IS_ZERO (x))
110 {
111 return pi_div_2ui (dest, 1, MPFR_IS_NEG (y), rnd_mode);
112 }
113 if (MPFR_IS_INF (y))
114 {
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 */
120 {
121 mpfr_t tmp2;
122 MPFR_ZIV_DECL (loop2);
123 mpfr_prec_t prec2 = MPFR_PREC (dest) + 10;
124
125 MPFR_SAVE_EXPO_MARK (expo);
126 mpfr_init2 (tmp2, prec2);
127 MPFR_ZIV_INIT (loop2, prec2);
128 for (;;)
129 {
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)))
136 break;
137 MPFR_ZIV_NEXT (loop2, prec2);
138 mpfr_set_prec (tmp2, prec2);
139 }
140 MPFR_ZIV_FREE (loop2);
141 if (MPFR_IS_NEG (y))
142 MPFR_CHANGE_SIGN (tmp2);
143 inexact = mpfr_set (dest, tmp2, rnd_mode);
144 mpfr_clear (tmp2);
145 MPFR_SAVE_EXPO_FREE (expo);
146 return mpfr_check_range (dest, inexact, rnd_mode);
147 }
148 }
149 MPFR_ASSERTD (MPFR_IS_INF (x));
150 if (MPFR_IS_NEG (x))
151 goto set_pi;
152 else
153 goto set_zero;
154 }
155
156 /* When x is a power of two, we call directly atan(y/x) since y/x is
157 exact. */
158 if (MPFR_UNLIKELY (MPFR_IS_POWER_OF_2 (x)))
159 {
160 int r;
161 mpfr_t yoverx;
162 unsigned int saved_flags = __gmpfr_flags;
163
164 mpfr_init2 (yoverx, MPFR_PREC (y));
165 if (MPFR_LIKELY (mpfr_div_2si (yoverx, y, MPFR_GET_EXP (x) - 1,
166 MPFR_RNDN) == 0))
167 {
168 /* Here the flags have not changed due to mpfr_div_2si. */
169 r = mpfr_atan (dest, yoverx, rnd_mode);
170 mpfr_clear (yoverx);
171 return r;
172 }
173 else
174 {
175 /* Division is inexact because of a small exponent range */
176 mpfr_clear (yoverx);
177 __gmpfr_flags = saved_flags;
178 }
179 }
180
181 MPFR_SAVE_EXPO_MARK (expo);
182
183 /* Set up initial prec */
184 prec = MPFR_PREC (dest) + 3 + MPFR_INT_CEIL_LOG2 (MPFR_PREC (dest));
185 mpfr_init2 (tmp, prec);
186
187 MPFR_ZIV_INIT (loop, prec);
188 if (MPFR_IS_POS (x))
189 /* use atan2(y,x) = atan(y/x) */
190 for (;;)
191 {
192 int div_inex;
193 MPFR_BLOCK_DECL (flags);
194
195 MPFR_BLOCK (flags, div_inex = mpfr_div (tmp, y, x, MPFR_RNDN));
196 if (div_inex == 0)
197 {
198 /* Result is exact. */
199 inexact = mpfr_atan (dest, tmp, rnd_mode);
200 goto end;
201 }
202
203 /* Error <= ulp (tmp) except in case of underflow or overflow. */
204
205 /* If the division underflowed, since |atan(z)/z| < 1, we have
206 an underflow. */
207 if (MPFR_UNDERFLOW (flags))
208 {
209 int sign;
210
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)).
217 > z - z^3 / 3.
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);
228 mpfr_clear (tmp);
229 MPFR_SAVE_EXPO_FREE (expo);
230 return mpfr_underflow (dest, rnd_mode, sign);
231 }
232
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),
238 rnd_mode)))
239 break;
240 MPFR_ZIV_NEXT (loop, prec);
241 mpfr_set_prec (tmp, prec);
242 }
243 else /* x < 0 */
244 /* Use sign(y)*(PI - atan (|y/x|)) */
245 {
246 mpfr_init2 (pi, prec);
247 for (;;)
248 {
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 */
258 if (MPFR_IS_NEG (y))
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),
262 -1)+2)*ulp(tmp) */
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),
266 rnd_mode)))
267 break;
268 MPFR_ZIV_NEXT (loop, prec);
269 mpfr_set_prec (tmp, prec);
270 mpfr_set_prec (pi, prec);
271 }
272 mpfr_clear (pi);
273 }
274 inexact = mpfr_set (dest, tmp, rnd_mode);
275
276 end:
277 MPFR_ZIV_FREE (loop);
278 mpfr_clear (tmp);
279 MPFR_SAVE_EXPO_FREE (expo);
280 return mpfr_check_range (dest, inexact, rnd_mode);
281 }
282