1 /* mpfr_atan -- arc-tangent of a floating-point number
2
3 Copyright 2001-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 #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 */
35 static 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
58 static void
59 set_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 */
91 static void
92 mpfr_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
285 int
286 mpfr_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[%Pd]=%.*Rg rnd=%d", mpfr_get_prec (x), mpfr_log_prec, x, rnd_mode),
302 ("atan[%Pd]=%.*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 }