1 /* mpfr_atan2u -- atan2u(y,x,u) = atan(|y/x|)*u/(2*pi) for x > 0
2 atan2u(y,x,u) = 1-atan(|y/x|)*u/(2*pi) for x < 0
3 mpfr_atan2pi -- atan2pi(x) = atan2u(u=2)
4
5 Copyright 2021-2023 Free Software Foundation, Inc.
6 Contributed by the AriC and Caramba projects, INRIA.
7
8 This file is part of the GNU MPFR Library.
9
10 The GNU MPFR Library is free software; you can redistribute it and/or modify
11 it under the terms of the GNU Lesser General Public License as published by
12 the Free Software Foundation; either version 3 of the License, or (at your
13 option) any later version.
14
15 The GNU MPFR Library is distributed in the hope that it will be useful, but
16 WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
17 or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public
18 License for more details.
19
20 You should have received a copy of the GNU Lesser General Public License
21 along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see
22 https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc.,
23 51 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */
24
25 #define MPFR_NEED_LONGLONG_H
26 #include "mpfr-impl.h"
27
28 #define ULSIZE (sizeof (unsigned long) * CHAR_BIT)
29
30 /* z <- s*u*2^e, with e between -3 and -1 */
31 static int
32 mpfr_atan2u_aux1 (mpfr_ptr z, unsigned long u, int e, int s,
33 mpfr_rnd_t rnd_mode)
34 {
35 if (s > 0)
36 return mpfr_set_ui_2exp (z, u, e, rnd_mode);
37 else
38 {
39 int inex;
40
41 rnd_mode = MPFR_INVERT_RND (rnd_mode);
42 inex = mpfr_set_ui_2exp (z, u, e, rnd_mode);
43 MPFR_CHANGE_SIGN (z);
44 return -inex;
45 }
46 }
47
48 /* z <- s*3*u*2^e, with e between -3 and -1 */
49 static int
50 mpfr_atan2u_aux2 (mpfr_ptr z, unsigned long u, int e, int s,
51 mpfr_rnd_t rnd_mode)
52 {
53 int inex;
54 mpfr_t t;
55 MPFR_SAVE_EXPO_DECL (expo);
56
57 MPFR_SAVE_EXPO_MARK (expo);
58 mpfr_init2 (t, ULSIZE + 2);
59 inex = mpfr_set_ui_2exp (t, u, e, MPFR_RNDZ); /* exact */
60 MPFR_ASSERTD (inex == 0);
61 inex = mpfr_mul_ui (t, t, 3, MPFR_RNDZ); /* exact */
62 MPFR_ASSERTD (inex == 0);
63 if (s < 0)
64 MPFR_CHANGE_SIGN (t);
65 inex = mpfr_set (z, t, rnd_mode);
66 mpfr_clear (t);
67 MPFR_SAVE_EXPO_FREE (expo);
68 return mpfr_check_range (z, inex, rnd_mode);
69 }
70
71 /* return round(sign(s)*(u/2-eps),rnd_mode), where eps < 1/2*ulp(u/2) */
72 static int
73 mpfr_atan2u_aux3 (mpfr_ptr z, unsigned long u, int s, mpfr_rnd_t rnd_mode)
74 {
75 mpfr_t t;
76 mpfr_prec_t prec;
77 int inex;
78
79 prec = (MPFR_PREC(z) + 2 > ULSIZE) ? MPFR_PREC(z) + 2 : ULSIZE;
80 /* prec >= PREC(z)+2 and prec >= ULSIZE */
81 mpfr_init2 (t, prec);
82 mpfr_set_ui_2exp (t, u, -1, MPFR_RNDN); /* exact since prec >= ULSIZE */
83 mpfr_nextbelow (t);
84 /* u/2 - 1/4*ulp_p(u/2) <= t <= u/2, where p = PREC(z),
85 which ensures round_p(t) = round_p(u/2-eps) */
86 if (s < 0)
87 mpfr_neg (t, t, MPFR_RNDN);
88 inex = mpfr_set (z, t, rnd_mode);
89 mpfr_clear (t);
90 return inex;
91 }
92
93 /* return round(sign(y)*(u/4-sign(x)*eps),rnd_mode),
94 where eps < 1/2*ulp(u/4) */
95 static int
96 mpfr_atan2u_aux4 (mpfr_ptr z, mpfr_srcptr y, mpfr_srcptr x, unsigned long u,
97 mpfr_rnd_t rnd_mode)
98 {
99 mpfr_t t;
100 mpfr_prec_t prec;
101 int inex;
102
103 prec = (MPFR_PREC(z) > ULSIZE) ? MPFR_PREC(z) + 2 : ULSIZE + 2;
104 /* prec >= PREC(z)+2 and prec >= ULSIZE + 2 */
105 mpfr_init2 (t, prec);
106 mpfr_set_ui_2exp (t, u, -2, MPFR_RNDN); /* exact */
107 if (MPFR_SIGN(x) > 0)
108 mpfr_nextbelow (t);
109 else
110 mpfr_nextabove (t);
111 if (MPFR_SIGN(y) < 0)
112 mpfr_neg (t, t, MPFR_RNDN);
113 inex = mpfr_set (z, t, rnd_mode);
114 mpfr_clear (t);
115 return inex;
116 }
117
118 /* deal with underflow case, i.e., when y/x rounds to zero */
119 static int
120 mpfr_atan2u_underflow (mpfr_ptr z, mpfr_srcptr y, mpfr_srcptr x,
121 unsigned long u, mpfr_rnd_t rnd_mode)
122 {
123 mpfr_exp_t e = MPFR_GET_EXP(y) - MPFR_GET_EXP(x) + (ULSIZE - 1);
124 /* Detect underflow: since |atan(|y/x|)| < |y/x| for |y/x| < 1,
125 |atan2u(y,x,u)| < |y/x|*u/(2*pi) < 2^(ULSIZE-2)*|y/x|
126 < 2^(EXP(y)-EXP(x)+(ULSIZE-1)).
127 For x > 0, we have underflow when
128 EXP(y)-EXP(x)+(ULSIZE-1) < emin.
129 For x < 0, we have underflow when
130 EXP(y)-EXP(x)+(ULSIZE-1) < EXP(u/2)-prec. */
131 if (MPFR_IS_POS(x))
132 {
133 MPFR_ASSERTN(e < __gmpfr_emin);
134 return mpfr_underflow (z,
135 (rnd_mode == MPFR_RNDN) ? MPFR_RNDZ : rnd_mode, MPFR_SIGN(y));
136 }
137 else
138 {
139 MPFR_ASSERTD (MPFR_IS_NEG(x));
140 MPFR_ASSERTN(e < (ULSIZE - 1) - MPFR_PREC(z));
141 return mpfr_atan2u_aux3 (z, u, MPFR_SIGN(y), rnd_mode);
142 }
143 }
144
145 /* deal with overflow case, i.e., when y/x rounds to infinity */
146 static int
147 mpfr_atan2u_overflow (mpfr_ptr z, mpfr_srcptr y, mpfr_srcptr x,
148 unsigned long u, mpfr_rnd_t rnd_mode)
149 {
150 /* When t goes to +Inf, pi/2 - 1/t < atan(t) < pi/2,
151 thus u/4 - u/(2*pi*t) < atanu(t) < u/4.
152 As soon as u/(2*pi*t) < 1/2*ulp(u/4), the result is either u/4
153 or the number just below.
154 Here t = y/x, thus 1/t <= x/y < 2^(EXP(x)-EXP(y)+1),
155 and u/(2*pi*t) < 2^(EXP(x)-EXP(y)+(ULSIZE-2)). */
156 mpfr_exp_t e = MPFR_GET_EXP(x) - MPFR_GET_EXP(y) + (ULSIZE - 2);
157 mpfr_exp_t ulpz = (ULSIZE - 2) - MPFR_PREC(z); /* ulp(u/4) <= 2^ulpz */
158 MPFR_ASSERTN (e < ulpz - 1);
159 return mpfr_atan2u_aux4 (z, y, x, u, rnd_mode);
160 }
161
162 /* put in z the correctly rounded value of atan2y(y,x,u) */
163 int
164 mpfr_atan2u (mpfr_ptr z, mpfr_srcptr y, mpfr_srcptr x, unsigned long u,
165 mpfr_rnd_t rnd_mode)
166 {
167 mpfr_t tmp;
168 mpfr_prec_t prec;
169 int inex;
170 MPFR_SAVE_EXPO_DECL (expo);
171 MPFR_ZIV_DECL (loop);
172
173 MPFR_LOG_FUNC
174 (("y[%Pd]=%.*Rg x[%Pd]=%.*Rg u=%lu rnd=%d",
175 mpfr_get_prec(y), mpfr_log_prec, y,
176 mpfr_get_prec(x), mpfr_log_prec, x,
177 u, rnd_mode),
178 ("z[%Pd]=%.*Rg inexact=%d", mpfr_get_prec (z), mpfr_log_prec, z,
179 inex));
180
181 /* Special cases */
182 if (MPFR_ARE_SINGULAR (x, y))
183 {
184 if (MPFR_IS_NAN (x) || MPFR_IS_NAN (y))
185 {
186 MPFR_SET_NAN (z);
187 MPFR_RET_NAN;
188 }
189 /* remains (y=Inf,x=Inf), (Inf,0), (Inf,regular), (0,Inf), (0,0),
190 (0,regular), (regular,Inf), (regular,0) */
191 if (MPFR_IS_INF (x))
192 {
193 /* cases (y=Inf,x=Inf), (0,Inf), (regular,Inf) */
194 if (MPFR_IS_INF (y))
195 {
196 if (MPFR_IS_NEG (x))
197 {
198 /* atan2Pi(+/-Inf,-Inf) = +/-3/4,
199 atan2u(+/-Inf,-Inf,u) = +/-3u/8 */
200 return mpfr_atan2u_aux2 (z, u, -3, MPFR_SIGN (y), rnd_mode);
201 }
202 else /* x > 0 */
203 {
204 /* atan2Pi(+/-Inf,-Inf) = +/-1/4,
205 atan2u(+/-Inf,-Inf,u) = +/-u/8 */
206 return mpfr_atan2u_aux1 (z, u, -3, MPFR_SIGN (y), rnd_mode);
207 }
208 }
209 /* remains (y,x) = (0,Inf) and (regular,Inf),
210 in which cases the IEEE 754-2019 rules for (y=0,x<>0)
211 and for (y<>0,x Inf) coincide */
212 if (MPFR_IS_NEG (x))
213 /* y>0: atan2Pi(+/-y,-Inf) = +/-1,
214 atan2u(+/-y,-Inf,u) = +/-u/2 */
215 return mpfr_atan2u_aux1 (z, u, -1, MPFR_SIGN (y), rnd_mode);
216 else
217 {
218 /* y>0: atan2Pi(+/-y,+Inf) = +/-0,
219 atan2u(+/-y,+Inf,u) = +/-0 */
220 MPFR_SET_ZERO (z);
221 MPFR_SET_SAME_SIGN (z, y);
222 MPFR_RET (0);
223 }
224 }
225 /* remains (y=Inf,x=0), (Inf,regular), (0,0), (0,regular), (regular,0) */
226 if (MPFR_IS_INF (y))
227 /* x is zero or regular */
228 /* atan2Pi(+/-Inf, x) = +/-1/2, atan2u(+/-Inf, x, u) = +/-u/4 */
229 return mpfr_atan2u_aux1 (z, u, -2, MPFR_SIGN(y), rnd_mode);
230 /* remains (y=0,x=0), (0,regular), (regular,0) */
231 if (MPFR_IS_ZERO (y))
232 {
233 if (MPFR_IS_NEG (x))
234 {
235 /* atan2Pi(+/-0, x) = +/-1, atan2u(+/-0, x, u) = +/-u/2 */
236 return mpfr_atan2u_aux1 (z, u, -1, MPFR_SIGN(y), rnd_mode);
237 }
238 else /* case x = +0 or x > 0 */
239 {
240 /* atan2Pi(+/-0, x) = +/-0, atan2u(+/-0, x, u) = +/-0 */
241 MPFR_SET_ZERO (z); /* does not modify sign, in case z=y */
242 MPFR_SET_SAME_SIGN (z, y);
243 MPFR_RET (0); /* exact result */
244 }
245 }
246 /* remains (regular,0) */
247 if (MPFR_IS_ZERO (x))
248 {
249 /* y<0: atan2Pi(y,+/-0) = -1/2, thus atan2u(y,+/-0,u) = -u/4 */
250 /* y>0: atan2Pi(y,+/-0) = 1/2, thus atan2u(y,+/-0,u) = u/4 */
251 return mpfr_atan2u_aux1 (z, u, -2, MPFR_SIGN(y), rnd_mode);
252 }
253 /* no special case should remain */
254 MPFR_RET_NEVER_GO_HERE ();
255 }
256
257 /* IEEE 754-2019 says that atan2Pi is odd with respect to y */
258
259 /* now both y and x are regular */
260 if (mpfr_cmpabs (y, x) == 0)
261 {
262 if (MPFR_IS_POS (x))
263 /* atan2u (+/-x,x,u) = +/u/8 for x > 0 */
264 return mpfr_atan2u_aux1 (z, u, -3, MPFR_SIGN(y), rnd_mode);
265 else /* x < 0 */
266 /* atan2u (+/-x,x,u) = -/+3u/8 for x > 0 */
267 return mpfr_atan2u_aux2 (z, u, -3, MPFR_SIGN(y), rnd_mode);
268 }
269
270 if (u == 0) /* return 0 with sign of y for x > 0,
271 and 1 with sign of y for x < 0 */
272 {
273 if (MPFR_SIGN(x) > 0)
274 {
275 MPFR_SET_ZERO (z);
276 MPFR_SET_SAME_SIGN (z, y);
277 MPFR_RET (0);
278 }
279 else
280 return mpfr_set_si (z, MPFR_SIGN(y) > 0 ? 1 : -1, rnd_mode);
281 }
282
283 MPFR_SAVE_EXPO_MARK (expo);
284 prec = MPFR_PREC (z);
285 prec += MPFR_INT_CEIL_LOG2(prec) + 10;
286 mpfr_init2 (tmp, prec);
287 MPFR_ZIV_INIT (loop, prec);
288 for (;;)
289 {
290 mpfr_exp_t expt, e;
291 /* atan2u(-y,x,u) = -atan(y,x,u)
292 atan2u(y,x,u) = atan(|y/x|)*u/(2*pi) for x > 0
293 atan2u(y,x,u) = u/2-atan(|y/x|)*u/(2*pi) for x < 0
294 atan2pi atan2u
295 First quadrant: [0,1/2] [0,u/4]
296 Second quadrant: [1/2,1] [u/4,u/2]
297 Third quadrant: [-1,-1/2] [-u/2,-u/4]
298 Fourth quadrant: [-1/2,0] [-u/4,0] */
299 inex = mpfr_div (tmp, y, x, MPFR_RNDN);
300 if (MPFR_IS_ZERO(tmp))
301 {
302 mpfr_clear (tmp);
303 MPFR_SAVE_EXPO_FREE (expo);
304 return mpfr_atan2u_underflow (z, y, x, u, rnd_mode);
305 }
306 if (MPFR_IS_INF(tmp))
307 {
308 mpfr_clear (tmp);
309 MPFR_SAVE_EXPO_FREE (expo);
310 return mpfr_atan2u_overflow (z, y, x, u, rnd_mode);
311 }
312 MPFR_SET_SIGN (tmp, 1);
313 expt = MPFR_GET_EXP(tmp);
314 /* |tmp - |y/x|| <= e1 := 1/2*ulp(tmp) = 2^(expt-prec-1) */
315 inex |= mpfr_atanu (tmp, tmp, u, MPFR_RNDN);
316 /* the derivative of atan(t) is 1/(1+t^2), thus that of atanu(t) is
317 u/(1+t^2)/(2*pi), and if tmp2 is the new value of tmp, we have
318 |tmp2 - atanu|y/x|| <= 1/2*ulp(tmp2) + e1*u/(1+tmp^2)/4 */
319 e = (expt < 1) ? 0 : expt - 1;
320 /* max(1,|tmp|) >= 2^e thus 1/(1+tmp^2) <= 2^(-2*e) */
321 e = expt - 2 * e + MPFR_INT_CEIL_LOG2(u) - 2;
322 /* now e1*u/(1+tmp^2)/4 <= 2^(e-prec-1) */
323 /* |tmp2 - atanu(y/x)| <= 1/2*ulp(tmp2) + 2^(e-prec-1) */
324 e = (e < MPFR_GET_EXP(tmp)) ? MPFR_GET_EXP(tmp) : e;
325 /* |tmp2 - atanu(y/x)| <= 2^(e-prec) */
326 if (MPFR_SIGN (x) < 0)
327 { /* compute u/2-tmp */
328 mpfr_mul_2ui (tmp, tmp, 1, MPFR_RNDN); /* error <= 2^(e+1-prec) */
329 mpfr_ui_sub (tmp, u, tmp, MPFR_RNDN);
330 expt = MPFR_GET_EXP(tmp);
331 /* error <= 2^(expt-prec-1) + 2^(e+1-prec) */
332 e = (expt - 1 > e + 1) ? expt - 1 : e + 1;
333 /* error <= 2^(e+1-prec) */
334 mpfr_div_2ui (tmp, tmp, 1, MPFR_RNDN);
335 /* error <= 2^(e-prec) */
336 }
337 /* both with x>0 and x<0, we have error <= 2^(e-prec),
338 now we want error <= 2^(expt-prec+err)
339 thus err = e-expt */
340 e -= MPFR_GET_EXP(tmp);
341 if (MPFR_SIGN(y) < 0) /* atan2u is odd with respect to y */
342 MPFR_CHANGE_SIGN(tmp);
343 if (MPFR_LIKELY (MPFR_CAN_ROUND (tmp, prec - e,
344 MPFR_PREC (z), rnd_mode)))
345 break;
346 MPFR_ZIV_NEXT (loop, prec);
347 mpfr_set_prec (tmp, prec);
348 }
349 MPFR_ZIV_FREE (loop);
350 inex = mpfr_set (z, tmp, rnd_mode);
351 mpfr_clear (tmp);
352 MPFR_SAVE_EXPO_FREE (expo);
353
354 return mpfr_check_range (z, inex, rnd_mode);
355 }
356
357 int
358 mpfr_atan2pi (mpfr_ptr z, mpfr_srcptr y, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
359 {
360 return mpfr_atan2u (z, y, x, 2, rnd_mode);
361 }