X-Git-Url: http://nsz.repo.hu/git/?a=blobdiff_plain;f=src%2Fmath%2Facosl.c;h=9e7b7fb3a146c0c978070cd4b9eb428d873a0d21;hb=ae4b0b96d63fe3cbd70008350f998570c9e91f7f;hp=21e6c95e8cba60bdeaa278a1c814d8418c1d10bd;hpb=b69f695acedd4ce2798ef9ea28d834ceccc789bd;p=musl diff --git a/src/math/acosl.c b/src/math/acosl.c index 21e6c95e..9e7b7fb3 100644 --- a/src/math/acosl.c +++ b/src/math/acosl.c @@ -24,68 +24,45 @@ long double acosl(long double x) #elif (LDBL_MANT_DIG == 64 || LDBL_MANT_DIG == 113) && LDBL_MAX_EXP == 16384 #include "__invtrigl.h" -static const long double -one = 1.00000000000000000000e+00; - -// FIXME -//#ifdef __i386__ -/* XXX Work around the fact that gcc truncates long double constants on i386 */ -static volatile double -pi1 = 3.14159265358979311600e+00, /* 0x1.921fb54442d18p+1 */ -pi2 = 1.22514845490862001043e-16; /* 0x1.1a80000000000p-53 */ -#define pi ((long double)pi1 + pi2) -//#else -#if 0 -static const long double -pi = 3.14159265358979323846264338327950280e+00L; -#endif - long double acosl(long double x) { union IEEEl2bits u; - long double z, p, q, r, w, s, c, df; + long double z, w, s, c, df; int16_t expsign, expt; u.e = x; expsign = u.xbits.expsign; expt = expsign & 0x7fff; - if (expt >= BIAS) { /* |x| >= 1 */ - if (expt == BIAS && + /* |x| >= 1 or nan */ + if (expt >= 0x3fff) { + if (expt == 0x3fff && ((u.bits.manh & ~LDBL_NBIT) | u.bits.manl) == 0) { if (expsign > 0) - return 0.0; /* acos(1) = 0 */ - else - return pi + 2.0 * pio2_lo; /* acos(-1)= pi */ + return 0; /* acos(1) = 0 */ + return 2*pio2_hi + 0x1p-120f; /* acos(-1)= pi */ } - return (x - x) / (x - x); /* acos(|x|>1) is NaN */ + return 0/(x-x); /* acos(|x|>1) is NaN */ } - if (expt < BIAS - 1) { /* |x| < 0.5 */ - if (expt < ACOS_CONST) - return pio2_hi + pio2_lo; /* x tiny: acosl=pi/2 */ - z = x * x; - p = P(z); - q = Q(z); - r = p / q; - return pio2_hi - (x - (pio2_lo - x * r)); - } else if (expsign < 0) { /* x < -0.5 */ - z = (one + x) * 0.5; - p = P(z); - q = Q(z); - s = sqrtl(z); - r = p / q; - w = r * s - pio2_lo; - return pi - 2.0 * (s + w); - } else { /* x > 0.5 */ - z = (one - x) * 0.5; + /* |x| < 0.5 */ + if (expt < 0x3fff - 1) { + if (expt < 0x3fff - 65) + return pio2_hi + 0x1p-120f; /* x < 0x1p-65: acosl(x)=pi/2 */ + return pio2_hi - (x - (pio2_lo - x * __invtrigl_R(x*x))); + } + /* x < -0.5 */ + if (expsign < 0) { + z = (1.0 + x) * 0.5; s = sqrtl(z); - u.e = s; - u.bits.manl = 0; - df = u.e; - c = (z - df * df) / (s + df); - p = P(z); - q = Q(z); - r = p / q; - w = r * s + c; - return 2.0 * (df + w); + w = __invtrigl_R(z) * s - pio2_lo; + return 2*(pio2_hi - (s + w)); } + /* x > 0.5 */ + z = (1.0 - x) * 0.5; + s = sqrtl(z); + u.e = s; + u.bits.manl = 0; + df = u.e; + c = (z - df * df) / (s + df); + w = __invtrigl_R(z) * s + c; + return 2*(df + w); } #endif