1 /* origin: FreeBSD /usr/src/lib/msun/src/e_asinl.c */ 2 /* 3 * ==================================================== 4 * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. 5 * 6 * Developed at SunSoft, a Sun Microsystems, Inc. business. 7 * Permission to use, copy, modify, and distribute this 8 * software is freely granted, provided that this notice 9 * is preserved. 10 * ==================================================== 11 */ 12 /* 13 * See comments in asin.c. 14 * Converted to long double by David Schultz <das@FreeBSD.ORG>. 15 */ 16 17 #include "libm.h" 18 19 #if LDBL_MANT_DIG == 53 && LDBL_MAX_EXP == 1024 20 long double asinl(long double x) 21 { 22 return asin(x); 23 } 24 #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 25 #include "__invtrigl.h" 26 #if LDBL_MANT_DIG == 64 27 #define CLOSETO1(u) (u.i.m>>56 >= 0xf7) 28 #define CLEARBOTTOM(u) (u.i.m &= -1ULL << 32) 29 #elif LDBL_MANT_DIG == 113 30 #define CLOSETO1(u) (u.i.top >= 0xee00) 31 #define CLEARBOTTOM(u) (u.i.lo = 0) 32 #endif 33 34 long double asinl(long double x) 35 { 36 union ldshape u = {x}; 37 long double z, r, s; 38 uint16_t e = u.i.se & 0x7fff; 39 int sign = u.i.se >> 15; 40 41 if (e >= 0x3fff) { /* |x| >= 1 or nan */ 42 /* asin(+-1)=+-pi/2 with inexact */ 43 if (x == 1 || x == -1) 44 return x*pio2_hi + 0x1p-120f; 45 return 0/(x-x); 46 } 47 if (e < 0x3fff - 1) { /* |x| < 0.5 */ 48 if (e < 0x3fff - (LDBL_MANT_DIG+1)/2) { 49 /* return x with inexact if x!=0 */ 50 FORCE_EVAL(x + 0x1p120f); 51 return x; 52 } 53 return x + x*__invtrigl_R(x*x); 54 } 55 /* 1 > |x| >= 0.5 */ 56 z = (1.0 - fabsl(x))*0.5; 57 s = sqrtl(z); 58 r = __invtrigl_R(z); 59 if (CLOSETO1(u)) { 60 x = pio2_hi - (2*(s+s*r)-pio2_lo); 61 } else { 62 long double f, c; 63 u.f = s; 64 CLEARBOTTOM(u); 65 f = u.f; 66 c = (z - f*f)/(s + f); 67 x = 0.5*pio2_hi-(2*s*r - (pio2_lo-2*c) - (0.5*pio2_hi-2*f)); 68 } 69 return sign ? -x : x; 70 } 71 #endif 72