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