diff options
Diffstat (limited to 'lib/msun/src/e_asinl.c')
-rw-r--r-- | lib/msun/src/e_asinl.c | 77 |
1 files changed, 77 insertions, 0 deletions
diff --git a/lib/msun/src/e_asinl.c b/lib/msun/src/e_asinl.c new file mode 100644 index 0000000..a85765f --- /dev/null +++ b/lib/msun/src/e_asinl.c @@ -0,0 +1,77 @@ + +/* @(#)e_asin.c 1.3 95/01/18 */ +/* FreeBSD: head/lib/msun/src/e_asin.c 176451 2008-02-22 02:30:36Z das */ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunSoft, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +#include <sys/cdefs.h> +__FBSDID("$FreeBSD$"); + +/* + * See comments in e_asin.c. + * Converted to long double by David Schultz <das@FreeBSD.ORG>. + */ + +#include <float.h> + +#include "invtrig.h" +#include "math.h" +#include "math_private.h" + +static const long double +one = 1.00000000000000000000e+00, +huge = 1.000e+300; + +long double +asinl(long double x) +{ + union IEEEl2bits u; + long double t=0.0,w,p,q,c,r,s; + int16_t expsign, expt; + u.e = x; + expsign = u.xbits.expsign; + expt = expsign & 0x7fff; + if(expt >= BIAS) { /* |x|>= 1 */ + if(expt==BIAS && ((u.bits.manh&~LDBL_NBIT)|u.bits.manl)==0) + /* asin(1)=+-pi/2 with inexact */ + return x*pio2_hi+x*pio2_lo; + return (x-x)/(x-x); /* asin(|x|>1) is NaN */ + } else if (expt<BIAS-1) { /* |x|<0.5 */ + if(expt<ASIN_LINEAR) { /* if |x| is small, asinl(x)=x */ + if(huge+x>one) return x;/* return x with inexact if x!=0*/ + } + t = x*x; + p = P(t); + q = Q(t); + w = p/q; + return x+x*w; + } + /* 1> |x|>= 0.5 */ + w = one-fabsl(x); + t = w*0.5; + p = P(t); + q = Q(t); + s = sqrtl(t); + if(u.bits.manh>=THRESH) { /* if |x| is close to 1 */ + w = p/q; + t = pio2_hi-(2.0*(s+s*w)-pio2_lo); + } else { + u.e = s; + u.bits.manl = 0; + w = u.e; + c = (t-w*w)/(s+w); + r = p/q; + p = 2.0*s*r-(pio2_lo-2.0*c); + q = pio4_hi-2.0*w; + t = pio4_hi-(p-q); + } + if(expsign>0) return t; else return -t; +} |