1 /* ef_asin.c -- float version of e_asin.c.
2  * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com.
3  */
4 
5 /*
6  * ====================================================
7  * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
8  *
9  * Developed at SunPro, a Sun Microsystems, Inc. business.
10  * Permission to use, copy, modify, and distribute this
11  * software is freely granted, provided that this notice
12  * is preserved.
13  * ====================================================
14  */
15 
16 #include "fdlibm.h"
17 
18 static const float one = 1.0000000000e+00, /* 0x3F800000 */
19     huge = 1.000e+30, pio2_hi = 1.57079637050628662109375f,
20                    pio2_lo = -4.37113900018624283e-8f,
21                    pio4_hi = 0.785398185253143310546875f,
22                    /* coefficient for R(x^2) */
23     pS0 = 1.6666667163e-01, /* 0x3e2aaaab */
24     pS1 = -3.2556581497e-01, /* 0xbea6b090 */
25     pS2 = 2.0121252537e-01, /* 0x3e4e0aa8 */
26     pS3 = -4.0055535734e-02, /* 0xbd241146 */
27     pS4 = 7.9153501429e-04, /* 0x3a4f7f04 */
28     pS5 = 3.4793309169e-05, /* 0x3811ef08 */
29     qS1 = -2.4033949375e+00, /* 0xc019d139 */
30     qS2 = 2.0209457874e+00, /* 0x4001572d */
31     qS3 = -6.8828397989e-01, /* 0xbf303361 */
32     qS4 = 7.7038154006e-02; /* 0x3d9dc62e */
33 
34 float
asinf(float x)35 asinf(float x)
36 {
37     float t, w, p, q, c, r, s;
38     __int32_t hx, ix;
39     GET_FLOAT_WORD(hx, x);
40     ix = hx & 0x7fffffff;
41     if (ix == 0x3f800000) {
42         /* asin(1)=+-pi/2 with inexact */
43         return x * pio2_hi + x * pio2_lo;
44     } else if (ix > 0x3f800000) { /* |x|>= 1 */
45         return __math_invalidf(x); /* asin(|x|>1) is NaN */
46     } else if (ix < 0x3f000000) { /* |x|<0.5 */
47         if (ix < 0x32000000) { /* if |x| < 2**-27 */
48             if (huge + x > one)
49                 return x; /* return x with inexact if x!=0*/
50         } else {
51             t = x * x;
52             p = t *
53                 (pS0 + t * (pS1 + t * (pS2 + t * (pS3 + t * (pS4 + t * pS5)))));
54             q = one + t * (qS1 + t * (qS2 + t * (qS3 + t * qS4)));
55             w = p / q;
56             return x + x * w;
57         }
58     }
59     /* 1> |x|>= 0.5 */
60     w = one - fabsf(x);
61     t = w * (float)0.5;
62     p = t * (pS0 + t * (pS1 + t * (pS2 + t * (pS3 + t * (pS4 + t * pS5)))));
63     q = one + t * (qS1 + t * (qS2 + t * (qS3 + t * qS4)));
64     s = sqrtf(t);
65     if (ix >= 0x3F79999A) { /* if |x| > 0.975 */
66         w = p / q;
67         t = pio2_hi - ((float)2.0 * (s + s * w) - pio2_lo);
68     } else {
69         __int32_t iw;
70         w = s;
71         GET_FLOAT_WORD(iw, w);
72         SET_FLOAT_WORD(w, iw & 0xfffff000);
73         c = (t - w * w) / (s + w);
74         r = p / q;
75         p = (float)2.0 * s * r - (pio2_lo - (float)2.0 * c);
76         q = pio4_hi - (float)2.0 * w;
77         t = pio4_hi - (p - q);
78     }
79     if (hx > 0)
80         return t;
81     else
82         return -t;
83 }
84 
85 _MATH_ALIAS_f_f(asin)
86