diff options
author | Ulrich Drepper <drepper@redhat.com> | 2001-04-17 06:51:57 +0000 |
---|---|---|
committer | Ulrich Drepper <drepper@redhat.com> | 2001-04-17 06:51:57 +0000 |
commit | 9b7ee67e0c61b046f86062ef8e2ba91f352ebcc1 (patch) | |
tree | dcecb58203979557aae009d73f6bae715029e740 /sysdeps/ieee754 | |
parent | c991a86a175546d30258447cc213b3b4cd240433 (diff) | |
download | glibc-9b7ee67e0c61b046f86062ef8e2ba91f352ebcc1.tar.gz glibc-9b7ee67e0c61b046f86062ef8e2ba91f352ebcc1.tar.xz glibc-9b7ee67e0c61b046f86062ef8e2ba91f352ebcc1.zip |
Update.
2001-04-16 Stephen L Moshier <moshier@mediaone.net> * sysdeps/ieee754/flt-32/e_asinf.c (pio2_hi, pio2_lo, pio4_hi): Correct the values. (pSx, qSx): Replace by shorter approximation. Use f suffix on float constants. * sysdeps/ieee754/ldbl-128/k_tanl.c: New file. Contributed by Stephen L Moshier <moshier@mediaone.net>.
Diffstat (limited to 'sysdeps/ieee754')
-rw-r--r-- | sysdeps/ieee754/flt-32/e_asinf.c | 51 | ||||
-rw-r--r-- | sysdeps/ieee754/ldbl-128/k_tanl.c | 147 |
2 files changed, 173 insertions, 25 deletions
diff --git a/sysdeps/ieee754/flt-32/e_asinf.c b/sysdeps/ieee754/flt-32/e_asinf.c index 5270f31144..5c3f0d500f 100644 --- a/sysdeps/ieee754/flt-32/e_asinf.c +++ b/sysdeps/ieee754/flt-32/e_asinf.c @@ -13,6 +13,12 @@ * ==================================================== */ +/* + Single precision expansion contributed by + Stephen L. Moshier <moshier@na-net.ornl.gov> +*/ + + #if defined(LIBM_SCCS) && !defined(lint) static char rcsid[] = "$NetBSD: e_asinf.c,v 1.5 1995/05/12 04:57:25 jtc Exp $"; #endif @@ -27,20 +33,19 @@ static float #endif one = 1.0000000000e+00, /* 0x3F800000 */ huge = 1.000e+30, -pio2_hi = 1.5707962513e+00, /* 0x3fc90fda */ -pio2_lo = 7.5497894159e-08, /* 0x33a22168 */ -pio4_hi = 7.8539818525e-01, /* 0x3f490fdb */ - /* coefficient for R(x^2) */ -pS0 = 1.6666667163e-01, /* 0x3e2aaaab */ -pS1 = -3.2556581497e-01, /* 0xbea6b090 */ -pS2 = 2.0121252537e-01, /* 0x3e4e0aa8 */ -pS3 = -4.0055535734e-02, /* 0xbd241146 */ -pS4 = 7.9153501429e-04, /* 0x3a4f7f04 */ -pS5 = 3.4793309169e-05, /* 0x3811ef08 */ -qS1 = -2.4033949375e+00, /* 0xc019d139 */ -qS2 = 2.0209457874e+00, /* 0x4001572d */ -qS3 = -6.8828397989e-01, /* 0xbf303361 */ -qS4 = 7.7038154006e-02; /* 0x3d9dc62e */ + +pio2_hi = 1.57079637050628662109375f, +pio2_lo = -4.37113900018624283e-8f, +pio4_hi = 0.785398185253143310546875f, + +/* asin x = x + x^3 p(x^2) + -0.5 <= x <= 0.5; + Peak relative error 4.8e-9 */ +p0 = 1.666675248e-1f, +p1 = 7.495297643e-2f, +p2 = 4.547037598e-2f, +p3 = 2.417951451e-2f, +p4 = 4.216630880e-2f; #ifdef __STDC__ float __ieee754_asinf(float x) @@ -63,30 +68,26 @@ qS4 = 7.7038154006e-02; /* 0x3d9dc62e */ if(huge+x>one) return x;/* return x with inexact if x!=0*/ } else { t = x*x; - p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5))))); - q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4))); - w = p/q; + w = t * (p0 + t * (p1 + t * (p2 + t * (p3 + t * p4)))); return x+x*w; } } /* 1> |x|>= 0.5 */ w = one-fabsf(x); - t = w*(float)0.5; - p = t*(pS0+t*(pS1+t*(pS2+t*(pS3+t*(pS4+t*pS5))))); - q = one+t*(qS1+t*(qS2+t*(qS3+t*qS4))); + t = w*0.5f; + p = t * (p0 + t * (p1 + t * (p2 + t * (p3 + t * p4)))); s = __ieee754_sqrtf(t); if(ix>=0x3F79999A) { /* if |x| > 0.975 */ - w = p/q; - t = pio2_hi-((float)2.0*(s+s*w)-pio2_lo); + t = pio2_hi-(2.0f*(s+s*p)-pio2_lo); } else { int32_t iw; w = s; GET_FLOAT_WORD(iw,w); SET_FLOAT_WORD(w,iw&0xfffff000); c = (t-w*w)/(s+w); - r = p/q; - p = (float)2.0*s*r-(pio2_lo-(float)2.0*c); - q = pio4_hi-(float)2.0*w; + r = p; + p = 2.0f*s*r-(pio2_lo-2.0f*c); + q = pio4_hi-2.0f*w; t = pio4_hi-(p-q); } if(hx>0) return t; else return -t; diff --git a/sysdeps/ieee754/ldbl-128/k_tanl.c b/sysdeps/ieee754/ldbl-128/k_tanl.c new file mode 100644 index 0000000000..3d8e0342e8 --- /dev/null +++ b/sysdeps/ieee754/ldbl-128/k_tanl.c @@ -0,0 +1,147 @@ +/* + * ==================================================== + * Copyright (C) 1993 by Sun Microsystems, Inc. All rights reserved. + * + * Developed at SunPro, a Sun Microsystems, Inc. business. + * Permission to use, copy, modify, and distribute this + * software is freely granted, provided that this notice + * is preserved. + * ==================================================== + */ + +/* + Long double expansions contributed by + Stephen L. Moshier <moshier@na-net.ornl.gov> +*/ + +/* __kernel_tanl( x, y, k ) + * kernel tan function on [-pi/4, pi/4], pi/4 ~ 0.7854 + * Input x is assumed to be bounded by ~pi/4 in magnitude. + * Input y is the tail of x. + * Input k indicates whether tan (if k=1) or + * -1/tan (if k= -1) is returned. + * + * Algorithm + * 1. Since tan(-x) = -tan(x), we need only to consider positive x. + * 2. if x < 2^-57, return x with inexact if x!=0. + * 3. tan(x) is approximated by a rational form x + x^3 / 3 + x^5 R(x^2) + * on [0,0.67433]. + * + * Note: tan(x+y) = tan(x) + tan'(x)*y + * ~ tan(x) + (1+x*x)*y + * Therefore, for better accuracy in computing tan(x+y), let + * r = x^3 * R(x^2) + * then + * tan(x+y) = x + (x^3 / 3 + (x^2 *(r+y)+y)) + * + * 4. For x in [0.67433,pi/4], let y = pi/4 - x, then + * tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y)) + * = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y))) + */ + +#include "math.h" +#include "math_private.h" +#ifdef __STDC__ +static const long double +#else +static long double +#endif + one = 1.0L, + pio4hi = 7.8539816339744830961566084581987569936977E-1L, + pio4lo = 2.1679525325309452561992610065108379921906E-35L, + + /* tan x = x + x^3 / 3 + x^5 T(x^2)/U(x^2) + 0 <= x <= 0.6743316650390625 + Peak relative error 8.0e-36 */ + TH = 3.333333333333333333333333333333333333333E-1L, + T0 = -1.813014711743583437742363284336855889393E7L, + T1 = 1.320767960008972224312740075083259247618E6L, + T2 = -2.626775478255838182468651821863299023956E4L, + T3 = 1.764573356488504935415411383687150199315E2L, + T4 = -3.333267763822178690794678978979803526092E-1L, + + U0 = -1.359761033807687578306772463253710042010E8L, + U1 = 6.494370630656893175666729313065113194784E7L, + U2 = -4.180787672237927475505536849168729386782E6L, + U3 = 8.031643765106170040139966622980914621521E4L, + U4 = -5.323131271912475695157127875560667378597E2L; + /* 1.000000000000000000000000000000000000000E0 */ + + +#ifdef __STDC__ +long double +__kernel_tanl (long double x, long double y, int iy) +#else +long double +__kernel_tanl (x, y, iy) + long double x, y; + int iy; +#endif +{ + long double z, r, v, w, s; + int32_t ix, sign; + ieee854_long_double_shape_type u, u1; + + u.value = x; + ix = u.parts32.w0 & 0x7fffffff; + if (ix < 0x3fc60000) /* x < 2**-57 */ + { + if ((int) x == 0) + { /* generate inexact */ + if ((ix | u.parts32.w1 | u.parts32.w2 | u.parts32.w3 + | (iy + 1)) == 0) + return one / fabs (x); + else + return (iy == 1) ? x : -one / x; + } + } + if (ix >= 0x3ffe5942) /* |x| >= 0.6743316650390625 */ + { + if ((u.parts32.w0 & 0x80000000) != 0) + { + x = -x; + y = -y; + sign = -1; + } + else + sign = 1; + z = pio4hi - x; + w = pio4lo - y; + x = z + w; + y = 0.0; + } + z = x * x; + r = T0 + z * (T1 + z * (T2 + z * (T3 + z * T4))); + v = U0 + z * (U1 + z * (U2 + z * (U3 + z * (U4 + z)))); + r = r / v; + + s = z * x; + r = y + z * (s * r + y); + r += TH * s; + w = x + r; + if (ix >= 0x3ffe5942) + { + v = (long double) iy; + w = (v - 2.0 * (x - (w * w / (w + v) - r))); + if (sign < 0) + w = -w; + return w; + } + if (iy == 1) + return w; + else + { /* if allow error up to 2 ulp, + simply return -1.0/(x+r) here */ + /* compute -1.0/(x+r) accurately */ + u1.value = w; + u1.parts32.w2 = 0; + u1.parts32.w3 = 0; + v = r - (u1.value - x); /* u1+v = r+x */ + z = -1.0 / w; + u.value = z; + u.parts32.w2 = 0; + u.parts32.w3 = 0; + s = 1.0 + u.value * u1.value; + return u.value + z * (s + u.value * v); + } +} |