This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
[PATCH 6/6] Remove slow paths from sin/cos
- From: Wilco Dijkstra <Wilco dot Dijkstra at arm dot com>
- To: "libc-alpha at sourceware dot org" <libc-alpha at sourceware dot org>
- Cc: nd <nd at arm dot com>
- Date: Fri, 9 Mar 2018 15:52:55 +0000
- Subject: [PATCH 6/6] Remove slow paths from sin/cos
- Authentication-results: sourceware.org; auth=none
- Authentication-results: spf=none (sender IP is ) smtp.mailfrom=Wilco dot Dijkstra at arm dot com;
- Nodisclaimer: True
- Spamdiagnosticmetadata: NSPM
- Spamdiagnosticoutput: 1:99
Restructure the sincos implementation - rather than rely on odd partial inlining
of preprocessed portions from sin and cos, explicitly write out the cases.
ChangeLog:
2018-03-09 Wilco Dijkstra <wdijkstr@arm.com>
* sysdeps/ieee754/dbl-64/s_sin.c (__sin): Cleanup ifdefs.
(__cos): Likewise.
* sysdeps/ieee754/dbl-64/s_sin.c (__sincos): Reimplement using the same
logic as sin and cos.
--
diff --git a/sysdeps/ieee754/dbl-64/s_sin.c b/sysdeps/ieee754/dbl-64/s_sin.c
index 91b0abc9c3ac21dae0e673576940ef97bfd20c23..8f804a42e6d94652a62f81b2e0b053135cf9f03a 100644
--- a/sysdeps/ieee754/dbl-64/s_sin.c
+++ b/sysdeps/ieee754/dbl-64/s_sin.c
@@ -208,12 +208,9 @@ do_sincos (double a, double da, int4 n)
/* An ultimate sin routine. Given an IEEE double machine number x */
/* it computes the correctly rounded (to nearest) value of sin(x) */
/*******************************************************************/
-#ifdef IN_SINCOS
-static double
-#else
+#ifndef IN_SINCOS
double
SECTION
-#endif
__sin (double x)
{
double xx, t, a, da;
@@ -221,9 +218,7 @@ __sin (double x)
int4 k, m, n;
double retval = 0;
-#ifndef IN_SINCOS
SET_RESTORE_ROUND_53BIT (FE_TONEAREST);
-#endif
u.x = x;
m = u.i[HIGH_HALF];
@@ -257,7 +252,6 @@ __sin (double x)
retval = __copysign (do_cos (t, hp1), x);
} /* else if (k < 0x400368fd) */
-#ifndef IN_SINCOS
/*-------------------------- 2.426265<|x|< 105414350 ----------------------*/
else if (k < 0x419921FB)
{
@@ -278,12 +272,6 @@ __sin (double x)
__set_errno (EDOM);
retval = x / x;
}
-#else
- /* Disable warning... */
- n = 0, n = n;
- a = 0, a = a;
- da = 0, da = da;
-#endif
return retval;
}
@@ -294,12 +282,8 @@ __sin (double x)
/* it computes the correctly rounded (to nearest) value of cos(x) */
/*******************************************************************/
-#ifdef IN_SINCOS
-static double
-#else
double
SECTION
-#endif
__cos (double x)
{
double y, xx, a, da;
@@ -308,9 +292,7 @@ __cos (double x)
double retval = 0;
-#ifndef IN_SINCOS
SET_RESTORE_ROUND_53BIT (FE_TONEAREST);
-#endif
u.x = x;
m = u.i[HIGH_HALF];
@@ -340,8 +322,6 @@ __cos (double x)
retval = __copysign (do_sin (a, da), a);
} /* else if (k < 0x400368fd) */
-
-#ifndef IN_SINCOS
else if (k < 0x419921FB)
{ /* 2.426265<|x|< 105414350 */
n = reduce_sincos (x, &a, &da);
@@ -361,10 +341,6 @@ __cos (double x)
__set_errno (EDOM);
retval = x / x; /* |x| > 2^1024 */
}
-#else
- /* Disable warning... */
- n = 0, n = n;
-#endif
return retval;
}
@@ -375,3 +351,5 @@ libm_alias_double (__cos, cos)
#ifndef __sin
libm_alias_double (__sin, sin)
#endif
+
+#endif
diff --git a/sysdeps/ieee754/dbl-64/s_sincos.c b/sysdeps/ieee754/dbl-64/s_sincos.c
index 4335ecbba3c9894e61c087ac970b392fa73abfab..c04972707b284e37b15e82933a00250cda959985 100644
--- a/sysdeps/ieee754/dbl-64/s_sincos.c
+++ b/sysdeps/ieee754/dbl-64/s_sincos.c
@@ -23,9 +23,7 @@
#include <math_private.h>
#include <libm-alias-double.h>
-#define __sin __sin_local
-#define __cos __cos_local
-#define IN_SINCOS 1
+#define IN_SINCOS
#include "s_sin.c"
void
@@ -37,31 +35,79 @@ __sincos (double x, double *sinx, double *cosx)
SET_RESTORE_ROUND_53BIT (FE_TONEAREST);
u.x = x;
- k = 0x7fffffff & u.i[HIGH_HALF];
+ k = u.i[HIGH_HALF] & 0x7fffffff;
if (k < 0x400368fd)
{
- *sinx = __sin_local (x);
- *cosx = __cos_local (x);
- return;
- }
- if (k < 0x419921FB)
- {
- double a, da;
- int4 n = reduce_sincos (x, &a, &da);
-
- *sinx = do_sincos (a, da, n);
- *cosx = do_sincos (a, da, n + 1);
+ double t, xx, a, da, y;
+ /* |x| < 2^-27 => cos (x) = 1, sin (x) = x. */
+ if (k < 0x3e400000)
+ {
+ if (k < 0x3e500000)
+ math_check_force_underflow (x);
+ *sinx = x;
+ *cosx = 1.0;
+ return;
+ }
+ /* |x| < 0.855469. */
+ else if (k < 0x3feb6000)
+ {
+ /* |x| < 0.25. */
+ if (k < 0x3fd00000)
+ {
+ xx = x * x;
+ t = POLYNOMIAL (xx) * (xx * x);
+ *sinx = x + t;
+ }
+ else
+ *sinx = __copysign (do_sin (x, 0), x);
+ *cosx = do_cos (x, 0);
+ return;
+ }
+ /* |x| < 2.426265. */
+ y = hp0 - fabs (x);
+ a = y + hp1;
+ da = (y - a) + hp1;
+ *sinx = __copysign (do_cos (a, da), x);
+ xx = a * a;
+ if (xx < 0.01588)
+ *cosx = TAYLOR_SIN (xx, a, da);
+ else
+ *cosx = __copysign (do_sin (a, da), a);
return;
}
+ /* |x| < 2^1024. */
if (k < 0x7ff00000)
{
- double a, da;
- int4 n = __branred (x, &a, &da);
+ double a, da, xx;
+ unsigned int n;
- *sinx = do_sincos (a, da, n);
- *cosx = do_sincos (a, da, n + 1);
+ /* If |x| < 105414350 use simple range reduction. */
+ n = k < 0x419921FB ? reduce_sincos (x, &a, &da) : __branred (x, &a, &da);
+ n = n & 3;
+
+ if (n == 1 || n == 2)
+ {
+ a = -a;
+ da = -da;
+ }
+
+ if (n & 1)
+ {
+ double *temp = cosx;
+ cosx = sinx;
+ sinx = temp;
+ }
+
+ xx = a * a;
+ if (xx < 0.01588)
+ *sinx = TAYLOR_SIN (xx, a, da);
+ else
+ *sinx = __copysign (do_sin (a, da), a);
+ xx = do_cos (a, da);
+ *cosx = (n & 2) ? -xx : xx;
+ return;
}
if (isinf (x))