1 /* $OpenBSD: e_asinl.c,v 1.1 2008/12/09 20:00:35 martynas Exp $ */
2 /* @(#)e_asin.c 1.3 95/01/18 */
3 /* FreeBSD: head/lib/msun/src/e_asin.c 176451 2008-02-22 02:30:36Z das */
5 * ====================================================
6 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved.
8 * Developed at SunSoft, a Sun Microsystems, Inc. business.
9 * Permission to use, copy, modify, and distribute this
10 * software is freely granted, provided that this notice
12 * ====================================================
16 * See comments in e_asin.c.
17 * Converted to long double by David Schultz <das@FreeBSD.ORG>.
18 * Adapted for OpenBSD by Martynas Venckus <martynas@openbsd.org>.
25 #include "math_private.h"
27 #ifdef EXT_IMPLICIT_NBIT
29 #else /* EXT_IMPLICIT_NBIT */
30 #define LDBL_NBIT 0x80000000
31 #endif /* EXT_IMPLICIT_NBIT */
33 static const long double
34 one = 1.00000000000000000000e+00,
44 long double t=0.0,w,p,q,c,r,s;
45 int16_t expsign, expt;
47 expsign = (u.bits.ext_sign << 15) | u.bits.ext_exp;
48 expt = expsign & 0x7fff;
49 if(expt >= BIAS) { /* |x|>= 1 */
50 if(expt==BIAS && ((u.bits.ext_frach&~LDBL_NBIT)
53 #endif /* EXT_FRACHMBITS */
56 #endif /* EXT_FRACLMBITS */
57 | u.bits.ext_fracl)==0)
58 /* asin(1)=+-pi/2 with inexact */
59 return x*pio2_hi+x*pio2_lo;
60 return (x-x)/(x-x); /* asin(|x|>1) is NaN */
61 } else if (expt<BIAS-1) { /* |x|<0.5 */
62 if(expt<ASIN_LINEAR) { /* if |x| is small, asinl(x)=x */
63 if(huge+x>one) return x;/* return x with inexact if x!=0*/
78 if((((uint64_t)u.bits.ext_frach << EXT_FRACHMBITS)
79 | u.bits.ext_frachm) >= THRESH) {
80 /* if |x| is close to 1 */
81 #else /* EXT_FRACHMBITS */
82 if(u.bits.ext_frach>=THRESH) { /* if |x| is close to 1 */
83 #endif /* EXT_FRACHMBITS */
85 t = pio2_hi-(2.0*(s+s*w)-pio2_lo);
90 u.bits.ext_fraclm = 0;
91 #endif /* EXT_FRACLMBITS */
95 p = 2.0*s*r-(pio2_lo-2.0*c);
99 if(expsign>0) return t; else return -t;