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 */ 4 /* 5 * ==================================================== 6 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 7 * 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 11 * is preserved. 12 * ==================================================== 13 */ 14 15 /* 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>. 19 */ 20 21 #include <float.h> 22 #include <math.h> 23 24 #include "invtrig.h" 25 #include "math_private.h" 26 27 #ifdef EXT_IMPLICIT_NBIT 28 #define LDBL_NBIT 0 29 #else /* EXT_IMPLICIT_NBIT */ 30 #define LDBL_NBIT 0x80000000 31 #endif /* EXT_IMPLICIT_NBIT */ 32 33 static const long double 34 one = 1.00000000000000000000e+00, 35 huge = 1.000e+300; 36 37 long double 38 asinl(long double x) 39 { 40 union { 41 long double e; 42 struct ieee_ext bits; 43 } u; 44 long double t=0.0,w,p,q,c,r,s; 45 int16_t expsign, expt; 46 u.e = x; 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) 51 #ifdef EXT_FRACHMBITS 52 | u.bits.ext_frachm 53 #endif /* EXT_FRACHMBITS */ 54 #ifdef EXT_FRACLMBITS 55 | u.bits.ext_fraclm 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*/ 64 } 65 t = x*x; 66 p = P(t); 67 q = Q(t); 68 w = p/q; 69 return x+x*w; 70 } 71 /* 1> |x|>= 0.5 */ 72 w = one-fabsl(x); 73 t = w*0.5; 74 p = P(t); 75 q = Q(t); 76 s = sqrtl(t); 77 #ifdef EXT_FRACHMBITS 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 */ 84 w = p/q; 85 t = pio2_hi-(2.0*(s+s*w)-pio2_lo); 86 } else { 87 u.e = s; 88 u.bits.ext_fracl = 0; 89 #ifdef EXT_FRACLMBITS 90 u.bits.ext_fraclm = 0; 91 #endif /* EXT_FRACLMBITS */ 92 w = u.e; 93 c = (t-w*w)/(s+w); 94 r = p/q; 95 p = 2.0*s*r-(pio2_lo-2.0*c); 96 q = pio4_hi-2.0*w; 97 t = pio4_hi-(p-q); 98 } 99 if(expsign>0) return t; else return -t; 100 } 101