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 <http://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 #include "quadmath-imp.h"
62
63 static const __float128
64 one = 1,
65 huge = 1.0e+4932Q,
66 pio2_hi = 1.5707963267948966192313216916397514420986Q,
67 pio2_lo = 4.3359050650618905123985220130216759843812E-35Q,
68 pio4_hi = 7.8539816339744830961566084581987569936977E-1Q,
69
70 /* coefficient for R(x^2) */
71
72 /* asin(x) = x + x^3 pS(x^2) / qS(x^2)
73 0 <= x <= 0.5
74 peak relative error 1.9e-35 */
75 pS0 = -8.358099012470680544198472400254596543711E2Q,
76 pS1 = 3.674973957689619490312782828051860366493E3Q,
77 pS2 = -6.730729094812979665807581609853656623219E3Q,
78 pS3 = 6.643843795209060298375552684423454077633E3Q,
79 pS4 = -3.817341990928606692235481812252049415993E3Q,
80 pS5 = 1.284635388402653715636722822195716476156E3Q,
81 pS6 = -2.410736125231549204856567737329112037867E2Q,
82 pS7 = 2.219191969382402856557594215833622156220E1Q,
83 pS8 = -7.249056260830627156600112195061001036533E-1Q,
84 pS9 = 1.055923570937755300061509030361395604448E-3Q,
85
86 qS0 = -5.014859407482408326519083440151745519205E3Q,
87 qS1 = 2.430653047950480068881028451580393430537E4Q,
88 qS2 = -4.997904737193653607449250593976069726962E4Q,
89 qS3 = 5.675712336110456923807959930107347511086E4Q,
90 qS4 = -3.881523118339661268482937768522572588022E4Q,
91 qS5 = 1.634202194895541569749717032234510811216E4Q,
92 qS6 = -4.151452662440709301601820849901296953752E3Q,
93 qS7 = 5.956050864057192019085175976175695342168E2Q,
94 qS8 = -4.175375777334867025769346564600396877176E1Q,
95 /* 1.000000000000000000000000000000000000000E0 */
96
97 /* asin(0.5625 + x) = asin(0.5625) + x rS(x) / sS(x)
98 -0.0625 <= x <= 0.0625
99 peak relative error 3.3e-35 */
100 rS0 = -5.619049346208901520945464704848780243887E0Q,
101 rS1 = 4.460504162777731472539175700169871920352E1Q,
102 rS2 = -1.317669505315409261479577040530751477488E2Q,
103 rS3 = 1.626532582423661989632442410808596009227E2Q,
104 rS4 = -3.144806644195158614904369445440583873264E1Q,
105 rS5 = -9.806674443470740708765165604769099559553E1Q,
106 rS6 = 5.708468492052010816555762842394927806920E1Q,
107 rS7 = 1.396540499232262112248553357962639431922E1Q,
108 rS8 = -1.126243289311910363001762058295832610344E1Q,
109 rS9 = -4.956179821329901954211277873774472383512E-1Q,
110 rS10 = 3.313227657082367169241333738391762525780E-1Q,
111
112 sS0 = -4.645814742084009935700221277307007679325E0Q,
113 sS1 = 3.879074822457694323970438316317961918430E1Q,
114 sS2 = -1.221986588013474694623973554726201001066E2Q,
115 sS3 = 1.658821150347718105012079876756201905822E2Q,
116 sS4 = -4.804379630977558197953176474426239748977E1Q,
117 sS5 = -1.004296417397316948114344573811562952793E2Q,
118 sS6 = 7.530281592861320234941101403870010111138E1Q,
119 sS7 = 1.270735595411673647119592092304357226607E1Q,
120 sS8 = -1.815144839646376500705105967064792930282E1Q,
121 sS9 = -7.821597334910963922204235247786840828217E-2Q,
122 /* 1.000000000000000000000000000000000000000E0 */
123
124 asinr5625 = 5.9740641664535021430381036628424864397707E-1Q;
125
126
127
128 __float128
129 asinq (__float128 x)
130 {
131 __float128 t, w, p, q, c, r, s;
132 int32_t ix, sign, flag;
133 ieee854_float128 u;
134
135 flag = 0;
136 u.value = x;
137 sign = u.words32.w0;
138 ix = sign & 0x7fffffff;
139 u.words32.w0 = ix; /* |x| */
140 if (ix >= 0x3fff0000) /* |x|>= 1 */
141 {
142 if (ix == 0x3fff0000
143 && (u.words32.w1 | u.words32.w2 | u.words32.w3) == 0)
144 /* asin(1)=+-pi/2 with inexact */
145 return x * pio2_hi + x * pio2_lo;
146 return (x - x) / (x - x); /* asin(|x|>1) is NaN */
147 }
148 else if (ix < 0x3ffe0000) /* |x| < 0.5 */
149 {
150 if (ix < 0x3fc60000) /* |x| < 2**-57 */
151 {
152 math_check_force_underflow (x);
153 __float128 force_inexact = huge + x;
154 math_force_eval (force_inexact);
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 }