aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
authorWilco Dijkstra <wdijkstr@arm.com>2018-04-03 16:49:33 +0100
committerWilco Dijkstra <wdijkstr@arm.com>2018-04-03 16:52:18 +0100
commite88ecbbfe836db5b6da809108759de4ca56be5e7 (patch)
tree0302cf9e5ae61fc04225f67e81423df3f869aaa7
parentaef3e2558a0ab0aff6d80f3e99ebe228321ab4b3 (diff)
downloadglibc-e88ecbbfe836db5b6da809108759de4ca56be5e7.tar
glibc-e88ecbbfe836db5b6da809108759de4ca56be5e7.tar.gz
glibc-e88ecbbfe836db5b6da809108759de4ca56be5e7.tar.bz2
glibc-e88ecbbfe836db5b6da809108759de4ca56be5e7.zip
[PATCH 7/7] sin/cos slow paths: refactor sincos implementation
Refactor the sincos implementation - rather than rely on odd partial inlining of preprocessed portions from sin and cos, explicitly write out the cases. This makes sincos much easier to maintain and provides an additional 16-20% speedup between 0 and 2^27. The overall speedup of sincos is 48% over this range. Between 0 and PI it is 66% faster. * sysdeps/ieee754/dbl-64/s_sin.c (__sin): Cleanup ifdefs. (__cos): Likewise. * sysdeps/ieee754/dbl-64/s_sin.c (__sincos): Refactor using the same logic as sin and cos.
-rw-r--r--ChangeLog7
-rw-r--r--sysdeps/ieee754/dbl-64/s_sin.c29
-rw-r--r--sysdeps/ieee754/dbl-64/s_sincos.c68
3 files changed, 59 insertions, 45 deletions
diff --git a/ChangeLog b/ChangeLog
index 2e5b570c26..36b022cb35 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,12 @@
2018-04-03 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): Refactor using the same
+ logic as sin and cos.
+
+2018-04-03 Wilco Dijkstra <wdijkstr@arm.com>
+
* sysdeps/ieee754/dbl-64/s_sin.c (do_sin): Use TAYLOR_SIN for small
inputs. Return correct sign.
(do_sincos): Remove small input check before do_sin, let do_sin set
diff --git a/sysdeps/ieee754/dbl-64/s_sin.c b/sysdeps/ieee754/dbl-64/s_sin.c
index 6e261f24bb..ba1dbe27b6 100644
--- a/sysdeps/ieee754/dbl-64/s_sin.c
+++ b/sysdeps/ieee754/dbl-64/s_sin.c
@@ -197,27 +197,17 @@ 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)
{
-#ifndef IN_SINCOS
double t, a, da;
mynumber u;
int4 k, m, n;
double retval = 0;
SET_RESTORE_ROUND_53BIT (FE_TONEAREST);
-#else
- double xx, t, cor;
- mynumber u;
- int4 k, m;
- double retval = 0;
-#endif
u.x = x;
m = u.i[HIGH_HALF];
@@ -242,7 +232,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)
{
@@ -263,7 +252,6 @@ __sin (double x)
__set_errno (EDOM);
retval = x / x;
}
-#endif
return retval;
}
@@ -274,27 +262,17 @@ __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, a, da;
mynumber u;
-#ifndef IN_SINCOS
int4 k, m, n;
-#else
- int4 k, m;
-#endif
double retval = 0;
-#ifndef IN_SINCOS
SET_RESTORE_ROUND_53BIT (FE_TONEAREST);
-#endif
u.x = x;
m = u.i[HIGH_HALF];
@@ -320,8 +298,6 @@ __cos (double x)
retval = do_sin (a, da);
} /* else if (k < 0x400368fd) */
-
-#ifndef IN_SINCOS
else if (k < 0x419921FB)
{ /* 2.426265<|x|< 105414350 */
n = reduce_sincos (x, &a, &da);
@@ -341,7 +317,6 @@ __cos (double x)
__set_errno (EDOM);
retval = x / x; /* |x| > 2^1024 */
}
-#endif
return retval;
}
@@ -352,3 +327,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 4335ecbba3..c7460371e4 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,63 @@ __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 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)
+ {
+ *sinx = do_sin (x, 0);
+ *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);
+ *cosx = do_sin (a, da);
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;
+ }
+
+ *sinx = do_sin (a, da);
+ xx = do_cos (a, da);
+ *cosx = (n & 2) ? -xx : xx;
+ return;
}
if (isinf (x))