X-Git-Url: http://erislabs.net/gitweb/?a=blobdiff_plain;f=lib%2Ftanl.c;h=07fbc09451a17ae4d4239b0d021083d63884b88b;hb=9d3af100c7a626ed2b501fb7d0a97936f72ce9ac;hp=e4b974e3d2523fcdfc2b7b0074e4efefd084a5a7;hpb=5657b2536d242966c41c268b9ab4fe8df727bb13;p=gnulib.git diff --git a/lib/tanl.c b/lib/tanl.c index e4b974e3d..07fbc0945 100644 --- a/lib/tanl.c +++ b/lib/tanl.c @@ -14,25 +14,30 @@ * ==================================================== */ +#include + +/* Specification. */ +#include + /* tanl(x) * Return tangent function of x. * * kernel function: - * __kernel_tanl ... tangent function on [-pi/4,pi/4] - * __ieee754_rem_pio2l ... argument reduction routine + * __kernel_tanl ... tangent function on [-pi/4,pi/4] + * __ieee754_rem_pio2l ... argument reduction routine * * Method. * Let S,C and T denote the sin, cos and tan respectively on - * [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2 - * in [-pi/4 , +pi/4], and let n = k mod 4. - * We have + * [-PI/4, +PI/4]. Reduce the argument x to y1+y2 = x-k*pi/2 + * in [-pi/4 , +pi/4], and let n = k mod 4. + * We have * * n sin(x) cos(x) tan(x) * ---------------------------------------------------------- - * 0 S C T - * 1 C -S -1/T - * 2 -S -C T - * 3 -C S -1/T + * 0 S C T + * 1 C -S -1/T + * 2 -S -C T + * 3 -C S -1/T * ---------------------------------------------------------- * * Special cases: @@ -41,19 +46,10 @@ * trig(NaN) is that NaN; * * Accuracy: - * TRIG(x) returns trig(x) nearly rounded + * TRIG(x) returns trig(x) nearly rounded */ -#include - -#include "mathl.h" - #include "trigl.h" -#ifdef HAVE_SINL -#ifdef HAVE_COSL -#include "trigl.c" -#endif -#endif /* * ==================================================== @@ -79,21 +75,21 @@ * -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) + * 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)) + * 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))) + * tan(x) = tan(pi/4-y) = (1-tan(y))/(1+tan(y)) + * = 1 - 2*(tan(y) - (tan(y)^2)/(1+tan(y))) */ @@ -119,33 +115,33 @@ static const long double /* 1.000000000000000000000000000000000000000E0 */ -long double +static long double kernel_tanl (long double x, long double y, int iy) { long double z, r, v, w, s, u, u1; - int flag, sign; + int invert = 0, sign; - sign = 1; - if (x < 0) - { - x = -x; - y = -y; - sign = -1; - } + sign = 1; + if (x < 0) + { + x = -x; + y = -y; + sign = -1; + } if (x < 0.000000000000000006938893903907228377647697925567626953125L) /* x < 2**-57 */ { if ((int) x == 0) - { /* generate inexact */ - if (iy == -1 && x == 0.0) - return 1.0L / fabs (x); - else - return (iy == 1) ? x : -1.0L / x; - } + { /* generate inexact */ + if (iy == -1 && x == 0.0) + return 1.0L / fabs (x); + else + return (iy == 1) ? x : -1.0L / x; + } } if (x >= 0.6743316650390625) /* |x| >= 0.6743316650390625 */ { - flag = 1; + invert = 1; z = pio4hi - x; w = pio4lo - y; @@ -161,19 +157,19 @@ kernel_tanl (long double x, long double y, int iy) r = y + z * (s * r + y); r += TH * s; w = x + r; - if (flag) + if (invert) { v = (long double) iy; w = (v - 2.0 * (x - (w * w / (w + v) - r))); if (sign < 0) - w = -w; + w = -w; return w; } if (iy == 1) return w; else - { /* if allow error up to 2 ulp, - simply return -1.0/(x+r) here */ + { /* if allow error up to 2 ulp, + simply return -1.0/(x+r) here */ /* compute -1.0/(x+r) accurately */ u1 = (double) w; v = r - (u1 - x); @@ -190,14 +186,18 @@ tanl (long double x) long double y[2], z = 0.0L; int n; + /* tanl(NaN) is NaN */ + if (isnanl (x)) + return x; + /* |x| ~< pi/4 */ if (x >= -0.7853981633974483096156608458198757210492 && x <= 0.7853981633974483096156608458198757210492) return kernel_tanl (x, z, 1); - /* tanl(Inf or NaN) is NaN, tanl(0) is 0 */ - else if (x + x == x || x != x) - return x - x; /* NaN */ + /* tanl(Inf) is NaN, tanl(0) is 0 */ + else if (x + x == x) + return x - x; /* NaN */ /* argument reduction needed */ else @@ -210,14 +210,14 @@ tanl (long double x) #if 0 int -main () +main (void) { - printf ("%.16Lg\n", tanl(0.7853981633974483096156608458198757210492)); - printf ("%.16Lg\n", tanl(-0.7853981633974483096156608458198757210492)); - printf ("%.16Lg\n", tanl(0.7853981633974483096156608458198757210492 *3)); - printf ("%.16Lg\n", tanl(-0.7853981633974483096156608458198757210492 *31)); - printf ("%.16Lg\n", tanl(0.7853981633974483096156608458198757210492 / 2)); - printf ("%.16Lg\n", tanl(0.7853981633974483096156608458198757210492 * 3/2)); - printf ("%.16Lg\n", tanl(0.7853981633974483096156608458198757210492 * 5/2)); + printf ("%.16Lg\n", tanl (0.7853981633974483096156608458198757210492)); + printf ("%.16Lg\n", tanl (-0.7853981633974483096156608458198757210492)); + printf ("%.16Lg\n", tanl (0.7853981633974483096156608458198757210492 *3)); + printf ("%.16Lg\n", tanl (-0.7853981633974483096156608458198757210492 *31)); + printf ("%.16Lg\n", tanl (0.7853981633974483096156608458198757210492 / 2)); + printf ("%.16Lg\n", tanl (0.7853981633974483096156608458198757210492 * 3/2)); + printf ("%.16Lg\n", tanl (0.7853981633974483096156608458198757210492 * 5/2)); } #endif