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 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[%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