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 math_check_force_underflow (x);
155 __float128 force_inexact = huge + x;
156 math_force_eval (force_inexact);
157 return x; /* return x with inexact if x!=0 */
158 }
159 else
160 {
161 t = x * x;
162 /* Mark to use pS, qS later on. */
163 flag = 1;
164 }
165 }
166 else if (ix < 0x3ffe4000) /* 0.625 */
167 {
168 t = u.value - 0.5625;
169 p = ((((((((((rS10 * t
170 + rS9) * t
171 + rS8) * t
172 + rS7) * t
173 + rS6) * t
174 + rS5) * t
175 + rS4) * t
176 + rS3) * t
177 + rS2) * t
178 + rS1) * t
179 + rS0) * t;
180
181 q = ((((((((( t
182 + sS9) * t
183 + sS8) * t
184 + sS7) * t
185 + sS6) * t
186 + sS5) * t
187 + sS4) * t
188 + sS3) * t
189 + sS2) * t
190 + sS1) * t
191 + sS0;
192 t = asinr5625 + p / q;
193 if ((sign & 0x80000000) == 0)
194 return t;
195 else
196 return -t;
197 }
198 else
199 {
200 /* 1 > |x| >= 0.625 */
201 w = one - u.value;
202 t = w * 0.5;
203 }
204
205 p = (((((((((pS9 * t
206 + pS8) * t
207 + pS7) * t
208 + pS6) * t
209 + pS5) * t
210 + pS4) * t
211 + pS3) * t
212 + pS2) * t
213 + pS1) * t
214 + pS0) * t;
215
216 q = (((((((( t
217 + qS8) * t
218 + qS7) * t
219 + qS6) * t
220 + qS5) * t
221 + qS4) * t
222 + qS3) * t
223 + qS2) * t
224 + qS1) * t
225 + qS0;
226
227 if (flag) /* 2^-57 < |x| < 0.5 */
228 {
229 w = p / q;
230 return x + x * w;
231 }
232
233 s = sqrtq (t);
234 if (ix >= 0x3ffef333) /* |x| > 0.975 */
235 {
236 w = p / q;
237 t = pio2_hi - (2.0 * (s + s * w) - pio2_lo);
238 }
239 else
240 {
241 u.value = s;
242 u.words32.w3 = 0;
243 u.words32.w2 = 0;
244 w = u.value;
245 c = (t - w * w) / (s + w);
246 r = p / q;
247 p = 2.0 * s * r - (pio2_lo - 2.0 * c);
248 q = pio4_hi - 2.0 * w;
249 t = pio4_hi - (p - q);
250 }
251
252 if ((sign & 0x80000000) == 0)
253 return t;
254 else
255 return -t;
256}
257