system/glibc
3
0
Fork 0
Code Issues Projects Releases Wiki Activity

[PATCH 4/7] sin/cos slow paths: remove slow paths from huge range reduction

Browse Source 
For huge inputs use the improved do_sincos function as well.  Now no cases use
the correction factor returned by do_sin, do_cos and TAYLOR_SIN, so remove it.

	* sysdeps/ieee754/dbl-64/s_sin.c (TAYLOR_SIN): Remove cor parameter.
	(do_cos): Remove corp parameter and calculations.
	(do_sin): Likewise.
	(do_sincos): Remove cor variable.
	(__sin): Use do_sincos for huge inputs.
	(__cos): Likewise.
	* sysdeps/ieee754/dbl-64/s_sincos.c (__sincos): Likewise.
	(reduce_and_compute_sincos): Remove unused function.
This commit is contained in:
Wilco Dijkstra 2018-04-03 16:41:36 +01:00
parent d9469deb14
commit 649095838b
3 changed files with 45 additions and 64 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

11
ChangeLog
View File

@ -1,3 +1,14 @@
2018-04-03 Wilco Dijkstra <wdijkstr@arm.com>
* sysdeps/ieee754/dbl-64/s_sin.c (TAYLOR_SIN): Remove cor parameter.
(do_cos): Remove corp parameter and calculations.
(do_sin): Likewise.
(do_sincos): Remove cor variable.
(__sin): Use do_sincos for huge inputs.
(__cos): Likewise.
* sysdeps/ieee754/dbl-64/s_sincos.c (__sincos): Likewise.
(reduce_and_compute_sincos): Remove unused function.
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_1): Rename to * sysdeps/ieee754/dbl-64/s_sin.c (reduce_sincos_1): Rename to

60
sysdeps/ieee754/dbl-64/s_sin.c
View File

@ -66,12 +66,11 @@
a - a^3/3! + a^5/5! - a^7/7! + a^9/9! + (1 - a^2) * da / 2 a - a^3/3! + a^5/5! - a^7/7! + a^9/9! + (1 - a^2) * da / 2
The constants s1, s2, s3, etc. are pre-computed values of 1/3!, 1/5! and so The constants s1, s2, s3, etc. are pre-computed values of 1/3!, 1/5! and so
on. The result is returned to LHS and correction in COR. */ on. The result is returned to LHS. */
#define TAYLOR_SIN(xx, a, da, cor) \ #define TAYLOR_SIN(xx, a, da) \
({ \ ({ \
double t = ((POLYNOMIAL (xx) * (a) - 0.5 * (da)) * (xx) + (da)); \ double t = ((POLYNOMIAL (xx) * (a) - 0.5 * (da)) * (xx) + (da)); \
double res = (a) + t; \ double res = (a) + t; \
(cor) = ((a) - res) + t; \
res; \ res; \
}) })
@ -145,10 +144,10 @@ static double cslow2 (double x);
/* Given a number partitioned into X and DX, this function computes the cosine /* Given a number partitioned into X and DX, this function computes the cosine
of the number by combining the sin and cos of X (as computed by a variation of the number by combining the sin and cos of X (as computed by a variation
of the Taylor series) with the values looked up from the sin/cos table to of the Taylor series) with the values looked up from the sin/cos table to
get the result in RES and a correction value in COR. */ get the result. */
static inline double static inline double
__always_inline __always_inline
do_cos (double x, double dx, double *corp) do_cos (double x, double dx)
{ {
mynumber u; mynumber u;
@ -158,16 +157,13 @@ do_cos (double x, double dx, double *corp)
u.x = big + fabs (x); u.x = big + fabs (x);
x = fabs (x) - (u.x - big) + dx; x = fabs (x) - (u.x - big) + dx;
double xx, s, sn, ssn, c, cs, ccs, res, cor; double xx, s, sn, ssn, c, cs, ccs, cor;
xx = x * x; xx = x * x;
s = x + x * xx * (sn3 + xx * sn5); s = x + x * xx * (sn3 + xx * sn5);
c = xx * (cs2 + xx * (cs4 + xx * cs6)); c = xx * (cs2 + xx * (cs4 + xx * cs6));
SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs); SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
cor = (ccs - s * ssn - cs * c) - sn * s; cor = (ccs - s * ssn - cs * c) - sn * s;
res = cs + cor; return cs + cor;
cor = (cs - res) + cor;
*corp = cor;
return res;
} }
/* A more precise variant of DO_COS. EPS is the adjustment to the correction /* A more precise variant of DO_COS. EPS is the adjustment to the correction
@ -207,10 +203,10 @@ do_cos_slow (double x, double dx, double eps, double *corp)
/* Given a number partitioned into X and DX, this function computes the sine of /* Given a number partitioned into X and DX, this function computes the sine of
the number by combining the sin and cos of X (as computed by a variation of the number by combining the sin and cos of X (as computed by a variation of
the Taylor series) with the values looked up from the sin/cos table to get the Taylor series) with the values looked up from the sin/cos table to get
the result in RES and a correction value in COR. */ the result. */
static inline double static inline double
__always_inline __always_inline
do_sin (double x, double dx, double *corp) do_sin (double x, double dx)
{ {
mynumber u; mynumber u;
@ -219,16 +215,13 @@ do_sin (double x, double dx, double *corp)
u.x = big + fabs (x); u.x = big + fabs (x);
x = fabs (x) - (u.x - big); x = fabs (x) - (u.x - big);
double xx, s, sn, ssn, c, cs, ccs, cor, res; double xx, s, sn, ssn, c, cs, ccs, cor;
xx = x * x; xx = x * x;
s = x + (dx + x * xx * (sn3 + xx * sn5)); s = x + (dx + x * xx * (sn3 + xx * sn5));
c = x * dx + xx * (cs2 + xx * (cs4 + xx * cs6)); c = x * dx + xx * (cs2 + xx * (cs4 + xx * cs6));
SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs); SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
cor = (ssn + s * ccs - sn * c) + cs * s; cor = (ssn + s * ccs - sn * c) + cs * s;
res = sn + cor; return sn + cor;
cor = (sn - res) + cor;
*corp = cor;
return res;
} }
/* A more precise variant of DO_SIN. EPS is the adjustment to the correction /* A more precise variant of DO_SIN. EPS is the adjustment to the correction
@ -330,19 +323,19 @@ static double
__always_inline __always_inline
do_sincos (double a, double da, int4 n) do_sincos (double a, double da, int4 n)
{ {
double retval, cor; double retval;
if (n & 1) if (n & 1)
/* Max ULP is 0.513. */ /* Max ULP is 0.513. */
retval = do_cos (a, da, &cor); retval = do_cos (a, da);
else else
{ {
double xx = a * a; double xx = a * a;
/* Max ULP is 0.501 if xx < 0.01588, otherwise ULP is 0.518. */ /* Max ULP is 0.501 if xx < 0.01588, otherwise ULP is 0.518. */
if (xx < 0.01588) if (xx < 0.01588)
retval = TAYLOR_SIN (xx, a, da, cor); retval = TAYLOR_SIN (xx, a, da);
else else
retval = __copysign (do_sin (a, da, &cor), a); retval = __copysign (do_sin (a, da), a);
} }
return (n & 2) ? -retval : retval; return (n & 2) ? -retval : retval;
@ -362,7 +355,7 @@ SECTION
__sin (double x) __sin (double x)
{ {
#ifndef IN_SINCOS #ifndef IN_SINCOS
double xx, t, a, da, cor; double xx, t, a, da;
mynumber u; mynumber u;
int4 k, m, n; int4 k, m, n;
double retval = 0; double retval = 0;
@ -396,7 +389,7 @@ __sin (double x)
else if (k < 0x3feb6000) else if (k < 0x3feb6000)
{ {
/* Max ULP is 0.548. */ /* Max ULP is 0.548. */
retval = __copysign (do_sin (x, 0, &cor), x); retval = __copysign (do_sin (x, 0), x);
} /* else if (k < 0x3feb6000) */ } /* else if (k < 0x3feb6000) */
/*----------------------- 0.855469 <|x|<2.426265 ----------------------*/ /*----------------------- 0.855469 <|x|<2.426265 ----------------------*/
@ -404,7 +397,7 @@ __sin (double x)
{ {
t = hp0 - fabs (x); t = hp0 - fabs (x);
/* Max ULP is 0.51. */ /* Max ULP is 0.51. */
retval = __copysign (do_cos (t, hp1, &cor), x); retval = __copysign (do_cos (t, hp1), x);
} /* else if (k < 0x400368fd) */ } /* else if (k < 0x400368fd) */
#ifndef IN_SINCOS #ifndef IN_SINCOS
@ -417,8 +410,10 @@ __sin (double x)
/* --------------------105414350 <|x| <2^1024------------------------------*/ /* --------------------105414350 <|x| <2^1024------------------------------*/
else if (k < 0x7ff00000) else if (k < 0x7ff00000)
retval = reduce_and_compute (x, false); {
n = __branred (x, &a, &da);
retval = do_sincos (a, da, n);
}
/*--------------------- |x| > 2^1024 ----------------------------------*/ /*--------------------- |x| > 2^1024 ----------------------------------*/
else else
{ {
@ -445,7 +440,7 @@ SECTION
#endif #endif
__cos (double x) __cos (double x)
{ {
double y, xx, cor, a, da; double y, xx, a, da;
mynumber u; mynumber u;
#ifndef IN_SINCOS #ifndef IN_SINCOS
int4 k, m, n; int4 k, m, n;
@ -470,7 +465,7 @@ __cos (double x)
else if (k < 0x3feb6000) else if (k < 0x3feb6000)
{ /* 2^-27 < |x| < 0.855469 */ { /* 2^-27 < |x| < 0.855469 */
/* Max ULP is 0.51. */ /* Max ULP is 0.51. */
retval = do_cos (x, 0, &cor); retval = do_cos (x, 0);
} /* else if (k < 0x3feb6000) */ } /* else if (k < 0x3feb6000) */
else if (k < 0x400368fd) else if (k < 0x400368fd)
@ -482,9 +477,9 @@ __cos (double x)
/* Max ULP is 0.501 if xx < 0.01588 or 0.518 otherwise. /* Max ULP is 0.501 if xx < 0.01588 or 0.518 otherwise.
Range reduction uses 106 bits here which is sufficient. */ Range reduction uses 106 bits here which is sufficient. */
if (xx < 0.01588) if (xx < 0.01588)
retval = TAYLOR_SIN (xx, a, da, cor); retval = TAYLOR_SIN (xx, a, da);
else else
retval = __copysign (do_sin (a, da, &cor), a); retval = __copysign (do_sin (a, da), a);
} /* else if (k < 0x400368fd) */ } /* else if (k < 0x400368fd) */
@ -497,7 +492,10 @@ __cos (double x)
/* 105414350 <|x| <2^1024 */ /* 105414350 <|x| <2^1024 */
else if (k < 0x7ff00000) else if (k < 0x7ff00000)
retval = reduce_and_compute (x, true); {
n = __branred (x, &a, &da);
retval = do_sincos (a, da, n + 1);
}
else else
{ {

38
sysdeps/ieee754/dbl-64/s_sincos.c
View File

@ -28,37 +28,6 @@
#define IN_SINCOS 1 #define IN_SINCOS 1
#include "s_sin.c" #include "s_sin.c"
/* Consolidated version of reduce_and_compute in s_sin.c that does range
reduction only once and computes sin and cos together. */
static inline void
__always_inline
reduce_and_compute_sincos (double x, double *sinx, double *cosx)
{
double a, da;
unsigned int n = __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;
}
if (a * a < 0.01588)
*sinx = bsloww (a, da, x, n);
else
*sinx = bsloww1 (a, da, x, n);
*cosx = bsloww2 (a, da, x, n);
}
void void
__sincos (double x, double *sinx, double *cosx) __sincos (double x, double *sinx, double *cosx)
{ {
@ -88,8 +57,11 @@ __sincos (double x, double *sinx, double *cosx)
} }
if (k < 0x7ff00000) if (k < 0x7ff00000)
{ {
reduce_and_compute_sincos (x, sinx, cosx); double a, da;
return; int4 n = __branred (x, &a, &da);
*sinx = do_sincos (a, da, n);
*cosx = do_sincos (a, da, n + 1);
} }
if (isinf (x)) if (isinf (x))