OLD | NEW |
(Empty) | |
| 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 |
OLD | NEW |