[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.
This commit is contained in:
Wilco Dijkstra 2018-04-03 16:49:33 +01:00 committed by Fangrui Song
parent c8aaaf67f6
commit 2d20ffe431
2 changed files with 53 additions and 46 deletions

View File

@ -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

View File

@ -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))