1/* mpfr_asin -- arc-sinus of a floating-point number
2
3Copyright 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
4Contributed by the Arenaire and Cacao projects, INRIA.
5
6This file is part of the GNU MPFR Library, and was contributed by Mathieu Dutour.
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#include "mpfr-impl.h"
24
25int
26mpfr_asin (mpfr_ptr asin, mpfr_srcptr x, mpfr_rnd_t rnd_mode)
27{
28  mpfr_t xp;
29  int compared, inexact;
30  mpfr_prec_t prec;
31  mpfr_exp_t xp_exp;
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                 ("asin[%#R]=%R inexact=%d", asin, asin, inexact));
37
38  /* Special 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 (asin);
44          MPFR_RET_NAN;
45        }
46      else /* x = 0 */
47        {
48          MPFR_ASSERTD (MPFR_IS_ZERO (x));
49          MPFR_SET_ZERO (asin);
50          MPFR_SET_SAME_SIGN (asin, x);
51          MPFR_RET (0); /* exact result */
52        }
53    }
54
55  /* asin(x) = x + x^3/6 + ... so the error is < 2^(3*EXP(x)-2) */
56  MPFR_FAST_COMPUTE_IF_SMALL_INPUT (asin, x, -2 * MPFR_GET_EXP (x), 2, 1,
57                                    rnd_mode, {});
58
59  /* Set x_p=|x| (x is a normal number) */
60  mpfr_init2 (xp, MPFR_PREC (x));
61  inexact = mpfr_abs (xp, x, MPFR_RNDN);
62  MPFR_ASSERTD (inexact == 0);
63
64  compared = mpfr_cmp_ui (xp, 1);
65
66  MPFR_SAVE_EXPO_MARK (expo);
67
68  if (MPFR_UNLIKELY (compared >= 0))
69    {
70      mpfr_clear (xp);
71      if (compared > 0)                  /* asin(x) = NaN for |x| > 1 */
72        {
73          MPFR_SAVE_EXPO_FREE (expo);
74          MPFR_SET_NAN (asin);
75          MPFR_RET_NAN;
76        }
77      else                              /* x = 1 or x = -1 */
78        {
79          if (MPFR_IS_POS (x)) /* asin(+1) = Pi/2 */
80            inexact = mpfr_const_pi (asin, rnd_mode);
81          else /* asin(-1) = -Pi/2 */
82            {
83              inexact = -mpfr_const_pi (asin, MPFR_INVERT_RND(rnd_mode));
84              MPFR_CHANGE_SIGN (asin);
85            }
86          mpfr_div_2ui (asin, asin, 1, rnd_mode);
87        }
88    }
89  else
90    {
91  /* Compute exponent of 1 - ABS(x) */
92  mpfr_ui_sub (xp, 1, xp, MPFR_RNDD);
93  MPFR_ASSERTD (MPFR_GET_EXP (xp) <= 0);
94  MPFR_ASSERTD (MPFR_GET_EXP (x) <= 0);
95  xp_exp = 2 - MPFR_GET_EXP (xp);
96
97  /* Set up initial prec */
98  prec = MPFR_PREC (asin) + 10 + xp_exp;
99
100  /* use asin(x) = atan(x/sqrt(1-x^2)) */
101  MPFR_ZIV_INIT (loop, prec);
102  for (;;)
103    {
104      mpfr_set_prec (xp, prec);
105      mpfr_sqr (xp, x, MPFR_RNDN);
106      mpfr_ui_sub (xp, 1, xp, MPFR_RNDN);
107      mpfr_sqrt (xp, xp, MPFR_RNDN);
108      mpfr_div (xp, x, xp, MPFR_RNDN);
109      mpfr_atan (xp, xp, MPFR_RNDN);
110      if (MPFR_LIKELY (MPFR_CAN_ROUND (xp, prec - xp_exp,
111                                       MPFR_PREC (asin), rnd_mode)))
112        break;
113      MPFR_ZIV_NEXT (loop, prec);
114    }
115  MPFR_ZIV_FREE (loop);
116  inexact = mpfr_set (asin, xp, rnd_mode);
117
118  mpfr_clear (xp);
119    }
120
121  MPFR_SAVE_EXPO_FREE (expo);
122  return mpfr_check_range (asin, inexact, rnd_mode);
123}
124