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  __float128 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 the
18  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/* asinq(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 *      Between .5 and .625 the approximation is
40 *              asin(0.5625 + x) = asin(0.5625) + x rS(x) / sS(x)
41 *	For x in [0.625,1]
42 *		asin(x) = pi/2-2*asin(sqrt((1-x)/2))
43 *	Let y = (1-x), z = y/2, s := sqrt(z), and pio2_hi+pio2_lo=pi/2;
44 *	then for x>0.98
45 *		asin(x) = pi/2 - 2*(s+s*z*R(z))
46 *			= pio2_hi - (2*(s+s*z*R(z)) - pio2_lo)
47 *	For x<=0.98, let pio4_hi = pio2_hi/2, then
48 *		f = hi part of s;
49 *		c = sqrt(z) - f = (z-f*f)/(s+f) 	...f+c=sqrt(z)
50 *	and
51 *		asin(x) = pi/2 - 2*(s+s*z*R(z))
52 *			= pio4_hi+(pio4-2s)-(2s*z*R(z)-pio2_lo)
53 *			= pio4_hi+(pio4-2f)-(2s*z*R(z)-(pio2_lo+2c))
54 *
55 * Special cases:
56 *	if x is NaN, return x itself;
57 *	if |x|>1, return NaN with invalid signal.
58 *
59 */
60
61
62#include "quadmath-imp.h"
63
64static const __float128
65  one = 1.0Q,
66  huge = 1.0e+4932Q,
67  pio2_hi = 1.5707963267948966192313216916397514420986Q,
68  pio2_lo = 4.3359050650618905123985220130216759843812E-35Q,
69  pio4_hi = 7.8539816339744830961566084581987569936977E-1Q,
70
71	/* coefficient for R(x^2) */
72
73  /* asin(x) = x + x^3 pS(x^2) / qS(x^2)
74     0 <= x <= 0.5
75     peak relative error 1.9e-35  */
76  pS0 = -8.358099012470680544198472400254596543711E2Q,
77  pS1 =  3.674973957689619490312782828051860366493E3Q,
78  pS2 = -6.730729094812979665807581609853656623219E3Q,
79  pS3 =  6.643843795209060298375552684423454077633E3Q,
80  pS4 = -3.817341990928606692235481812252049415993E3Q,
81  pS5 =  1.284635388402653715636722822195716476156E3Q,
82  pS6 = -2.410736125231549204856567737329112037867E2Q,
83  pS7 =  2.219191969382402856557594215833622156220E1Q,
84  pS8 = -7.249056260830627156600112195061001036533E-1Q,
85  pS9 =  1.055923570937755300061509030361395604448E-3Q,
86
87  qS0 = -5.014859407482408326519083440151745519205E3Q,
88  qS1 =  2.430653047950480068881028451580393430537E4Q,
89  qS2 = -4.997904737193653607449250593976069726962E4Q,
90  qS3 =  5.675712336110456923807959930107347511086E4Q,
91  qS4 = -3.881523118339661268482937768522572588022E4Q,
92  qS5 =  1.634202194895541569749717032234510811216E4Q,
93  qS6 = -4.151452662440709301601820849901296953752E3Q,
94  qS7 =  5.956050864057192019085175976175695342168E2Q,
95  qS8 = -4.175375777334867025769346564600396877176E1Q,
96  /* 1.000000000000000000000000000000000000000E0 */
97
98  /* asin(0.5625 + x) = asin(0.5625) + x rS(x) / sS(x)
99     -0.0625 <= x <= 0.0625
100     peak relative error 3.3e-35  */
101  rS0 = -5.619049346208901520945464704848780243887E0Q,
102  rS1 =  4.460504162777731472539175700169871920352E1Q,
103  rS2 = -1.317669505315409261479577040530751477488E2Q,
104  rS3 =  1.626532582423661989632442410808596009227E2Q,
105  rS4 = -3.144806644195158614904369445440583873264E1Q,
106  rS5 = -9.806674443470740708765165604769099559553E1Q,
107  rS6 =  5.708468492052010816555762842394927806920E1Q,
108  rS7 =  1.396540499232262112248553357962639431922E1Q,
109  rS8 = -1.126243289311910363001762058295832610344E1Q,
110  rS9 = -4.956179821329901954211277873774472383512E-1Q,
111  rS10 =  3.313227657082367169241333738391762525780E-1Q,
112
113  sS0 = -4.645814742084009935700221277307007679325E0Q,
114  sS1 =  3.879074822457694323970438316317961918430E1Q,
115  sS2 = -1.221986588013474694623973554726201001066E2Q,
116  sS3 =  1.658821150347718105012079876756201905822E2Q,
117  sS4 = -4.804379630977558197953176474426239748977E1Q,
118  sS5 = -1.004296417397316948114344573811562952793E2Q,
119  sS6 =  7.530281592861320234941101403870010111138E1Q,
120  sS7 =  1.270735595411673647119592092304357226607E1Q,
121  sS8 = -1.815144839646376500705105967064792930282E1Q,
122  sS9 = -7.821597334910963922204235247786840828217E-2Q,
123  /*  1.000000000000000000000000000000000000000E0 */
124
125 asinr5625 =  5.9740641664535021430381036628424864397707E-1Q;
126
127
128
129__float128
130asinq (__float128 x)
131{
132  __float128 t = 0;
133  __float128 w, p, q, c, r, s;
134  int32_t ix, sign, flag;
135  ieee854_float128 u;
136
137  flag = 0;
138  u.value = x;
139  sign = u.words32.w0;
140  ix = sign & 0x7fffffff;
141  u.words32.w0 = ix;    /* |x| */
142  if (ix >= 0x3fff0000)	/* |x|>= 1 */
143    {
144      if (ix == 0x3fff0000
145	  && (u.words32.w1 | u.words32.w2 | u.words32.w3) == 0)
146	/* asin(1)=+-pi/2 with inexact */
147	return x * pio2_hi + x * pio2_lo;
148      return (x - x) / (x - x);	/* asin(|x|>1) is NaN */
149    }
150  else if (ix < 0x3ffe0000) /* |x| < 0.5 */
151    {
152      if (ix < 0x3fc60000) /* |x| < 2**-57 */
153	{
154	  if (huge + x > one)
155	    return x;		/* return x with inexact if x!=0 */
156	}
157      else
158	{
159	  t = x * x;
160	  /* Mark to use pS, qS later on.  */
161	  flag = 1;
162	}
163    }
164  else if (ix < 0x3ffe4000) /* 0.625 */
165    {
166      t = u.value - 0.5625;
167      p = ((((((((((rS10 * t
168		    + rS9) * t
169		   + rS8) * t
170		  + rS7) * t
171		 + rS6) * t
172		+ rS5) * t
173	       + rS4) * t
174	      + rS3) * t
175	     + rS2) * t
176	    + rS1) * t
177	   + rS0) * t;
178
179      q = ((((((((( t
180		    + sS9) * t
181		  + sS8) * t
182		 + sS7) * t
183		+ sS6) * t
184	       + sS5) * t
185	      + sS4) * t
186	     + sS3) * t
187	    + sS2) * t
188	   + sS1) * t
189	+ sS0;
190      t = asinr5625 + p / q;
191      if ((sign & 0x80000000) == 0)
192	return t;
193      else
194	return -t;
195    }
196  else
197    {
198      /* 1 > |x| >= 0.625 */
199      w = one - u.value;
200      t = w * 0.5;
201    }
202
203  p = (((((((((pS9 * t
204	       + pS8) * t
205	      + pS7) * t
206	     + pS6) * t
207	    + pS5) * t
208	   + pS4) * t
209	  + pS3) * t
210	 + pS2) * t
211	+ pS1) * t
212       + pS0) * t;
213
214  q = (((((((( t
215	      + qS8) * t
216	     + qS7) * t
217	    + qS6) * t
218	   + qS5) * t
219	  + qS4) * t
220	 + qS3) * t
221	+ qS2) * t
222       + qS1) * t
223    + qS0;
224
225  if (flag) /* 2^-57 < |x| < 0.5 */
226    {
227      w = p / q;
228      return x + x * w;
229    }
230
231  s = sqrtq (t);
232  if (ix >= 0x3ffef333) /* |x| > 0.975 */
233    {
234      w = p / q;
235      t = pio2_hi - (2.0 * (s + s * w) - pio2_lo);
236    }
237  else
238    {
239      u.value = s;
240      u.words32.w3 = 0;
241      u.words32.w2 = 0;
242      w = u.value;
243      c = (t - w * w) / (s + w);
244      r = p / q;
245      p = 2.0 * s * r - (pio2_lo - 2.0 * c);
246      q = pio4_hi - 2.0 * w;
247      t = pio4_hi - (p - q);
248    }
249
250  if ((sign & 0x80000000) == 0)
251    return t;
252  else
253    return -t;
254}
255