1/*
2 * ====================================================
3 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
4 *
5 * Developed at SunPro, a Sun Microsystems, Inc. business.
6 * Permission to use, copy, modify, and distribute this
7 * software is freely granted, provided that this notice
8 * is preserved.
9 * ====================================================
10 */
11
12/*
13  Long double expansions are
14  Copyright (C) 2001 Stephen L. Moshier <moshier@na-net.ornl.gov>
15  and are incorporated herein by permission of the author.  The author
16  reserves the right to distribute this material elsewhere under different
17  copying permissions.  These modifications are distributed here under
18  the following terms:
19
20    This library is free software; you can redistribute it and/or
21    modify it under the terms of the GNU Lesser General Public
22    License as published by the Free Software Foundation; either
23    version 2.1 of the License, or (at your option) any later version.
24
25    This library is distributed in the hope that it will be useful,
26    but WITHOUT ANY WARRANTY; without even the implied warranty of
27    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
28    Lesser General Public License for more details.
29
30    You should have received a copy of the GNU Lesser General Public
31    License along with this library; if not, write to the Free Software
32    Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307  USA */
33
34/* __ieee754_asin(x)
35 * Method :
36 *	Since  asin(x) = x + x^3/6 + x^5*3/40 + x^7*15/336 + ...
37 *	we approximate asin(x) on [0,0.5] by
38 *		asin(x) = x + x*x^2*R(x^2)
39 *
40 *	For x in [0.5,1]
41 *		asin(x) = pi/2-2*asin(sqrt((1-x)/2))
42 *	Let y = (1-x), z = y/2, s := sqrt(z), and pio2_hi+pio2_lo=pi/2;
43 *	then for x>0.98
44 *		asin(x) = pi/2 - 2*(s+s*z*R(z))
45 *			= pio2_hi - (2*(s+s*z*R(z)) - pio2_lo)
46 *	For x<=0.98, let pio4_hi = pio2_hi/2, then
47 *		f = hi part of s;
48 *		c = sqrt(z) - f = (z-f*f)/(s+f) 	...f+c=sqrt(z)
49 *	and
50 *		asin(x) = pi/2 - 2*(s+s*z*R(z))
51 *			= pio4_hi+(pio4-2s)-(2s*z*R(z)-pio2_lo)
52 *			= pio4_hi+(pio4-2f)-(2s*z*R(z)-(pio2_lo+2c))
53 *
54 * Special cases:
55 *	if x is NaN, return x itself;
56 *	if |x|>1, return NaN with invalid signal.
57 *
58 */
59
60
61#include "math.h"
62#include "math_private.h"
63
64#ifdef __STDC__
65static const long double
66#else
67static long double
68#endif
69  one = 1.0L,
70  huge = 1.0e+4932L,
71 pio2_hi = 1.5707963267948966192021943710788178805159986950457096099853515625L,
72  pio2_lo = 2.9127320560933561582586004641843300502121E-20L,
73  pio4_hi = 7.8539816339744830960109718553940894025800E-1L,
74
75	/* coefficient for R(x^2) */
76
77  /* asin(x) = x + x^3 pS(x^2) / qS(x^2)
78     0 <= x <= 0.5
79     peak relative error 1.9e-21  */
80  pS0 =  -1.008714657938491626019651170502036851607E1L,
81  pS1 =   2.331460313214179572063441834101394865259E1L,
82  pS2 =  -1.863169762159016144159202387315381830227E1L,
83  pS3 =   5.930399351579141771077475766877674661747E0L,
84  pS4 =  -6.121291917696920296944056882932695185001E-1L,
85  pS5 =   3.776934006243367487161248678019350338383E-3L,
86
87  qS0 =  -6.052287947630949712886794360635592886517E1L,
88  qS1 =   1.671229145571899593737596543114258558503E2L,
89  qS2 =  -1.707840117062586426144397688315411324388E2L,
90  qS3 =   7.870295154902110425886636075950077640623E1L,
91  qS4 =  -1.568433562487314651121702982333303458814E1L;
92    /* 1.000000000000000000000000000000000000000E0 */
93
94#ifdef __STDC__
95long double
96__ieee754_asinl (long double x)
97#else
98double
99__ieee754_asinl (x)
100     long double x;
101#endif
102{
103  long double t, w, p, q, c, r, s;
104  int32_t ix;
105  u_int32_t se, i0, i1, k;
106
107  GET_LDOUBLE_WORDS (se, i0, i1, x);
108  ix = se & 0x7fff;
109  ix = (ix << 16) | (i0 >> 16);
110  if (ix >= 0x3fff8000)
111    {				/* |x|>= 1 */
112      if (ix == 0x3fff8000 && ((i0 - 0x80000000) | i1) == 0)
113	/* asin(1)=+-pi/2 with inexact */
114	return x * pio2_hi + x * pio2_lo;
115      return (x - x) / (x - x);	/* asin(|x|>1) is NaN */
116    }
117  else if (ix < 0x3ffe8000)
118    {				/* |x|<0.5 */
119      if (ix < 0x3fde8000)
120	{			/* if |x| < 2**-33 */
121	  if (huge + x > one)
122	    return x;		/* return x with inexact if x!=0 */
123	}
124      else
125	{
126	  t = x * x;
127	  p =
128	    t * (pS0 +
129		 t * (pS1 + t * (pS2 + t * (pS3 + t * (pS4 + t * pS5)))));
130	  q = qS0 + t * (qS1 + t * (qS2 + t * (qS3 + t * (qS4 + t))));
131	  w = p / q;
132	  return x + x * w;
133	}
134    }
135  /* 1> |x|>= 0.5 */
136  w = one - fabsl (x);
137  t = w * 0.5;
138  p = t * (pS0 + t * (pS1 + t * (pS2 + t * (pS3 + t * (pS4 + t * pS5)))));
139  q = qS0 + t * (qS1 + t * (qS2 + t * (qS3 + t * (qS4 + t))));
140  s = __ieee754_sqrtl (t);
141  if (ix >= 0x3ffef999)
142    {				/* if |x| > 0.975 */
143      w = p / q;
144      t = pio2_hi - (2.0 * (s + s * w) - pio2_lo);
145    }
146  else
147    {
148      GET_LDOUBLE_WORDS (k, i0, i1, s);
149      i1 = 0;
150      SET_LDOUBLE_WORDS (w,k,i0,i1);
151      c = (t - w * w) / (s + w);
152      r = p / q;
153      p = 2.0 * s * r - (pio2_lo - 2.0 * c);
154      q = pio4_hi - 2.0 * w;
155      t = pio4_hi - (p - q);
156    }
157  if ((se & 0x8000) == 0)
158    return t;
159  else
160    return -t;
161}
162