1 /* mpfr_atan2 -- arc-tan 2 of a floating-point number
2
3 Copyright 2005-2023 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramba 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 https://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
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
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[%Pd]=%.*Rg x[%Pd]=%.*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[%Pd]=%.*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_CAN_ROUND (tmp2, MPFR_PREC (tmp2) - 2,
134 MPFR_PREC (dest), rnd_mode))
135 break;
136 MPFR_ZIV_NEXT (loop2, prec2);
137 mpfr_set_prec (tmp2, prec2);
138 }
139 MPFR_ZIV_FREE (loop2);
140 if (MPFR_IS_NEG (y))
141 MPFR_CHANGE_SIGN (tmp2);
142 inexact = mpfr_set (dest, tmp2, rnd_mode);
143 mpfr_clear (tmp2);
144 MPFR_SAVE_EXPO_FREE (expo);
145 return mpfr_check_range (dest, inexact, rnd_mode);
146 }
147 }
148 MPFR_ASSERTD (MPFR_IS_INF (x));
149 if (MPFR_IS_NEG (x))
150 goto set_pi;
151 else
152 goto set_zero;
153 }
154
155 /* When x is a power of two, we call directly atan(y/x) since y/x is
156 exact. */
157 if (MPFR_UNLIKELY (MPFR_IS_POS (x) && mpfr_powerof2_raw (x)))
158 {
159 int r;
160 mpfr_t yoverx;
161 mpfr_flags_t saved_flags = __gmpfr_flags;
162
163 mpfr_init2 (yoverx, MPFR_PREC (y));
164 if (MPFR_LIKELY (mpfr_div_2si (yoverx, y, MPFR_GET_EXP (x) - 1,
165 MPFR_RNDN) == 0))
166 {
167 /* Here the flags have not changed due to mpfr_div_2si. */
168 r = mpfr_atan (dest, yoverx, rnd_mode);
169 mpfr_clear (yoverx);
170 return r;
171 }
172 else
173 {
174 /* Division is inexact because of a small exponent range */
175 mpfr_clear (yoverx);
176 __gmpfr_flags = saved_flags;
177 }
178 }
179
180 MPFR_SAVE_EXPO_MARK (expo);
181
182 /* Set up initial prec */
183 prec = MPFR_PREC (dest) + 3 + MPFR_INT_CEIL_LOG2 (MPFR_PREC (dest));
184 mpfr_init2 (tmp, prec);
185
186 MPFR_ZIV_INIT (loop, prec);
187 if (MPFR_IS_POS (x))
188 /* use atan2(y,x) = atan(y/x) */
189 for (;;)
190 {
191 int div_inex;
192 MPFR_BLOCK_DECL (flags);
193
194 MPFR_BLOCK (flags, div_inex = mpfr_div (tmp, y, x, MPFR_RNDN));
195 if (div_inex == 0)
196 {
197 /* Result is exact. */
198 inexact = mpfr_atan (dest, tmp, rnd_mode);
199 goto end;
200 }
201
202 /* Error <= ulp (tmp) except in case of underflow or overflow. */
203
204 /* If the division underflowed, since |atan(z)/z| < 1, we have
205 an underflow. */
206 if (MPFR_UNDERFLOW (flags))
207 {
208 int sign;
209
210 /* In the case MPFR_RNDN with 2^(emin-2) < |y/x| < 2^(emin-1):
211 The smallest significand value S > 1 of |y/x| is:
212 * 1 / (1 - 2^(-px)) if py <= px,
213 * (1 - 2^(-px) + 2^(-py)) / (1 - 2^(-px)) if py >= px.
214 Therefore S - 1 > 2^(-pz), where pz = max(px,py). We have:
215 atan(|y/x|) > atan(z), where z = 2^(emin-2) * (1 + 2^(-pz)).
216 > z - z^3 / 3.
217 > 2^(emin-2) * (1 + 2^(-pz) - 2^(2 emin - 5))
218 Assuming pz <= -2 emin + 5, we can round away from zero
219 (this is what mpfr_underflow always does on MPFR_RNDN).
220 In the case MPFR_RNDN with |y/x| <= 2^(emin-2), we round
221 toward zero, as |atan(z)/z| < 1. */
222 MPFR_ASSERTN (MPFR_PREC_MAX <=
223 2 * (mpfr_uexp_t) - MPFR_EMIN_MIN + 5);
224 if (rnd_mode == MPFR_RNDN && MPFR_IS_ZERO (tmp))
225 rnd_mode = MPFR_RNDZ;
226 sign = MPFR_SIGN (tmp);
227 mpfr_clear (tmp);
228 MPFR_SAVE_EXPO_FREE (expo);
229 return mpfr_underflow (dest, rnd_mode, sign);
230 }
231
232 mpfr_atan (tmp, tmp, MPFR_RNDN); /* Error <= 2*ulp (tmp) since
233 abs(D(arctan)) <= 1 */
234 /* TODO: check that the error bound is correct in case of overflow. */
235 /* FIXME: Error <= ulp(tmp) ? */
236 if (MPFR_LIKELY (MPFR_CAN_ROUND (tmp, prec - 2, MPFR_PREC (dest),
237 rnd_mode)))
238 break;
239 MPFR_ZIV_NEXT (loop, prec);
240 mpfr_set_prec (tmp, prec);
241 }
242 else /* x < 0 */
243 /* Use sign(y)*(PI - atan (|y/x|)) */
244 {
245 mpfr_init2 (pi, prec);
246 for (;;)
247 {
248 mpfr_div (tmp, y, x, MPFR_RNDN); /* Error <= ulp (tmp) */
249 /* If tmp is 0, we have |y/x| <= 2^(-emin-2), thus
250 atan|y/x| < 2^(-emin-2). */
251 MPFR_SET_POS (tmp); /* no error */
252 mpfr_atan (tmp, tmp, MPFR_RNDN); /* Error <= 2*ulp (tmp) since
253 abs(D(arctan)) <= 1 */
254 mpfr_const_pi (pi, MPFR_RNDN); /* Error <= ulp(pi) /2 */
255 e = MPFR_NOTZERO(tmp) ? MPFR_GET_EXP (tmp) : __gmpfr_emin - 1;
256 mpfr_sub (tmp, pi, tmp, MPFR_RNDN); /* see above */
257 if (MPFR_IS_NEG (y))
258 MPFR_CHANGE_SIGN (tmp);
259 /* Error(tmp) <= (1/2+2^(EXP(pi)-EXP(tmp)-1)+2^(e-EXP(tmp)+1))*ulp
260 <= 2^(MAX (MAX (EXP(PI)-EXP(tmp)-1, e-EXP(tmp)+1),
261 -1)+2)*ulp(tmp) */
262 e = MAX (MAX (MPFR_GET_EXP (pi)-MPFR_GET_EXP (tmp) - 1,
263 e - MPFR_GET_EXP (tmp) + 1), -1) + 2;
264 if (MPFR_LIKELY (MPFR_CAN_ROUND (tmp, prec - e, MPFR_PREC (dest),
265 rnd_mode)))
266 break;
267 MPFR_ZIV_NEXT (loop, prec);
268 mpfr_set_prec (tmp, prec);
269 mpfr_set_prec (pi, prec);
270 }
271 mpfr_clear (pi);
272 }
273 inexact = mpfr_set (dest, tmp, rnd_mode);
274
275 end:
276 MPFR_ZIV_FREE (loop);
277 mpfr_clear (tmp);
278 MPFR_SAVE_EXPO_FREE (expo);
279 return mpfr_check_range (dest, inexact, rnd_mode);
280 }