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 5Copyright 2021-2023 Free Software Foundation, Inc. 6Contributed by the AriC and Caramba projects, INRIA. 7 8This file is part of the GNU MPFR Library. 9 10The GNU MPFR Library is free software; you can redistribute it and/or modify 11it under the terms of the GNU Lesser General Public License as published by 12the Free Software Foundation; either version 3 of the License, or (at your 13option) any later version. 14 15The GNU MPFR Library is distributed in the hope that it will be useful, but 16WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 17or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 18License for more details. 19 20You should have received a copy of the GNU Lesser General Public License 21along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 22https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2351 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 */ 31static int 32mpfr_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 */ 49static int 50mpfr_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) */ 72static int 73mpfr_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) */ 95static int 96mpfr_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 */ 119static int 120mpfr_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 */ 146static int 147mpfr_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) */ 163int 164mpfr_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 357int 358mpfr_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} 362