1/* mpfr_atan -- arc-tangent of a floating-point number 2 3Copyright 2001-2023 Free Software Foundation, Inc. 4Contributed by the AriC and Caramba projects, INRIA. 5 6This file is part of the GNU MPFR Library. 7 8The GNU MPFR Library is free software; you can redistribute it and/or modify 9it under the terms of the GNU Lesser General Public License as published by 10the Free Software Foundation; either version 3 of the License, or (at your 11option) any later version. 12 13The GNU MPFR Library is distributed in the hope that it will be useful, but 14WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY 15or FITNESS FOR A PARTICULAR PURPOSE. See the GNU Lesser General Public 16License for more details. 17 18You should have received a copy of the GNU Lesser General Public License 19along with the GNU MPFR Library; see the file COPYING.LESSER. If not, see 20https://www.gnu.org/licenses/ or write to the Free Software Foundation, Inc., 2151 Franklin St, Fifth Floor, Boston, MA 02110-1301, USA. */ 22 23#define MPFR_NEED_LONGLONG_H 24#include "mpfr-impl.h" 25 26#if GMP_NUMB_BITS == 64 27/* for each pair (r,p), we store a 192-bit approximation of atan(x)/x for 28 x=p/2^r, with lowest limb first. 29 Sage code: 30 for p in range(1,2^ceil(r/2)): 31 x=p/2^r 32 l=floor(2^192*n(atan(x)/x, 300)).digits(2^64) 33 print ("{0x%x, 0x%x, 0x%x}, /"+"* (%d,%d) *"+"/") % (l[0],l[1],l[2],r,p) 34*/ 35static const mp_limb_t atan_table[][3] = { 36 {0x6e141587261cdf00, 0x6fe445ecbc3a8d03, 0xed63382b0dda7b45}, /* (1,1) */ 37 {0xaa7fa90388b3836b, 0x6dc79ef5f7a217e5, 0xfadbafc96406eb15}, /* (2,1) */ 38 {0x319c12cf59d4b2dc, 0xcb2792dc0e2e0d51, 0xffaaddb967ef4e36}, /* (4,1) */ 39 {0x8b3957d95d9ad922, 0xc897989f3e888ef7, 0xfeadd4d5617b6e32}, /* (4,2) */ 40 {0xc4e6abc8af62e439, 0x4eb9bf602625f0b4, 0xfd0fcdd343cac19b}, /* (4,3) */ 41 {0x7c18baeb9bc95789, 0xb12afb6b6d4f7e16, 0xffffaaaaddddb94b}, /* (8,1) */ 42 {0x6856a0171a2f001a, 0x62351fbbe60af47, 0xfffeaaadddd4b968}, /* (8,2) */ 43 {0x69164c094f49da06, 0xd517294f7373d07a, 0xfffd001032cb1179}, /* (8,3) */ 44 {0x20ef65c10deef460, 0xe78c564015f76048, 0xfffaaadddb94d5bb}, /* (8,4) */ 45 {0x3ce233aa002f0344, 0x9dd8ea342a65d4cc, 0xfff7ab27a1f32f95}, /* (8,5) */ 46 {0xa37f403c7279c5cb, 0x13ab53a1c8db8497, 0xfff40103192ce74d}, /* (8,6) */ 47 {0xe5a85657103c1aa8, 0xb8409e6c914191d3, 0xffefac8a9c40a26b}, /* (8,7) */ 48 {0x806d0294c0db8816, 0x779d776dda8c6213, 0xffeaaddd4bb12542}, /* (8,8) */ 49 {0x5545d1914ef21478, 0x3aea58d6660f5a12, 0xffe5051f0aebf73a}, /* (8,9) */ 50 {0x6e47a91d015f4133, 0xc085ab6b490b7f02, 0xffdeb2787d4adac1}, /* (8,10) */ 51 {0x4efc1f931f7ec9b3, 0xb7f43cd16195ef4b, 0xffd7b61702b09aad}, /* (8,11) */ 52 {0xd27d1dbf55fed60d, 0xd812c11d7d473e5e, 0xffd0102cb3c1bfbe}, /* (8,12) */ 53 {0xca629e927383fe97, 0x8c61aedf58e42206, 0xffc7c0f05db9d1b6}, /* (8,13) */ 54 {0x4eff0b53d4e905b7, 0x28ac1e800ca31e9d, 0xffbec89d7dddd7e9}, /* (8,14) */ 55 {0xb0a7931deec6fe60, 0xb46feea78588554b, 0xffb527743c8cdd8f} /* (8,15) */ 56 }; 57 58static void 59set_table (mpfr_ptr y, const mp_limb_t x[3]) 60{ 61 mpfr_prec_t p = MPFR_PREC(y); 62 mp_size_t n = MPFR_PREC2LIMBS(p); 63 mpfr_prec_t sh; 64 mp_limb_t *yp = MPFR_MANT(y); 65 66 MPFR_UNSIGNED_MINUS_MODULO (sh, p); 67 MPFR_ASSERTD (n >= 1 && n <= 3); 68 mpn_copyi (yp, x + 3 - n, n); 69 yp[0] &= ~MPFR_LIMB_MASK(sh); 70 MPFR_SET_EXP(y, 0); 71} 72#endif 73 74/* If x = p/2^r, put in y an approximation to atan(x)/x using 2^m terms 75 for the series expansion, with an error of at most 1 ulp. 76 Assumes 0 < x < 1, thus 1 <= p < 2^r. 77 More precisely, p consists of the floor(r/2) bits of the binary expansion 78 of a number 0 < s < 1: 79 * the bit of weight 2^-1 is for r=1, thus p <= 1 80 * the bit of weight 2^-2 is for r=2, thus p <= 1 81 * the two bits of weight 2^-3 and 2^-4 are for r=4, thus p <= 3 82 * more generally p < 2^(r/2). 83 84 If X=x^2, we want 1 - X/3 + X^2/5 - ... + (-1)^k*X^k/(2k+1) + ... 85 86 When we sum terms up to x^k/(2k+1), the denominator Q[0] is 87 3*5*7*...*(2k+1) ~ (2k/e)^k. 88 89 The tab[] array should have at least 3*(m+1) entries. 90*/ 91static void 92mpfr_atan_aux (mpfr_ptr y, mpz_ptr p, unsigned long r, int m, mpz_t *tab) 93{ 94 mpz_t *S, *Q, *ptoj; 95 mp_bitcnt_t n, h, j; /* unsigned type, which is >= unsigned long */ 96 mpfr_exp_t diff, expo; 97 int im, i, k, l, done; 98 mpfr_prec_t mult; 99 mpfr_prec_t accu[MPFR_PREC_BITS], log2_nb_terms[MPFR_PREC_BITS]; 100 mpfr_prec_t precy = MPFR_PREC(y); 101 102 MPFR_ASSERTD (mpz_sgn (p) > 0); 103 MPFR_ASSERTD (m > 0); 104 MPFR_ASSERTD (m <= MPFR_PREC_BITS - 1); 105 106#if GMP_NUMB_BITS == 64 107 /* tabulate values for small precision and small value of r (which are the 108 most expensive to compute) */ 109 if (precy <= 192) 110 { 111 unsigned long u; 112 113 switch (r) 114 { 115 case 1: 116 /* p has 1 bit: necessarily p=1 */ 117 MPFR_ASSERTD(mpz_cmp_ui (p, 1) == 0); 118 set_table (y, atan_table[0]); 119 return; 120 case 2: 121 /* p has 1 bit: necessarily p=1 too */ 122 MPFR_ASSERTD(mpz_cmp_ui (p, 1) == 0); 123 set_table (y, atan_table[1]); 124 return; 125 case 4: 126 /* p has at most 2 bits: 1 <= p <= 3 */ 127 u = mpz_get_ui (p); 128 MPFR_ASSERTD(1 <= u && u <= 3); 129 set_table (y, atan_table[1 + u]); 130 return; 131 case 8: 132 /* p has at most 4 bits: 1 <= p <= 15 */ 133 u = mpz_get_ui (p); 134 MPFR_ASSERTD(1 <= u && u <= 15); 135 set_table (y, atan_table[4 + u]); 136 return; 137 } 138 } 139#endif 140 141 /* Set Tables */ 142 S = tab; /* S */ 143 ptoj = S + 1*(m+1); /* p^2^j Precomputed table */ 144 Q = S + 2*(m+1); /* Product of Odd integer table */ 145 146 /* From p to p^2, and r to 2r */ 147 mpz_mul (p, p, p); 148 MPFR_ASSERTD (2 * r > r); 149 r = 2 * r; 150 151 /* Normalize p */ 152 n = mpz_scan1 (p, 0); 153 if (n > 0) 154 { 155 mpz_tdiv_q_2exp (p, p, n); /* exact */ 156 MPFR_ASSERTD (r > n); 157 r -= n; 158 } 159 160 /* Since |p/2^r| < 1, and p is a non-zero integer, necessarily r > 0. */ 161 MPFR_ASSERTD (mpz_sgn (p) > 0); 162 MPFR_ASSERTD (m > 0); 163 MPFR_ASSERTD (r > 0); 164 165 /* check if p=1 (special case) */ 166 l = 0; 167 /* 168 We compute by binary splitting, with X = x^2 = p/2^r: 169 P(a,b) = p if a+1=b, P(a,c)*P(c,b) otherwise 170 Q(a,b) = (2a+1)*2^r if a+1=b [except Q(0,1)=1], Q(a,c)*Q(c,b) otherwise 171 S(a,b) = p*(2a+1) if a+1=b, Q(c,b)*S(a,c)+Q(a,c)*P(a,c)*S(c,b) otherwise 172 Then atan(x)/x ~ S(0,i)/Q(0,i) for i so that (p/2^r)^i/i is small enough. 173 The factor 2^(r*(b-a)) in Q(a,b) is implicit, thus we have to take it 174 into account when we compute with Q. 175 */ 176 accu[0] = 0; /* accu[k] = Mult[0] + ... + Mult[k], where Mult[j] is the 177 number of bits of the corresponding term S[j]/Q[j] */ 178 if (mpz_cmp_ui (p, 1) != 0) 179 { 180 /* p <> 1: precompute ptoj table */ 181 mpz_set (ptoj[0], p); 182 for (im = 1 ; im <= m ; im ++) 183 mpz_mul (ptoj[im], ptoj[im - 1], ptoj[im - 1]); 184 /* main loop */ 185 n = 1UL << m; 186 MPFR_ASSERTN (n != 0); /* no overflow */ 187 /* the i-th term being X^i/(2i+1) with X=p/2^r, we can stop when 188 p^i/2^(r*i) < 2^(-precy), i.e. r*i > precy + log2(p^i) */ 189 for (i = k = done = 0; (i < n) && (done == 0); i += 2, k ++) 190 { 191 /* initialize both S[k],Q[k] and S[k+1],Q[k+1] */ 192 mpz_set_ui (Q[k+1], 2 * i + 3); /* Q(i+1,i+2) */ 193 mpz_mul_ui (S[k+1], p, 2 * i + 1); /* S(i+1,i+2) */ 194 mpz_mul_2exp (S[k], Q[k+1], r); 195 mpz_sub (S[k], S[k], S[k+1]); /* S(i,i+2) */ 196 mpz_mul_ui (Q[k], Q[k+1], 2 * i + 1); /* Q(i,i+2) */ 197 log2_nb_terms[k] = 1; /* S[k]/Q[k] corresponds to 2 terms */ 198 for (j = (i + 2) >> 1, l = 1; (j & 1) == 0; l ++, j >>= 1, k --) 199 { 200 /* invariant: S[k-1]/Q[k-1] and S[k]/Q[k] correspond 201 to 2^l terms each. We combine them into S[k-1]/Q[k-1] */ 202 MPFR_ASSERTD (k > 0); 203 mpz_mul (S[k], S[k], Q[k-1]); 204 mpz_mul (S[k], S[k], ptoj[l]); 205 mpz_mul (S[k-1], S[k-1], Q[k]); 206 mpz_mul_2exp (S[k-1], S[k-1], r << l); 207 mpz_add (S[k-1], S[k-1], S[k]); 208 mpz_mul (Q[k-1], Q[k-1], Q[k]); 209 log2_nb_terms[k-1] = l + 1; 210 /* now S[k-1]/Q[k-1] corresponds to 2^(l+1) terms */ 211 MPFR_MPZ_SIZEINBASE2(mult, ptoj[l+1]); 212 mult = (r << (l + 1)) - mult - 1; 213 accu[k-1] = (k == 1) ? mult : accu[k-2] + mult; 214 if (accu[k-1] > precy) 215 done = 1; 216 } 217 } 218 } 219 else /* special case p=1: the i-th term being X^i/(2i+1) with X=1/2^r, 220 we can stop when r*i > precy i.e. i > precy/r */ 221 { 222 n = 1UL << m; 223 if (precy / r <= n) 224 n = (precy / r) + 1; 225 MPFR_ASSERTN (n != 0); /* no overflow */ 226 for (i = k = 0; i < n; i += 2, k ++) 227 { 228 mpz_set_ui (Q[k + 1], 2 * i + 3); 229 mpz_mul_2exp (S[k], Q[k+1], r); 230 mpz_sub_ui (S[k], S[k], 1 + 2 * i); 231 mpz_mul_ui (Q[k], Q[k + 1], 1 + 2 * i); 232 log2_nb_terms[k] = 1; /* S[k]/Q[k] corresponds to 2 terms */ 233 for (j = (i + 2) >> 1, l = 1; (j & 1) == 0; l++, j >>= 1, k --) 234 { 235 MPFR_ASSERTD (k > 0); 236 mpz_mul (S[k], S[k], Q[k-1]); 237 mpz_mul (S[k-1], S[k-1], Q[k]); 238 mpz_mul_2exp (S[k-1], S[k-1], r << l); 239 mpz_add (S[k-1], S[k-1], S[k]); 240 mpz_mul (Q[k-1], Q[k-1], Q[k]); 241 log2_nb_terms[k-1] = l + 1; 242 } 243 } 244 } 245 246 /* we need to combine S[0]/Q[0]...S[k-1]/Q[k-1] */ 247 h = 0; /* number of terms accumulated in S[k]/Q[k] */ 248 while (k > 1) 249 { 250 k --; 251 /* combine S[k-1]/Q[k-1] and S[k]/Q[k] */ 252 mpz_mul (S[k], S[k], Q[k-1]); 253 if (mpz_cmp_ui (p, 1) != 0) 254 mpz_mul (S[k], S[k], ptoj[log2_nb_terms[k-1]]); 255 mpz_mul (S[k-1], S[k-1], Q[k]); 256 h += (mp_bitcnt_t) 1 << log2_nb_terms[k]; 257 mpz_mul_2exp (S[k-1], S[k-1], r * h); 258 mpz_add (S[k-1], S[k-1], S[k]); 259 mpz_mul (Q[k-1], Q[k-1], Q[k]); 260 } 261 262 MPFR_MPZ_SIZEINBASE2 (diff, S[0]); 263 diff -= 2 * precy; 264 expo = diff; 265 if (diff >= 0) 266 mpz_tdiv_q_2exp (S[0], S[0], diff); 267 else 268 mpz_mul_2exp (S[0], S[0], -diff); 269 270 MPFR_MPZ_SIZEINBASE2 (diff, Q[0]); 271 diff -= precy; 272 expo -= diff; 273 if (diff >= 0) 274 mpz_tdiv_q_2exp (Q[0], Q[0], diff); 275 else 276 mpz_mul_2exp (Q[0], Q[0], -diff); 277 278 mpz_tdiv_q (S[0], S[0], Q[0]); 279 mpfr_set_z (y, S[0], MPFR_RNDD); 280 /* TODO: Check/prove that the following expression doesn't overflow. */ 281 expo = MPFR_GET_EXP (y) + expo - r * (i - 1); 282 MPFR_SET_EXP (y, expo); 283} 284 285int 286mpfr_atan (mpfr_ptr atan, mpfr_srcptr x, mpfr_rnd_t rnd_mode) 287{ 288 mpfr_t xp, arctgt, sk, tmp, tmp2; 289 mpz_t ukz; 290 mpz_t tabz[3*(MPFR_PREC_BITS+1)]; 291 mpfr_exp_t exptol; 292 mpfr_prec_t prec, realprec, est_lost, lost; 293 unsigned long twopoweri, log2p, red; 294 int comparison, inexact; 295 int i, n0, oldn0; 296 MPFR_GROUP_DECL (group); 297 MPFR_SAVE_EXPO_DECL (expo); 298 MPFR_ZIV_DECL (loop); 299 300 MPFR_LOG_FUNC 301 (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd_mode), 302 ("atan[%Pu]=%.*Rg inexact=%d", 303 mpfr_get_prec (atan), mpfr_log_prec, atan, inexact)); 304 305 /* Singular cases */ 306 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) 307 { 308 if (MPFR_IS_NAN (x)) 309 { 310 MPFR_SET_NAN (atan); 311 MPFR_RET_NAN; 312 } 313 else if (MPFR_IS_INF (x)) 314 { 315 MPFR_SAVE_EXPO_MARK (expo); 316 if (MPFR_IS_POS (x)) /* arctan(+inf) = Pi/2 */ 317 inexact = mpfr_const_pi (atan, rnd_mode); 318 else /* arctan(-inf) = -Pi/2 */ 319 { 320 inexact = -mpfr_const_pi (atan, 321 MPFR_INVERT_RND (rnd_mode)); 322 MPFR_CHANGE_SIGN (atan); 323 } 324 mpfr_div_2ui (atan, atan, 1, rnd_mode); /* exact (no exceptions) */ 325 MPFR_SAVE_EXPO_FREE (expo); 326 return mpfr_check_range (atan, inexact, rnd_mode); 327 } 328 else /* x is necessarily 0 */ 329 { 330 MPFR_ASSERTD (MPFR_IS_ZERO (x)); 331 MPFR_SET_ZERO (atan); 332 MPFR_SET_SAME_SIGN (atan, x); 333 MPFR_RET (0); 334 } 335 } 336 337 /* atan(x) = x - x^3/3 + x^5/5... 338 so the error is < 2^(3*EXP(x)-1) 339 so `EXP(x)-(3*EXP(x)-1)` = -2*EXP(x)+1 */ 340 MPFR_FAST_COMPUTE_IF_SMALL_INPUT (atan, x, -2 * MPFR_GET_EXP (x), 1, 0, 341 rnd_mode, {}); 342 343 /* Set x_p=|x| */ 344 MPFR_TMP_INIT_ABS (xp, x); 345 346 MPFR_SAVE_EXPO_MARK (expo); 347 348 /* Other simple case arctan(-+1)=-+pi/4 */ 349 comparison = mpfr_cmp_ui (xp, 1); 350 if (MPFR_UNLIKELY (comparison == 0)) 351 { 352 int neg = MPFR_IS_NEG (x); 353 inexact = mpfr_const_pi (atan, MPFR_IS_POS (x) ? rnd_mode 354 : MPFR_INVERT_RND (rnd_mode)); 355 if (neg) 356 { 357 inexact = -inexact; 358 MPFR_CHANGE_SIGN (atan); 359 } 360 mpfr_div_2ui (atan, atan, 2, rnd_mode); /* exact (no exceptions) */ 361 MPFR_SAVE_EXPO_FREE (expo); 362 return mpfr_check_range (atan, inexact, rnd_mode); 363 } 364 365 realprec = MPFR_PREC (atan) + MPFR_INT_CEIL_LOG2 (MPFR_PREC (atan)) + 4; 366 prec = realprec + GMP_NUMB_BITS; 367 368 /* Initialisation */ 369 mpz_init2 (ukz, prec); /* ukz will need 'prec' bits below */ 370 MPFR_GROUP_INIT_4 (group, prec, sk, tmp, tmp2, arctgt); 371 oldn0 = 0; 372 373 MPFR_ZIV_INIT (loop, prec); 374 for (;;) 375 { 376 /* First, if |x| < 1, we need to have more prec to be able to round (sup) 377 n0 = ceil(log(prec_requested + 2 + 1+ln(2.4)/ln(2))/log(2)) */ 378 mpfr_prec_t sup; 379 sup = MPFR_GET_EXP (xp) < 0 ? 2 - MPFR_GET_EXP (xp) : 1; /* sup >= 1 */ 380 381 n0 = MPFR_INT_CEIL_LOG2 ((realprec + sup) + 3); 382 /* since realprec >= 4, n0 >= ceil(log2(8)) >= 3, thus 3*n0 > 2 */ 383 prec = (realprec + sup) + 1 + MPFR_INT_CEIL_LOG2 (3*n0-2); 384 385 /* the number of lost bits due to argument reduction is 386 9 - 2 * EXP(sk), which we estimate by 9 + 2*ceil(log2(p)) 387 since we manage that sk < 1/p */ 388 if (MPFR_PREC (atan) > 100) 389 { 390 log2p = MPFR_INT_CEIL_LOG2(prec) / 2 - 3; 391 est_lost = 9 + 2 * log2p; 392 prec += est_lost; 393 } 394 else 395 log2p = est_lost = 0; /* don't reduce the argument */ 396 397 /* Initialisation */ 398 MPFR_GROUP_REPREC_4 (group, prec, sk, tmp, tmp2, arctgt); 399 MPFR_ASSERTD (n0 <= MPFR_PREC_BITS); 400 /* Note: the tabz[] entries are used to get a rational approximation 401 of atan(x) to precision 'prec', thus allocating them to 'prec' bits 402 should be good enough. */ 403 for (i = oldn0; i < 3 * (n0 + 1); i++) 404 mpz_init2 (tabz[i], prec); 405 oldn0 = 3 * (n0 + 1); 406 407 /* The mpfr_ui_div below mustn't underflow. This is guaranteed by 408 MPFR_SAVE_EXPO_MARK, but let's check that for maintainability. */ 409 MPFR_ASSERTD (__gmpfr_emax <= 1 - __gmpfr_emin); 410 411 if (comparison > 0) /* use atan(xp) = Pi/2 - atan(1/xp) */ 412 mpfr_ui_div (sk, 1, xp, MPFR_RNDN); 413 else 414 mpfr_set (sk, xp, MPFR_RNDN); 415 416 /* now 0 < sk <= 1 */ 417 418 /* Argument reduction: atan(x) = 2 atan((sqrt(1+x^2)-1)/x). 419 We want |sk| < k/sqrt(p) where p is the target precision. */ 420 lost = 0; 421 for (red = 0; MPFR_GET_EXP(sk) > - (mpfr_exp_t) log2p; red ++) 422 { 423 lost = 9 - 2 * MPFR_EXP(sk); 424 mpfr_sqr (tmp, sk, MPFR_RNDN); 425 mpfr_add_ui (tmp, tmp, 1, MPFR_RNDN); 426 mpfr_sqrt (tmp, tmp, MPFR_RNDN); 427 mpfr_sub_ui (tmp, tmp, 1, MPFR_RNDN); 428 if (red == 0 && comparison > 0) 429 /* use xp = 1/sk */ 430 mpfr_mul (sk, tmp, xp, MPFR_RNDN); 431 else 432 mpfr_div (sk, tmp, sk, MPFR_RNDN); 433 } 434 435 /* We started from x0 = 1/|x| if |x| > 1, and |x| otherwise, thus 436 we had x0 = min(|x|, 1/|x|) <= 1, and applied 'red' times the 437 argument reduction x -> (sqrt(1+x^2)-1)/x, which keeps 0 < x <= 1 */ 438 439 /* We first show that if the for-loop is executed at least once, then 440 sk < 1 after the loop. Indeed for 1/2 <= x <= 1, interval 441 arithmetic with precision 5 shows that (sqrt(1+x^2)-1)/x, 442 when evaluated with rounding to nearest, gives a value <= 0.875. 443 Now assume 2^(-k-1) <= x <= 2^(-k) for k >= 1. 444 Then o(x^2) <= 2^(-2k), o(1+x^2) <= 1+2^(-2k), 445 o(sqrt(1+x^2)) <= 1+2^(-2k-1), o(sqrt(1+x^2)-1) <= 2^(-2k-1), 446 and o((sqrt(1+x^2)-1)/x) <= 2^(-k) <= 1/2. 447 448 Now if sk=1 before the loop, then EXP(sk)=1 and since log2p >= 0, 449 the loop is performed at least once, thus the case sk=1 cannot 450 happen below. 451 */ 452 453 MPFR_ASSERTD(mpfr_cmp_ui (sk, 1) < 0); 454 455 /* Assignation */ 456 MPFR_SET_ZERO (arctgt); 457 twopoweri = 1 << 0; 458 MPFR_ASSERTD (n0 >= 4); 459 for (i = 0 ; i < n0; i++) 460 { 461 if (MPFR_UNLIKELY (MPFR_IS_ZERO (sk))) 462 break; 463 /* Calculation of trunc(tmp) --> mpz */ 464 mpfr_mul_2ui (tmp, sk, twopoweri, MPFR_RNDN); 465 mpfr_trunc (tmp, tmp); 466 if (!MPFR_IS_ZERO (tmp)) 467 { 468 /* tmp = ukz*2^exptol */ 469 exptol = mpfr_get_z_2exp (ukz, tmp); 470 /* since the s_k are decreasing (see algorithms.tex), 471 and s_0 = min(|x|, 1/|x|) < 1, we have sk < 1, 472 thus exptol < 0 */ 473 MPFR_ASSERTD (exptol < 0); 474 mpz_tdiv_q_2exp (ukz, ukz, (unsigned long int) (-exptol)); 475 /* since tmp is a non-zero integer, and tmp = ukzold*2^exptol, 476 we now have ukz = tmp, thus ukz is non-zero */ 477 /* Calculation of arctan(Ak) */ 478 mpfr_set_z (tmp, ukz, MPFR_RNDN); 479 mpfr_div_2ui (tmp, tmp, twopoweri, MPFR_RNDN); 480 mpfr_atan_aux (tmp2, ukz, twopoweri, n0 - i, tabz); 481 mpfr_mul (tmp2, tmp2, tmp, MPFR_RNDN); 482 /* Addition */ 483 mpfr_add (arctgt, arctgt, tmp2, MPFR_RNDN); 484 /* Next iteration */ 485 mpfr_sub (tmp2, sk, tmp, MPFR_RNDN); 486 mpfr_mul (sk, sk, tmp, MPFR_RNDN); 487 mpfr_add_ui (sk, sk, 1, MPFR_RNDN); 488 mpfr_div (sk, tmp2, sk, MPFR_RNDN); 489 } 490 twopoweri <<= 1; 491 } 492 /* Add last step (Arctan(sk) ~= sk */ 493 mpfr_add (arctgt, arctgt, sk, MPFR_RNDN); 494 495 /* argument reduction */ 496 mpfr_mul_2ui (arctgt, arctgt, red, MPFR_RNDN); 497 498 if (comparison > 0) 499 { /* atan(x) = Pi/2-atan(1/x) for x > 0 */ 500 mpfr_const_pi (tmp, MPFR_RNDN); 501 mpfr_div_2ui (tmp, tmp, 1, MPFR_RNDN); 502 mpfr_sub (arctgt, tmp, arctgt, MPFR_RNDN); 503 } 504 MPFR_SET_POS (arctgt); 505 506 if (MPFR_LIKELY (MPFR_CAN_ROUND (arctgt, realprec + est_lost - lost, 507 MPFR_PREC (atan), rnd_mode))) 508 break; 509 MPFR_ZIV_NEXT (loop, realprec); 510 } 511 MPFR_ZIV_FREE (loop); 512 513 inexact = mpfr_set4 (atan, arctgt, rnd_mode, MPFR_SIGN (x)); 514 515 for (i = 0 ; i < oldn0 ; i++) 516 mpz_clear (tabz[i]); 517 mpz_clear (ukz); 518 MPFR_GROUP_CLEAR (group); 519 520 MPFR_SAVE_EXPO_FREE (expo); 521 return mpfr_check_range (atan, inexact, rnd_mode); 522} 523