1 /* origin: FreeBSD /usr/src/lib/msun/src/e_asinl.c */
3 * ====================================================
4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
6 * Developed at SunSoft, a Sun Microsystems, Inc. business.
7 * Permission to use, copy, modify, and distribute this
8 * software is freely granted, provided that this notice
10 * ====================================================
13 * See comments in asin.c.
14 * Converted to long double by David Schultz <das@FreeBSD.ORG>.
19 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
20 long double asinl(long double x
)
24 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
25 #include "__invtrigl.h"
26 #if LDBL_MANT_DIG == 64
27 #define CLOSETO1(u) (u.i.m>>56 >= 0xf7)
28 #define CLEARBOTTOM(u) (u.i.m &= -1ULL << 32)
29 #elif LDBL_MANT_DIG == 113
30 #define CLOSETO1(u) (u.i.top >= 0xee00)
31 #define CLEARBOTTOM(u) (u.i.lo = 0)
34 long double asinl(long double x
)
36 union ldshape u
= {x
};
38 uint16_t e
= u
.i
.se
& 0x7fff;
39 int sign
= u
.i
.se
>> 15;
41 if (e
>= 0x3fff) { /* |x| >= 1 or nan */
42 /* asin(+-1)=+-pi/2 with inexact */
43 if (x
== 1 || x
== -1)
44 return x
*pio2_hi
+ 0x1p
-120f
;
47 if (e
< 0x3fff - 1) { /* |x| < 0.5 */
48 if (e
< 0x3fff - (LDBL_MANT_DIG
+1)/2) {
49 /* return x with inexact if x!=0 */
50 FORCE_EVAL(x
+ 0x1p
120f
);
53 return x
+ x
*__invtrigl_R(x
*x
);
56 z
= (1.0 - fabsl(x
))*0.5;
60 x
= pio2_hi
- (2*(s
+s
*r
)-pio2_lo
);
66 c
= (z
- f
*f
)/(s
+ f
);
67 x
= 0.5*pio2_hi
-(2*s
*r
- (pio2_lo
-2*c
) - (0.5*pio2_hi
-2*f
));