From b12a73d5bf595b7fbb73db30c6bb144078e86ef5 Mon Sep 17 00:00:00 2001 From: Szabolcs Nagy Date: Tue, 11 Dec 2012 23:56:59 +0100 Subject: math: clean up inverse trigonometric functions modifications: * avoid unsigned->signed conversions * removed various volatile hacks * use FORCE_EVAL when evaluating only for side-effects * factor out R() rational approximation instead of manual inline * __invtrigl.h now only provides __invtrigl_R, __pio2_hi and __pio2_lo * use 2*pio2_hi, 2*pio2_lo instead of pi_hi, pi_lo otherwise the logic is not changed, long double versions will need a revisit when a genaral long double cleanup happens --- src/math/asinf.c | 51 +++++++++++++++++++++++++-------------------------- 1 file changed, 25 insertions(+), 26 deletions(-) (limited to 'src/math/asinf.c') diff --git a/src/math/asinf.c b/src/math/asinf.c index c1c42c7b..462bf043 100644 --- a/src/math/asinf.c +++ b/src/math/asinf.c @@ -12,52 +12,51 @@ * is preserved. * ==================================================== */ - #include "libm.h" +static const double +pio2 = 1.570796326794896558e+00; + static const float -huge = 1.000e+30, /* coefficients for R(x^2) */ pS0 = 1.6666586697e-01, pS1 = -4.2743422091e-02, pS2 = -8.6563630030e-03, qS1 = -7.0662963390e-01; -static const double -pio2 = 1.570796326794896558e+00; +static float R(float z) +{ + float p, q; + p = z*(pS0+z*(pS1+z*pS2)); + q = 1.0f+z*qS1; + return p/q; +} float asinf(float x) { double s; - float t,w,p,q; - int32_t hx,ix; + float z; + uint32_t hx,ix; GET_FLOAT_WORD(hx, x); ix = hx & 0x7fffffff; if (ix >= 0x3f800000) { /* |x| >= 1 */ if (ix == 0x3f800000) /* |x| == 1 */ - return x*pio2; /* asin(+-1) = +-pi/2 with inexact */ - return (x-x)/(x-x); /* asin(|x|>1) is NaN */ - } else if (ix < 0x3f000000) { /* |x|<0.5 */ + return x*pio2 + 0x1p-120f; /* asin(+-1) = +-pi/2 with inexact */ + return 0/(x-x); /* asin(|x|>1) is NaN */ + } + if (ix < 0x3f000000) { /* |x| < 0.5 */ if (ix < 0x39800000) { /* |x| < 2**-12 */ - if (huge+x > 1.0f) - return x; /* return x with inexact if x!=0 */ + FORCE_EVAL(x + 0x1p120f); + return x; /* return x with inexact if x!=0 */ } - t = x*x; - p = t*(pS0+t*(pS1+t*pS2)); - q = 1.0f+t*qS1; - w = p/q; - return x + x*w; + return x + x*R(x*x); } /* 1 > |x| >= 0.5 */ - w = 1.0f - fabsf(x); - t = w*0.5f; - p = t*(pS0+t*(pS1+t*pS2)); - q = 1.0f+t*qS1; - s = sqrt(t); - w = p/q; - t = pio2-2.0*(s+s*w); - if (hx > 0) - return t; - return -t; + z = (1 - fabsf(x))*0.5f; + s = sqrt(z); + x = pio2 - 2*(s+s*R(z)); + if (hx >> 31) + return -x; + return x; } -- cgit v1.2.1