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 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, see
32 <https://www.gnu.org/licenses/>. */
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 * 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 <float.h>
63 #include <math.h>
64 #include <math-barriers.h>
65 #include <math_private.h>
66 #include <math-underflow.h>
67 #include <libm-alias-finite.h>
68
69 static const _Float128
70 one = 1,
71 huge = L(1.0e+4932),
72 pio2_hi = L(1.5707963267948966192313216916397514420986),
73 pio2_lo = L(4.3359050650618905123985220130216759843812E-35),
74 pio4_hi = L(7.8539816339744830961566084581987569936977E-1),
75
76 /* coefficient for R(x^2) */
77
78 /* asin(x) = x + x^3 pS(x^2) / qS(x^2)
79 0 <= x <= 0.5
80 peak relative error 1.9e-35 */
81 pS0 = L(-8.358099012470680544198472400254596543711E2),
82 pS1 = L(3.674973957689619490312782828051860366493E3),
83 pS2 = L(-6.730729094812979665807581609853656623219E3),
84 pS3 = L(6.643843795209060298375552684423454077633E3),
85 pS4 = L(-3.817341990928606692235481812252049415993E3),
86 pS5 = L(1.284635388402653715636722822195716476156E3),
87 pS6 = L(-2.410736125231549204856567737329112037867E2),
88 pS7 = L(2.219191969382402856557594215833622156220E1),
89 pS8 = L(-7.249056260830627156600112195061001036533E-1),
90 pS9 = L(1.055923570937755300061509030361395604448E-3),
91
92 qS0 = L(-5.014859407482408326519083440151745519205E3),
93 qS1 = L(2.430653047950480068881028451580393430537E4),
94 qS2 = L(-4.997904737193653607449250593976069726962E4),
95 qS3 = L(5.675712336110456923807959930107347511086E4),
96 qS4 = L(-3.881523118339661268482937768522572588022E4),
97 qS5 = L(1.634202194895541569749717032234510811216E4),
98 qS6 = L(-4.151452662440709301601820849901296953752E3),
99 qS7 = L(5.956050864057192019085175976175695342168E2),
100 qS8 = L(-4.175375777334867025769346564600396877176E1),
101 /* 1.000000000000000000000000000000000000000E0 */
102
103 /* asin(0.5625 + x) = asin(0.5625) + x rS(x) / sS(x)
104 -0.0625 <= x <= 0.0625
105 peak relative error 3.3e-35 */
106 rS0 = L(-5.619049346208901520945464704848780243887E0),
107 rS1 = L(4.460504162777731472539175700169871920352E1),
108 rS2 = L(-1.317669505315409261479577040530751477488E2),
109 rS3 = L(1.626532582423661989632442410808596009227E2),
110 rS4 = L(-3.144806644195158614904369445440583873264E1),
111 rS5 = L(-9.806674443470740708765165604769099559553E1),
112 rS6 = L(5.708468492052010816555762842394927806920E1),
113 rS7 = L(1.396540499232262112248553357962639431922E1),
114 rS8 = L(-1.126243289311910363001762058295832610344E1),
115 rS9 = L(-4.956179821329901954211277873774472383512E-1),
116 rS10 = L(3.313227657082367169241333738391762525780E-1),
117
118 sS0 = L(-4.645814742084009935700221277307007679325E0),
119 sS1 = L(3.879074822457694323970438316317961918430E1),
120 sS2 = L(-1.221986588013474694623973554726201001066E2),
121 sS3 = L(1.658821150347718105012079876756201905822E2),
122 sS4 = L(-4.804379630977558197953176474426239748977E1),
123 sS5 = L(-1.004296417397316948114344573811562952793E2),
124 sS6 = L(7.530281592861320234941101403870010111138E1),
125 sS7 = L(1.270735595411673647119592092304357226607E1),
126 sS8 = L(-1.815144839646376500705105967064792930282E1),
127 sS9 = L(-7.821597334910963922204235247786840828217E-2),
128 /* 1.000000000000000000000000000000000000000E0 */
129
130 asinr5625 = L(5.9740641664535021430381036628424864397707E-1);
131
132
133
134 _Float128
135 __ieee754_asinl (_Float128 x)
136 {
137 _Float128 t, w, p, q, c, r, s;
138 int32_t ix, sign, flag;
139 ieee854_long_double_shape_type u;
140
141 flag = 0;
142 u.value = x;
143 sign = u.parts32.w0;
144 ix = sign & 0x7fffffff;
145 u.parts32.w0 = ix; /* |x| */
146 if (ix >= 0x3fff0000) /* |x|>= 1 */
147 {
148 if (ix == 0x3fff0000
149 && (u.parts32.w1 | u.parts32.w2 | u.parts32.w3) == 0)
150 /* asin(1)=+-pi/2 with inexact */
151 return x * pio2_hi + x * pio2_lo;
152 return (x - x) / (x - x); /* asin(|x|>1) is NaN */
153 }
154 else if (ix < 0x3ffe0000) /* |x| < 0.5 */
155 {
156 if (ix < 0x3fc60000) /* |x| < 2**-57 */
157 {
158 math_check_force_underflow (x);
159 _Float128 force_inexact = huge + x;
160 math_force_eval (force_inexact);
161 return x; /* return x with inexact if x!=0 */
162 }
163 else
164 {
165 t = x * x;
166 /* Mark to use pS, qS later on. */
167 flag = 1;
168 }
169 }
170 else if (ix < 0x3ffe4000) /* 0.625 */
171 {
172 t = u.value - 0.5625;
173 p = ((((((((((rS10 * t
174 + rS9) * t
175 + rS8) * t
176 + rS7) * t
177 + rS6) * t
178 + rS5) * t
179 + rS4) * t
180 + rS3) * t
181 + rS2) * t
182 + rS1) * t
183 + rS0) * t;
184
185 q = ((((((((( t
186 + sS9) * t
187 + sS8) * t
188 + sS7) * t
189 + sS6) * t
190 + sS5) * t
191 + sS4) * t
192 + sS3) * t
193 + sS2) * t
194 + sS1) * t
195 + sS0;
196 t = asinr5625 + p / q;
197 if ((sign & 0x80000000) == 0)
198 return t;
199 else
200 return -t;
201 }
202 else
203 {
204 /* 1 > |x| >= 0.625 */
205 w = one - u.value;
206 t = w * 0.5;
207 }
208
209 p = (((((((((pS9 * t
210 + pS8) * t
211 + pS7) * t
212 + pS6) * t
213 + pS5) * t
214 + pS4) * t
215 + pS3) * t
216 + pS2) * t
217 + pS1) * t
218 + pS0) * t;
219
220 q = (((((((( t
221 + qS8) * t
222 + qS7) * t
223 + qS6) * t
224 + qS5) * t
225 + qS4) * t
226 + qS3) * t
227 + qS2) * t
228 + qS1) * t
229 + qS0;
230
231 if (flag) /* 2^-57 < |x| < 0.5 */
232 {
233 w = p / q;
234 return x + x * w;
235 }
236
237 s = sqrtl (t);
238 if (ix >= 0x3ffef333) /* |x| > 0.975 */
239 {
240 w = p / q;
241 t = pio2_hi - (2.0 * (s + s * w) - pio2_lo);
242 }
243 else
244 {
245 u.value = s;
246 u.parts32.w3 = 0;
247 u.parts32.w2 = 0;
248 w = u.value;
249 c = (t - w * w) / (s + w);
250 r = p / q;
251 p = 2.0 * s * r - (pio2_lo - 2.0 * c);
252 q = pio4_hi - 2.0 * w;
253 t = pio4_hi - (p - q);
254 }
255
256 if ((sign & 0x80000000) == 0)
257 return t;
258 else
259 return -t;
260 }
261 libm_alias_finite (__ieee754_asinl, __asinl)