2 * ====================================================
3 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
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
9 * ====================================================
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
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.
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.
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/>. */
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)
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;
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
49 * c = sqrt(z) - f = (z-f*f)/(s+f) ...f+c=sqrt(z)
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))
56 * if x is NaN, return x itself;
57 * if |x|>1, return NaN with invalid signal.
64 #include <math-barriers.h>
65 #include <math_private.h>
66 #include <math-underflow.h>
67 #include <libm-alias-finite.h>
69 static const long double
72 pio2_hi
= 1.5707963267948966192313216916397514420986L,
73 pio2_lo
= 4.3359050650618905123985220130216759843812E-35L,
74 pio4_hi
= 7.8539816339744830961566084581987569936977E-1L,
76 /* coefficient for R(x^2) */
78 /* asin(x) = x + x^3 pS(x^2) / qS(x^2)
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,
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 */
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,
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 */
130 asinr5625
= 5.9740641664535021430381036628424864397707E-1L;
135 __ieee754_asinl (long double x
)
137 long double a
, t
, w
, p
, q
, c
, r
, s
;
140 if (__glibc_unlikely (isnan (x
)))
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 */
147 return (x
- x
) / (x
- x
); /* asin(|x|>1) is NaN */
150 if (a
< 6.938893903907228e-18L) /* |x| < 2**-57 */
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 */
160 /* Mark to use pS, qS later on. */
167 p
= ((((((((((rS10
* t
190 t
= asinr5625
+ p
/ q
;
198 /* 1 > |x| >= 0.625 */
225 if (flag
) /* 2^-57 < |x| < 0.5 */
235 t
= pio2_hi
- (2.0 * (s
+ s
* w
) - pio2_lo
);
240 c
= (t
- w
* w
) / (s
+ w
);
242 p
= 2.0 * s
* r
- (pio2_lo
- 2.0 * c
);
243 q
= pio4_hi
- 2.0 * w
;
244 t
= pio4_hi
- (p
- q
);
252 libm_alias_finite (__ieee754_asinl
, __asinl
)