[PATCH 3/7] sin/cos slow paths: remove slow paths from small range reduction

This patch improves the accuracy of the range reduction.  When the input is
large (2^27) and very close to a multiple of PI/2, using 110 bits of PI is not
enough.  Improve range reduction accuracy to 136 bits.  As a result the special
checks for results close to zero can be removed.  The ULP of the polynomials is
at worst 0.55ULP, so there is no reason for the slow functions, and they can be
removed.

	* sysdeps/ieee754/dbl-64/s_sin.c (reduce_sincos_1): Rename to
	reduce_sincos, improve accuracy to 136 bits.
	(do_sincos_1): Rename to do_sincos, remove fallbacks to slow functions.
	(__sin): Use improved reduction and simplified do_sincos calculation.
	(__cos): Likewise.
	* sysdeps/ieee754/dbl-64/s_sincos.c (__sincos): Likewise.
This commit is contained in:
Wilco Dijkstra 2018-04-03 16:33:13 +01:00
parent 7a5640f23a
commit d9469deb14
3 changed files with 56 additions and 53 deletions

View File

@ -1,3 +1,12 @@
2018-04-03 Wilco Dijkstra <wdijkstr@arm.com>
* sysdeps/ieee754/dbl-64/s_sin.c (reduce_sincos_1): Rename to
reduce_sincos, improve accuracy to 136 bits.
(do_sincos_1): Rename to do_sincos, remove fallbacks to slow functions.
(__sin): Use improved reduction and simplified do_sincos calculation.
(__cos): Likewise.
* sysdeps/ieee754/dbl-64/s_sincos.c (__sincos): Likewise.
2018-04-03 Wilco Dijkstra <wdijkstr@arm.com> 2018-04-03 Wilco Dijkstra <wdijkstr@arm.com>
* sysdeps/ieee754/dbl-64/s_sin.c (reduce_sincos_2): Remove function. * sysdeps/ieee754/dbl-64/s_sin.c (reduce_sincos_2): Remove function.

View File

@ -295,9 +295,13 @@ reduce_and_compute (double x, bool shift_quadrant)
return retval; return retval;
} }
/* Reduce range of x to within PI/2 with abs (x) < 105414350. The high part
is written to *a, the low part to *da. Range reduction is accurate to 136
bits so that when x is large and *a very close to zero, all 53 bits of *a
are correct. */
static inline int4 static inline int4
__always_inline __always_inline
reduce_sincos_1 (double x, double *a, double *da) reduce_sincos (double x, double *a, double *da)
{ {
mynumber v; mynumber v;
@ -306,62 +310,45 @@ reduce_sincos_1 (double x, double *a, double *da)
v.x = t; v.x = t;
double y = (x - xn * mp1) - xn * mp2; double y = (x - xn * mp1) - xn * mp2;
int4 n = v.i[LOW_HALF] & 3; int4 n = v.i[LOW_HALF] & 3;
double db = xn * mp3;
double b = y - db; double b, db, t1, t2;
db = (y - b) - db; t1 = xn * pp3;
t2 = y - t1;
db = (y - t2) - t1;
t1 = xn * pp4;
b = t2 - t1;
db += (t2 - b) - t1;
*a = b; *a = b;
*da = db; *da = db;
return n; return n;
} }
/* Compute sin (A + DA). cos can be computed by passing SHIFT_QUADRANT as /* Compute sin or cos (A + DA) for the given quadrant N. */
true, which results in shifting the quadrant N clockwise. */
static double static double
__always_inline __always_inline
do_sincos_1 (double a, double da, double x, int4 n, bool shift_quadrant) do_sincos (double a, double da, int4 n)
{ {
double xx, retval, res, cor; double retval, cor;
double eps = fabs (x) * 1.2e-30;
int k1 = (n + shift_quadrant) & 3; if (n & 1)
switch (k1) /* Max ULP is 0.513. */
{ /* quarter of unit circle */ retval = do_cos (a, da, &cor);
case 2:
a = -a;
da = -da;
/* Fall through. */
case 0:
xx = a * a;
if (xx < 0.01588)
{
/* Taylor series. */
res = TAYLOR_SIN (xx, a, da, cor);
cor = 1.02 * cor + __copysign (eps, cor);
retval = (res == res + cor) ? res : sloww (a, da, x, shift_quadrant);
}
else else
{ {
res = do_sin (a, da, &cor); double xx = a * a;
cor = 1.035 * cor + __copysign (eps, cor); /* Max ULP is 0.501 if xx < 0.01588, otherwise ULP is 0.518. */
retval = ((res == res + cor) ? __copysign (res, a) if (xx < 0.01588)
: sloww1 (a, da, x, shift_quadrant)); retval = TAYLOR_SIN (xx, a, da, cor);
} else
break; retval = __copysign (do_sin (a, da, &cor), a);
case 1:
case 3:
res = do_cos (a, da, &cor);
cor = 1.025 * cor + __copysign (eps, cor);
retval = ((res == res + cor) ? ((n & 2) ? -res : res)
: sloww2 (a, da, x, n));
break;
} }
return retval; return (n & 2) ? -retval : retval;
} }
/*******************************************************************/ /*******************************************************************/
/* An ultimate sin routine. Given an IEEE double machine number x */ /* An ultimate sin routine. Given an IEEE double machine number x */
/* it computes the correctly rounded (to nearest) value of sin(x) */ /* it computes the correctly rounded (to nearest) value of sin(x) */
@ -374,13 +361,18 @@ SECTION
#endif #endif
__sin (double x) __sin (double x)
{ {
#ifndef IN_SINCOS
double xx, t, a, da, cor;
mynumber u;
int4 k, m, n;
double retval = 0;
SET_RESTORE_ROUND_53BIT (FE_TONEAREST);
#else
double xx, t, cor; double xx, t, cor;
mynumber u; mynumber u;
int4 k, m; int4 k, m;
double retval = 0; double retval = 0;
#ifndef IN_SINCOS
SET_RESTORE_ROUND_53BIT (FE_TONEAREST);
#endif #endif
u.x = x; u.x = x;
@ -419,9 +411,8 @@ __sin (double x)
/*-------------------------- 2.426265<|x|< 105414350 ----------------------*/ /*-------------------------- 2.426265<|x|< 105414350 ----------------------*/
else if (k < 0x419921FB) else if (k < 0x419921FB)
{ {
double a, da; n = reduce_sincos (x, &a, &da);
int4 n = reduce_sincos_1 (x, &a, &da); retval = do_sincos (a, da, n);
retval = do_sincos_1 (a, da, x, n, false);
} /* else if (k < 0x419921FB ) */ } /* else if (k < 0x419921FB ) */
/* --------------------105414350 <|x| <2^1024------------------------------*/ /* --------------------105414350 <|x| <2^1024------------------------------*/
@ -456,7 +447,11 @@ __cos (double x)
{ {
double y, xx, cor, a, da; double y, xx, cor, a, da;
mynumber u; mynumber u;
#ifndef IN_SINCOS
int4 k, m, n;
#else
int4 k, m; int4 k, m;
#endif
double retval = 0; double retval = 0;
@ -496,9 +491,8 @@ __cos (double x)
#ifndef IN_SINCOS #ifndef IN_SINCOS
else if (k < 0x419921FB) else if (k < 0x419921FB)
{ /* 2.426265<|x|< 105414350 */ { /* 2.426265<|x|< 105414350 */
double a, da; n = reduce_sincos (x, &a, &da);
int4 n = reduce_sincos_1 (x, &a, &da); retval = do_sincos (a, da, n + 1);
retval = do_sincos_1 (a, da, x, n, true);
} /* else if (k < 0x419921FB ) */ } /* else if (k < 0x419921FB ) */
/* 105414350 <|x| <2^1024 */ /* 105414350 <|x| <2^1024 */

View File

@ -79,10 +79,10 @@ __sincos (double x, double *sinx, double *cosx)
if (k < 0x419921FB) if (k < 0x419921FB)
{ {
double a, da; double a, da;
int4 n = reduce_sincos_1 (x, &a, &da); int4 n = reduce_sincos (x, &a, &da);
*sinx = do_sincos_1 (a, da, x, n, false); *sinx = do_sincos (a, da, n);
*cosx = do_sincos_1 (a, da, x, n, true); *cosx = do_sincos (a, da, n + 1);
return; return;
} }