1 /* mpfr_atan -- arc-tangent of a floating-point number
3 Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011, 2012, 2013 Free Software Foundation, Inc.
4 Contributed by the AriC and Caramel projects, INRIA.
6 This file is part of the GNU MPFR Library.
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.
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.
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. */
23 #define MPFR_NEED_LONGLONG_H
24 #include "mpfr-impl.h"
26 /* If x = p/2^r, put in y an approximation of atan(x)/x using 2^m terms
27 for the series expansion, with an error of at most 1 ulp.
30 If X=x^2, we want 1 - X/3 + X^2/5 - ... + (-1)^k*X^k/(2k+1) + ...
34 When we sum terms up to x^k/(2k+1), the denominator Q[0] is
35 3*5*7*...*(2k+1) ~ (2k/e)^k.
38 mpfr_atan_aux (mpfr_ptr y
, mpz_ptr p
, long r
, int m
, mpz_t
*tab
)
41 unsigned long n
, i
, k
, j
, l
;
42 mpfr_exp_t diff
, expo
;
44 mpfr_prec_t mult
, *accu
, *log2_nb_terms
;
45 mpfr_prec_t precy
= MPFR_PREC(y
);
47 MPFR_ASSERTD(mpz_cmp_ui (p
, 0) != 0);
49 accu
= (mpfr_prec_t
*) (*__gmp_allocate_func
) ((2 * m
+ 2) * sizeof (mpfr_prec_t
));
50 log2_nb_terms
= accu
+ m
+ 1;
54 ptoj
= S
+ 1*(m
+1); /* p^2^j Precomputed table */
55 Q
= S
+ 2*(m
+1); /* Product of Odd integer table */
57 /* From p to p^2, and r to 2r */
59 MPFR_ASSERTD (2 * r
> r
);
64 mpz_tdiv_q_2exp (p
, p
, n
); /* exact */
67 /* since |p/2^r| < 1, and p is a non-zero integer, necessarily r > 0 */
69 MPFR_ASSERTD (mpz_sgn (p
) > 0);
72 /* check if p=1 (special case) */
75 We compute by binary splitting, with X = x^2 = p/2^r:
76 P(a,b) = p if a+1=b, P(a,c)*P(c,b) otherwise
77 Q(a,b) = (2a+1)*2^r if a+1=b [except Q(0,1)=1], Q(a,c)*Q(c,b) otherwise
78 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
79 Then atan(x)/x ~ S(0,i)/Q(0,i) for i so that (p/2^r)^i/i is small enough.
80 The factor 2^(r*(b-a)) in Q(a,b) is implicit, thus we have to take it
81 into account when we compute with Q.
83 accu
[0] = 0; /* accu[k] = Mult[0] + ... + Mult[k], where Mult[j] is the
84 number of bits of the corresponding term S[j]/Q[j] */
85 if (mpz_cmp_ui (p
, 1) != 0)
87 /* p <> 1: precompute ptoj table */
89 for (im
= 1 ; im
<= m
; im
++)
90 mpz_mul (ptoj
[im
], ptoj
[im
- 1], ptoj
[im
- 1]);
93 /* the ith term being X^i/(2i+1) with X=p/2^r, we can stop when
94 p^i/2^(r*i) < 2^(-precy), i.e. r*i > precy + log2(p^i) */
95 for (i
= k
= done
= 0; (i
< n
) && (done
== 0); i
+= 2, k
++)
97 /* initialize both S[k],Q[k] and S[k+1],Q[k+1] */
98 mpz_set_ui (Q
[k
+1], 2 * i
+ 3); /* Q(i+1,i+2) */
99 mpz_mul_ui (S
[k
+1], p
, 2 * i
+ 1); /* S(i+1,i+2) */
100 mpz_mul_2exp (S
[k
], Q
[k
+1], r
);
101 mpz_sub (S
[k
], S
[k
], S
[k
+1]); /* S(i,i+2) */
102 mpz_mul_ui (Q
[k
], Q
[k
+1], 2 * i
+ 1); /* Q(i,i+2) */
103 log2_nb_terms
[k
] = 1; /* S[k]/Q[k] corresponds to 2 terms */
104 for (j
= (i
+ 2) >> 1, l
= 1; (j
& 1) == 0; l
++, j
>>= 1, k
--)
106 /* invariant: S[k-1]/Q[k-1] and S[k]/Q[k] correspond
107 to 2^l terms each. We combine them into S[k-1]/Q[k-1] */
108 MPFR_ASSERTD (k
> 0);
109 mpz_mul (S
[k
], S
[k
], Q
[k
-1]);
110 mpz_mul (S
[k
], S
[k
], ptoj
[l
]);
111 mpz_mul (S
[k
-1], S
[k
-1], Q
[k
]);
112 mpz_mul_2exp (S
[k
-1], S
[k
-1], r
<< l
);
113 mpz_add (S
[k
-1], S
[k
-1], S
[k
]);
114 mpz_mul (Q
[k
-1], Q
[k
-1], Q
[k
]);
115 log2_nb_terms
[k
-1] = l
+ 1;
116 /* now S[k-1]/Q[k-1] corresponds to 2^(l+1) terms */
117 MPFR_MPZ_SIZEINBASE2(mult
, ptoj
[l
+1]);
118 /* FIXME: precompute bits(ptoj[l+1]) outside the loop? */
119 mult
= (r
<< (l
+ 1)) - mult
- 1;
120 accu
[k
-1] = (k
== 1) ? mult
: accu
[k
-2] + mult
;
121 if (accu
[k
-1] > precy
)
126 else /* special case p=1: the ith term being X^i/(2i+1) with X=1/2^r,
127 we can stop when r*i > precy i.e. i > precy/r */
130 for (i
= k
= 0; (i
< n
) && (i
<= precy
/ r
); i
+= 2, k
++)
132 mpz_set_ui (Q
[k
+ 1], 2 * i
+ 3);
133 mpz_mul_2exp (S
[k
], Q
[k
+1], r
);
134 mpz_sub_ui (S
[k
], S
[k
], 1 + 2 * i
);
135 mpz_mul_ui (Q
[k
], Q
[k
+ 1], 1 + 2 * i
);
136 log2_nb_terms
[k
] = 1; /* S[k]/Q[k] corresponds to 2 terms */
137 for (j
= (i
+ 2) >> 1, l
= 1; (j
& 1) == 0; l
++, j
>>= 1, k
--)
139 MPFR_ASSERTD (k
> 0);
140 mpz_mul (S
[k
], S
[k
], Q
[k
-1]);
141 mpz_mul (S
[k
-1], S
[k
-1], Q
[k
]);
142 mpz_mul_2exp (S
[k
-1], S
[k
-1], r
<< l
);
143 mpz_add (S
[k
-1], S
[k
-1], S
[k
]);
144 mpz_mul (Q
[k
-1], Q
[k
-1], Q
[k
]);
145 log2_nb_terms
[k
-1] = l
+ 1;
150 /* we need to combine S[0]/Q[0]...S[k-1]/Q[k-1] */
151 l
= 0; /* number of terms accumulated in S[k]/Q[k] */
155 /* combine S[k-1]/Q[k-1] and S[k]/Q[k] */
156 j
= log2_nb_terms
[k
-1];
157 mpz_mul (S
[k
], S
[k
], Q
[k
-1]);
158 if (mpz_cmp_ui (p
, 1) != 0)
159 mpz_mul (S
[k
], S
[k
], ptoj
[j
]);
160 mpz_mul (S
[k
-1], S
[k
-1], Q
[k
]);
161 l
+= 1 << log2_nb_terms
[k
];
162 mpz_mul_2exp (S
[k
-1], S
[k
-1], r
* l
);
163 mpz_add (S
[k
-1], S
[k
-1], S
[k
]);
164 mpz_mul (Q
[k
-1], Q
[k
-1], Q
[k
]);
166 (*__gmp_free_func
) (accu
, (2 * m
+ 2) * sizeof (mpfr_prec_t
));
168 MPFR_MPZ_SIZEINBASE2 (diff
, S
[0]);
172 mpz_tdiv_q_2exp (S
[0], S
[0], diff
);
174 mpz_mul_2exp (S
[0], S
[0], -diff
);
176 MPFR_MPZ_SIZEINBASE2 (diff
, Q
[0]);
180 mpz_tdiv_q_2exp (Q
[0], Q
[0], diff
);
182 mpz_mul_2exp (Q
[0], Q
[0], -diff
);
184 mpz_tdiv_q (S
[0], S
[0], Q
[0]);
185 mpfr_set_z (y
, S
[0], MPFR_RNDD
);
186 MPFR_SET_EXP (y
, MPFR_EXP(y
) + expo
- r
* (i
- 1));
190 mpfr_atan (mpfr_ptr atan
, mpfr_srcptr x
, mpfr_rnd_t rnd_mode
)
192 mpfr_t xp
, arctgt
, sk
, tmp
, tmp2
;
196 mpfr_prec_t prec
, realprec
, est_lost
, lost
;
197 unsigned long twopoweri
, log2p
, red
;
198 int comparaison
, inexact
;
200 MPFR_GROUP_DECL (group
);
201 MPFR_SAVE_EXPO_DECL (expo
);
202 MPFR_ZIV_DECL (loop
);
205 (("x[%Pu]=%.*Rg rnd=%d", mpfr_get_prec (x
), mpfr_log_prec
, x
, rnd_mode
),
206 ("atan[%Pu]=%.*Rg inexact=%d",
207 mpfr_get_prec (atan
), mpfr_log_prec
, atan
, inexact
));
210 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x
)))
217 else if (MPFR_IS_INF (x
))
219 MPFR_SAVE_EXPO_MARK (expo
);
220 if (MPFR_IS_POS (x
)) /* arctan(+inf) = Pi/2 */
221 inexact
= mpfr_const_pi (atan
, rnd_mode
);
222 else /* arctan(-inf) = -Pi/2 */
224 inexact
= -mpfr_const_pi (atan
,
225 MPFR_INVERT_RND (rnd_mode
));
226 MPFR_CHANGE_SIGN (atan
);
228 mpfr_div_2ui (atan
, atan
, 1, rnd_mode
); /* exact (no exceptions) */
229 MPFR_SAVE_EXPO_FREE (expo
);
230 return mpfr_check_range (atan
, inexact
, rnd_mode
);
232 else /* x is necessarily 0 */
234 MPFR_ASSERTD (MPFR_IS_ZERO (x
));
235 MPFR_SET_ZERO (atan
);
236 MPFR_SET_SAME_SIGN (atan
, x
);
241 /* atan(x) = x - x^3/3 + x^5/5...
242 so the error is < 2^(3*EXP(x)-1)
243 so `EXP(x)-(3*EXP(x)-1)` = -2*EXP(x)+1 */
244 MPFR_FAST_COMPUTE_IF_SMALL_INPUT (atan
, x
, -2 * MPFR_GET_EXP (x
), 1, 0,
248 MPFR_TMP_INIT_ABS (xp
, x
);
250 MPFR_SAVE_EXPO_MARK (expo
);
252 /* Other simple case arctan(-+1)=-+pi/4 */
253 comparaison
= mpfr_cmp_ui (xp
, 1);
254 if (MPFR_UNLIKELY (comparaison
== 0))
256 int neg
= MPFR_IS_NEG (x
);
257 inexact
= mpfr_const_pi (atan
, MPFR_IS_POS (x
) ? rnd_mode
258 : MPFR_INVERT_RND (rnd_mode
));
262 MPFR_CHANGE_SIGN (atan
);
264 mpfr_div_2ui (atan
, atan
, 2, rnd_mode
); /* exact (no exceptions) */
265 MPFR_SAVE_EXPO_FREE (expo
);
266 return mpfr_check_range (atan
, inexact
, rnd_mode
);
269 realprec
= MPFR_PREC (atan
) + MPFR_INT_CEIL_LOG2 (MPFR_PREC (atan
)) + 4;
270 prec
= realprec
+ GMP_NUMB_BITS
;
274 MPFR_GROUP_INIT_4 (group
, prec
, sk
, tmp
, tmp2
, arctgt
);
278 MPFR_ZIV_INIT (loop
, prec
);
281 /* First, if |x| < 1, we need to have more prec to be able to round (sup)
282 n0 = ceil(log(prec_requested + 2 + 1+ln(2.4)/ln(2))/log(2)) */
284 sup
= MPFR_GET_EXP (xp
) < 0 ? 2 - MPFR_GET_EXP (xp
) : 1; /* sup >= 1 */
286 n0
= MPFR_INT_CEIL_LOG2 ((realprec
+ sup
) + 3);
287 /* since realprec >= 4, n0 >= ceil(log2(8)) >= 3, thus 3*n0 > 2 */
288 prec
= (realprec
+ sup
) + 1 + MPFR_INT_CEIL_LOG2 (3*n0
-2);
290 /* the number of lost bits due to argument reduction is
291 9 - 2 * EXP(sk), which we estimate by 9 + 2*ceil(log2(p))
292 since we manage that sk < 1/p */
293 if (MPFR_PREC (atan
) > 100)
295 log2p
= MPFR_INT_CEIL_LOG2(prec
) / 2 - 3;
296 est_lost
= 9 + 2 * log2p
;
300 log2p
= est_lost
= 0; /* don't reduce the argument */
303 MPFR_GROUP_REPREC_4 (group
, prec
, sk
, tmp
, tmp2
, arctgt
);
304 if (MPFR_LIKELY (oldn0
== 0))
306 oldn0
= 3 * (n0
+ 1);
307 tabz
= (mpz_t
*) (*__gmp_allocate_func
) (oldn0
* sizeof (mpz_t
));
308 for (i
= 0; i
< oldn0
; i
++)
311 else if (MPFR_UNLIKELY (oldn0
< 3 * (n0
+ 1)))
313 tabz
= (mpz_t
*) (*__gmp_reallocate_func
)
314 (tabz
, oldn0
* sizeof (mpz_t
), 3 * (n0
+ 1)*sizeof (mpz_t
));
315 for (i
= oldn0
; i
< 3 * (n0
+ 1); i
++)
317 oldn0
= 3 * (n0
+ 1);
320 /* The mpfr_ui_div below mustn't underflow. This is guaranteed by
321 MPFR_SAVE_EXPO_MARK, but let's check that for maintainability. */
322 MPFR_ASSERTD (__gmpfr_emax
<= 1 - __gmpfr_emin
);
324 if (comparaison
> 0) /* use atan(xp) = Pi/2 - atan(1/xp) */
325 mpfr_ui_div (sk
, 1, xp
, MPFR_RNDN
);
327 mpfr_set (sk
, xp
, MPFR_RNDN
);
329 /* now 0 < sk <= 1 */
331 /* Argument reduction: atan(x) = 2 atan((sqrt(1+x^2)-1)/x).
332 We want |sk| < k/sqrt(p) where p is the target precision. */
334 for (red
= 0; MPFR_GET_EXP(sk
) > - (mpfr_exp_t
) log2p
; red
++)
336 lost
= 9 - 2 * MPFR_EXP(sk
);
337 mpfr_mul (tmp
, sk
, sk
, MPFR_RNDN
);
338 mpfr_add_ui (tmp
, tmp
, 1, MPFR_RNDN
);
339 mpfr_sqrt (tmp
, tmp
, MPFR_RNDN
);
340 mpfr_sub_ui (tmp
, tmp
, 1, MPFR_RNDN
);
341 if (red
== 0 && comparaison
> 0)
343 mpfr_mul (sk
, tmp
, xp
, MPFR_RNDN
);
345 mpfr_div (sk
, tmp
, sk
, MPFR_RNDN
);
348 /* we started from x0 = 1/|x| if |x| > 1, and |x| otherwise, thus
349 we had x0 = min(|x|, 1/|x|) <= 1, and applied 'red' times the
350 argument reduction x -> (sqrt(1+x^2)-1)/x, which keeps 0 < x < 1,
351 thus 0 < sk <= 1, and sk=1 can occur only if red=0 */
353 /* If sk=1, then if |x| < 1, we have 1 - 2^(-prec-1) <= |x| < 1,
354 or if |x| > 1, we have 1 - 2^(-prec-1) <= 1/|x| < 1, thus in all
355 cases ||x| - 1| <= 2^(-prec), from which it follows
356 |atan|x| - Pi/4| <= 2^(-prec), given the Taylor expansion
357 atan(1+x) = Pi/4 + x/2 - x^2/4 + ...
358 Since Pi/4 = 0.785..., the error is at most one ulp.
360 if (MPFR_UNLIKELY(mpfr_cmp_ui (sk
, 1) == 0))
362 mpfr_const_pi (arctgt
, MPFR_RNDN
); /* 1/2 ulp extra error */
363 mpfr_div_2ui (arctgt
, arctgt
, 2, MPFR_RNDN
); /* exact */
369 MPFR_SET_ZERO (arctgt
);
371 MPFR_ASSERTD (n0
>= 4);
372 for (i
= 0 ; i
< n0
; i
++)
374 if (MPFR_UNLIKELY (MPFR_IS_ZERO (sk
)))
376 /* Calculation of trunc(tmp) --> mpz */
377 mpfr_mul_2ui (tmp
, sk
, twopoweri
, MPFR_RNDN
);
378 mpfr_trunc (tmp
, tmp
);
379 if (!MPFR_IS_ZERO (tmp
))
381 /* tmp = ukz*2^exptol */
382 exptol
= mpfr_get_z_2exp (ukz
, tmp
);
383 /* since the s_k are decreasing (see algorithms.tex),
384 and s_0 = min(|x|, 1/|x|) < 1, we have sk < 1,
386 MPFR_ASSERTD (exptol
< 0);
387 mpz_tdiv_q_2exp (ukz
, ukz
, (unsigned long int) (-exptol
));
388 /* since tmp is a non-zero integer, and tmp = ukzold*2^exptol,
389 we now have ukz = tmp, thus ukz is non-zero */
390 /* Calculation of arctan(Ak) */
391 mpfr_set_z (tmp
, ukz
, MPFR_RNDN
);
392 mpfr_div_2ui (tmp
, tmp
, twopoweri
, MPFR_RNDN
);
393 mpfr_atan_aux (tmp2
, ukz
, twopoweri
, n0
- i
, tabz
);
394 mpfr_mul (tmp2
, tmp2
, tmp
, MPFR_RNDN
);
396 mpfr_add (arctgt
, arctgt
, tmp2
, MPFR_RNDN
);
398 mpfr_sub (tmp2
, sk
, tmp
, MPFR_RNDN
);
399 mpfr_mul (sk
, sk
, tmp
, MPFR_RNDN
);
400 mpfr_add_ui (sk
, sk
, 1, MPFR_RNDN
);
401 mpfr_div (sk
, tmp2
, sk
, MPFR_RNDN
);
405 /* Add last step (Arctan(sk) ~= sk */
406 mpfr_add (arctgt
, arctgt
, sk
, MPFR_RNDN
);
408 /* argument reduction */
409 mpfr_mul_2exp (arctgt
, arctgt
, red
, MPFR_RNDN
);
412 { /* atan(x) = Pi/2-atan(1/x) for x > 0 */
413 mpfr_const_pi (tmp
, MPFR_RNDN
);
414 mpfr_div_2ui (tmp
, tmp
, 1, MPFR_RNDN
);
415 mpfr_sub (arctgt
, tmp
, arctgt
, MPFR_RNDN
);
417 MPFR_SET_POS (arctgt
);
420 if (MPFR_LIKELY (MPFR_CAN_ROUND (arctgt
, realprec
+ est_lost
- lost
,
421 MPFR_PREC (atan
), rnd_mode
)))
423 MPFR_ZIV_NEXT (loop
, realprec
);
425 MPFR_ZIV_FREE (loop
);
427 inexact
= mpfr_set4 (atan
, arctgt
, rnd_mode
, MPFR_SIGN (x
));
429 for (i
= 0 ; i
< oldn0
; i
++)
432 (*__gmp_free_func
) (tabz
, oldn0
* sizeof (mpz_t
));
433 MPFR_GROUP_CLEAR (group
);
435 MPFR_SAVE_EXPO_FREE (expo
);
436 return mpfr_check_range (atan
, inexact
, rnd_mode
);