| author | Rich Felker <dalias@aerifal.cx> | 2012年03月13日 01:17:53 -0400 |
|---|---|---|
| committer | Rich Felker <dalias@aerifal.cx> | 2012年03月13日 01:17:53 -0400 |
| commit | b69f695acedd4ce2798ef9ea28d834ceccc789bd (patch) | |
| tree | eafd98b9b75160210f3295ac074d699f863d958e /src/math/asinf.c | |
| parent | d46cf2e14cc4df7cc75e77d7009fcb6df1f48a33 (diff) | |
| download | musl-b69f695acedd4ce2798ef9ea28d834ceccc789bd.tar.gz | |
| -rw-r--r-- | src/math/asinf.c | 64 |
diff --git a/src/math/asinf.c b/src/math/asinf.c new file mode 100644 index 00000000..729dd37f --- /dev/null +++ b/src/math/asinf.c @@ -0,0 +1,64 @@ +/* origin: FreeBSD /usr/src/lib/msun/src/e_asinf.c */ +/* + * Conversion to float by Ian Lance Taylor, Cygnus Support, ian@cygnus.com. + */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#include "libm.h" + +static const float +one = 1.0000000000e+00, /* 0x3F800000 */ +huge = 1.000e+30, +/* coefficients for R(x^2) */ +pS0 = 1.6666586697e-01, +pS1 = -4.2743422091e-02, +pS2 = -8.6563630030e-03, +qS1 = -7.0662963390e-01; + +static const double +pio2 = 1.570796326794896558e+00; + +float asinf(float x) +{ + double s; + float t,w,p,q; + int32_t hx,ix; + + GET_FLOAT_WORD(hx, x); + ix = hx & 0x7fffffff; + if (ix >= 0x3f800000) { /* |x| >= 1 */ + if (ix == 0x3f800000) /* |x| == 1 */ + return x*pio2; /* asin(+-1) = +-pi/2 with inexact */ + return (x-x)/(x-x); /* asin(|x|>1) is NaN */ + } else if (ix < 0x3f000000) { /* |x|<0.5 */ + if (ix < 0x39800000) { /* |x| < 2**-12 */ + if (huge+x > one) + return x; /* return x with inexact if x!=0 */ + } + t = x*x; + p = t*(pS0+t*(pS1+t*pS2)); + q = one+t*qS1; + w = p/q; + return x + x*w; + } + /* 1 > |x| >= 0.5 */ + w = one - fabsf(x); + t = w*(float)0.5; + p = t*(pS0+t*(pS1+t*pS2)); + q = one+t*qS1; + s = sqrt(t); + w = p/q; + t = pio2-2.0*(s+s*w); + if (hx > 0) + return t; + return -t; +} |