1 /* origin: FreeBSD /usr/src/lib/msun/src/e_acosl.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 acos.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 acosl(long double x) 21 { 22 return acos(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 CLEARBOTTOM(u) (u.i.m &= -1ULL << 32) 28 #elif LDBL_MANT_DIG == 113 29 #define CLEARBOTTOM(u) (u.i.lo = 0) 30 #endif 31 32 long double acosl(long double x) 33 { 34 union ldshape u = {x}; 35 long double z, s, c, f; 36 uint16_t e = u.i.se & 0x7fff; 37 38 /* |x| >= 1 or nan */ 39 if (e >= 0x3fff) { 40 if (x == 1) 41 return 0; 42 if (x == -1) 43 return 2*pio2_hi + 0x1p-120f; 44 return 0/(x-x); 45 } 46 /* |x| < 0.5 */ 47 if (e < 0x3fff - 1) { 48 if (e < 0x3fff - LDBL_MANT_DIG - 1) 49 return pio2_hi + 0x1p-120f; 50 return pio2_hi - (__invtrigl_R(x*x)*x - pio2_lo + x); 51 } 52 /* x < -0.5 */ 53 if (u.i.se >> 15) { 54 z = (1 + x)*0.5; 55 s = sqrtl(z); 56 return 2*(pio2_hi - (__invtrigl_R(z)*s - pio2_lo + s)); 57 } 58 /* x > 0.5 */ 59 z = (1 - x)*0.5; 60 s = sqrtl(z); 61 u.f = s; 62 CLEARBOTTOM(u); 63 f = u.f; 64 c = (z - f*f)/(s + f); 65 return 2*(__invtrigl_R(z)*s + c + f); 66 } 67 #endif 68