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 long double
70 one = 1.0L,
71 huge = 1.0e+300L,
72 pio2_hi = 1.5707963267948966192313216916397514420986L,
73 pio2_lo = 4.3359050650618905123985220130216759843812E-35L,
74 pio4_hi = 7.8539816339744830961566084581987569936977E-1L,
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 = -8.358099012470680544198472400254596543711E2L,
82 pS1 = 3.674973957689619490312782828051860366493E3L,
83 pS2 = -6.730729094812979665807581609853656623219E3L,
84 pS3 = 6.643843795209060298375552684423454077633E3L,
85 pS4 = -3.817341990928606692235481812252049415993E3L,
86 pS5 = 1.284635388402653715636722822195716476156E3L,
87 pS6 = -2.410736125231549204856567737329112037867E2L,
88 pS7 = 2.219191969382402856557594215833622156220E1L,
89 pS8 = -7.249056260830627156600112195061001036533E-1L,
90 pS9 = 1.055923570937755300061509030361395604448E-3L,
91
92 qS0 = -5.014859407482408326519083440151745519205E3L,
93 qS1 = 2.430653047950480068881028451580393430537E4L,
94 qS2 = -4.997904737193653607449250593976069726962E4L,
95 qS3 = 5.675712336110456923807959930107347511086E4L,
96 qS4 = -3.881523118339661268482937768522572588022E4L,
97 qS5 = 1.634202194895541569749717032234510811216E4L,
98 qS6 = -4.151452662440709301601820849901296953752E3L,
99 qS7 = 5.956050864057192019085175976175695342168E2L,
100 qS8 = -4.175375777334867025769346564600396877176E1L,
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 = -5.619049346208901520945464704848780243887E0L,
107 rS1 = 4.460504162777731472539175700169871920352E1L,
108 rS2 = -1.317669505315409261479577040530751477488E2L,
109 rS3 = 1.626532582423661989632442410808596009227E2L,
110 rS4 = -3.144806644195158614904369445440583873264E1L,
111 rS5 = -9.806674443470740708765165604769099559553E1L,
112 rS6 = 5.708468492052010816555762842394927806920E1L,
113 rS7 = 1.396540499232262112248553357962639431922E1L,
114 rS8 = -1.126243289311910363001762058295832610344E1L,
115 rS9 = -4.956179821329901954211277873774472383512E-1L,
116 rS10 = 3.313227657082367169241333738391762525780E-1L,
117
118 sS0 = -4.645814742084009935700221277307007679325E0L,
119 sS1 = 3.879074822457694323970438316317961918430E1L,
120 sS2 = -1.221986588013474694623973554726201001066E2L,
121 sS3 = 1.658821150347718105012079876756201905822E2L,
122 sS4 = -4.804379630977558197953176474426239748977E1L,
123 sS5 = -1.004296417397316948114344573811562952793E2L,
124 sS6 = 7.530281592861320234941101403870010111138E1L,
125 sS7 = 1.270735595411673647119592092304357226607E1L,
126 sS8 = -1.815144839646376500705105967064792930282E1L,
127 sS9 = -7.821597334910963922204235247786840828217E-2L,
128 /* 1.000000000000000000000000000000000000000E0 */
129
130 asinr5625 = 5.9740641664535021430381036628424864397707E-1L;
131
132
133
134 long double
__ieee754_asinl(long double x)135 __ieee754_asinl (long double x)
136 {
137 long double a, t, w, p, q, c, r, s;
138 int flag;
139
140 if (__glibc_unlikely (isnan (x)))
141 return x + x;
142 flag = 0;
143 a = __builtin_fabsl (x);
144 if (a == 1.0L) /* |x|>= 1 */
145 return x * pio2_hi + x * pio2_lo; /* asin(1)=+-pi/2 with inexact */
146 else if (a >= 1.0L)
147 return (x - x) / (x - x); /* asin(|x|>1) is NaN */
148 else if (a < 0.5L)
149 {
150 if (a < 6.938893903907228e-18L) /* |x| < 2**-57 */
151 {
152 math_check_force_underflow (x);
153 long double 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 (a < 0.625L)
165 {
166 t = a - 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 (x > 0.0L)
192 return t;
193 else
194 return -t;
195 }
196 else
197 {
198 /* 1 > |x| >= 0.625 */
199 w = one - a;
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 = sqrtl (t);
232 if (a > 0.975L)
233 {
234 w = p / q;
235 t = pio2_hi - (2.0 * (s + s * w) - pio2_lo);
236 }
237 else
238 {
239 w = ldbl_high (s);
240 c = (t - w * w) / (s + w);
241 r = p / q;
242 p = 2.0 * s * r - (pio2_lo - 2.0 * c);
243 q = pio4_hi - 2.0 * w;
244 t = pio4_hi - (p - q);
245 }
246
247 if (x > 0.0L)
248 return t;
249 else
250 return -t;
251 }
252 libm_alias_finite (__ieee754_asinl, __asinl)
253