math: add fma TODO comments about the underflow issue
[oweals/musl.git] / src / math / asinl.c
1 /* origin: FreeBSD /usr/src/lib/msun/src/e_asinl.c */
2 /*
3  * ====================================================
4  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
5  *
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
9  * is preserved.
10  * ====================================================
11  */
12 /*
13  * See comments in asin.c.
14  * Converted to long double by David Schultz <das@FreeBSD.ORG>.
15  */
16
17 #include "libm.h"
18
19 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024
20 long double asinl(long double x)
21 {
22         return asin(x);
23 }
24 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384
25 #include "__invtrigl.h"
26 /* 0.95 */
27 #define THRESH  ((0xe666666666666666ULL>>(64-(LDBL_MANH_SIZE-1)))|LDBL_NBIT)
28
29 long double asinl(long double x)
30 {
31         union IEEEl2bits u;
32         long double z,r,s;
33         uint16_t expsign, expt;
34
35         u.e = x;
36         expsign = u.xbits.expsign;
37         expt = expsign & 0x7fff;
38         if (expt >= 0x3fff) {   /* |x| >= 1 or nan */
39                 if (expt == 0x3fff &&
40                     ((u.bits.manh&~LDBL_NBIT)|u.bits.manl) == 0)
41                         /* asin(+-1)=+-pi/2 with inexact */
42                         return x*pio2_hi + 0x1p-120f;
43                 return 0/(x-x);
44         }
45         if (expt < 0x3fff - 1) {  /* |x| < 0.5 */
46                 if (expt < 0x3fff - 32) {  /* |x|<0x1p-32, asinl(x)=x */
47                         /* return x with inexact if x!=0 */
48                         FORCE_EVAL(x + 0x1p120f);
49                         return x;
50                 }
51                 return x + x*__invtrigl_R(x*x);
52         }
53         /* 1 > |x| >= 0.5 */
54         z = (1.0 - fabsl(x))*0.5;
55         s = sqrtl(z);
56         r = __invtrigl_R(z);
57         if (u.bits.manh >= THRESH) { /* if |x| is close to 1 */
58                 x = pio2_hi - (2*(s+s*r)-pio2_lo);
59         } else {
60                 long double f, c;
61                 u.e = s;
62                 u.bits.manl = 0;
63                 f = u.e;
64                 c = (z-f*f)/(s+f);
65                 x = 0.5*pio2_hi-(2*s*r - (pio2_lo-2*c) - (0.5*pio2_hi-2*f));
66         }
67         if (expsign>>15)
68                 return -x;
69         return x;
70 }
71 #endif