OLD | NEW |
| (Empty) |
1 /* mpfr_acos -- arc-cosinus of a floating-point number | |
2 | |
3 Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009 Free Software Fou
ndation, Inc. | |
4 Contributed by the Arenaire and Cacao projects, INRIA. | |
5 | |
6 This file is part of the GNU MPFR Library, and was contributed by Mathieu Dutour
. | |
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 2.1 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.LIB. If not, write to | |
20 the Free Software Foundation, Inc., 51 Franklin St, Fifth Floor, Boston, | |
21 MA 02110-1301, USA. */ | |
22 | |
23 #include "mpfr-impl.h" | |
24 | |
25 int | |
26 mpfr_acos (mpfr_ptr acos, mpfr_srcptr x, mp_rnd_t rnd_mode) | |
27 { | |
28 mpfr_t xp, arcc, tmp; | |
29 mp_exp_t supplement; | |
30 mp_prec_t prec; | |
31 int sign, compared, inexact; | |
32 MPFR_SAVE_EXPO_DECL (expo); | |
33 MPFR_ZIV_DECL (loop); | |
34 | |
35 MPFR_LOG_FUNC (("x[%#R]=%R rnd=%d", x, x, rnd_mode), | |
36 ("acos[%#R]=%R inexact=%d", acos, acos, inexact)); | |
37 | |
38 /* Singular cases */ | |
39 if (MPFR_UNLIKELY (MPFR_IS_SINGULAR (x))) | |
40 { | |
41 if (MPFR_IS_NAN (x) || MPFR_IS_INF (x)) | |
42 { | |
43 MPFR_SET_NAN (acos); | |
44 MPFR_RET_NAN; | |
45 } | |
46 else /* necessarily x=0 */ | |
47 { | |
48 MPFR_ASSERTD(MPFR_IS_ZERO(x)); | |
49 /* acos(0)=Pi/2 */ | |
50 inexact = mpfr_const_pi (acos, rnd_mode); | |
51 mpfr_div_2ui (acos, acos, 1, rnd_mode); /* exact */ | |
52 MPFR_RET (inexact); | |
53 } | |
54 } | |
55 | |
56 /* Set x_p=|x| */ | |
57 sign = MPFR_SIGN (x); | |
58 mpfr_init2 (xp, MPFR_PREC (x)); | |
59 mpfr_abs (xp, x, GMP_RNDN); /* Exact */ | |
60 | |
61 compared = mpfr_cmp_ui (xp, 1); | |
62 | |
63 if (MPFR_UNLIKELY (compared >= 0)) | |
64 { | |
65 mpfr_clear (xp); | |
66 if (compared > 0) /* acos(x) = NaN for x > 1 */ | |
67 { | |
68 MPFR_SET_NAN(acos); | |
69 MPFR_RET_NAN; | |
70 } | |
71 else | |
72 { | |
73 if (MPFR_IS_POS_SIGN (sign)) /* acos(+1) = 0 */ | |
74 return mpfr_set_ui (acos, 0, rnd_mode); | |
75 else /* acos(-1) = Pi */ | |
76 return mpfr_const_pi (acos, rnd_mode); | |
77 } | |
78 } | |
79 | |
80 MPFR_SAVE_EXPO_MARK (expo); | |
81 | |
82 /* Compute the supplement */ | |
83 mpfr_ui_sub (xp, 1, xp, GMP_RNDD); | |
84 if (MPFR_IS_POS_SIGN (sign)) | |
85 supplement = 2 - 2 * MPFR_GET_EXP (xp); | |
86 else | |
87 supplement = 2 - MPFR_GET_EXP (xp); | |
88 mpfr_clear (xp); | |
89 | |
90 prec = MPFR_PREC (acos) + 10 + supplement; | |
91 | |
92 /* VL: The following change concerning prec comes from r3145 | |
93 "Optimize mpfr_acos by choosing a better initial precision." | |
94 but it doesn't seem to be correct and leads to problems (assertion | |
95 failure or very important inefficiency) with tiny arguments. | |
96 Therefore, I've disabled it. */ | |
97 /* If x ~ 2^-N, acos(x) ~ PI/2 - x - x^3/6 | |
98 If Prec < 2*N, we can't round since x^3/6 won't be counted. */ | |
99 #if 0 | |
100 if (MPFR_PREC (acos) >= MPFR_PREC (x) && MPFR_GET_EXP (x) < 0) | |
101 { | |
102 mpfr_uexp_t pmin = (mpfr_uexp_t) (-2 * MPFR_GET_EXP (x)) + 5; | |
103 MPFR_ASSERTN (pmin <= MPFR_PREC_MAX); | |
104 if (prec < pmin) | |
105 prec = pmin; | |
106 } | |
107 #endif | |
108 | |
109 mpfr_init2 (tmp, prec); | |
110 mpfr_init2 (arcc, prec); | |
111 | |
112 MPFR_ZIV_INIT (loop, prec); | |
113 for (;;) | |
114 { | |
115 /* acos(x) = Pi/2 - asin(x) = Pi/2 - atan(x/sqrt(1-x^2)) */ | |
116 mpfr_sqr (tmp, x, GMP_RNDN); | |
117 mpfr_ui_sub (tmp, 1, tmp, GMP_RNDN); | |
118 mpfr_sqrt (tmp, tmp, GMP_RNDN); | |
119 mpfr_div (tmp, x, tmp, GMP_RNDN); | |
120 mpfr_atan (arcc, tmp, GMP_RNDN); | |
121 mpfr_const_pi (tmp, GMP_RNDN); | |
122 mpfr_div_2ui (tmp, tmp, 1, GMP_RNDN); | |
123 mpfr_sub (arcc, tmp, arcc, GMP_RNDN); | |
124 | |
125 if (MPFR_LIKELY (MPFR_CAN_ROUND (arcc, prec-supplement, | |
126 MPFR_PREC (acos), rnd_mode))) | |
127 break; | |
128 MPFR_ZIV_NEXT (loop, prec); | |
129 mpfr_set_prec (tmp, prec); | |
130 mpfr_set_prec (arcc, prec); | |
131 } | |
132 MPFR_ZIV_FREE (loop); | |
133 | |
134 inexact = mpfr_set (acos, arcc, rnd_mode); | |
135 mpfr_clear (tmp); | |
136 mpfr_clear (arcc); | |
137 | |
138 MPFR_SAVE_EXPO_FREE (expo); | |
139 return mpfr_check_range (acos, inexact, rnd_mode); | |
140 } | |
OLD | NEW |