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 20long 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 32long 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