This is the mail archive of the
libc-alpha@sourceware.org
mailing list for the glibc project.
[PATCH 3/5] Consolidate input partitioning into do_cos and do_sin
- From: Siddhesh Poyarekar <siddhesh at sourceware dot org>
- To: libc-alpha at sourceware dot org
- Date: Tue, 23 Aug 2016 23:52:43 +0530
- Subject: [PATCH 3/5] Consolidate input partitioning into do_cos and do_sin
- Authentication-results: sourceware.org; auth=none
- References: <1471976565-3576-1-git-send-email-siddhesh@sourceware.org>
All calls to do_cos are preceded by code that partitions x into a
larger double that gives an offset into the sincos table and a smaller
double that is used in a polynomial computation. Consolidate all of
them into do_cos and do_sin to reduce code duplication.
* sysdeps/ieee754/dbl-64/s_sin.c (do_cos): Accept X and DX as input
arguments. Consolidate input partitioning from callers here.
(do_cos_slow): Likewise.
(do_sin): Likewise.
(do_sin_slow): Likewise.
(do_sincos_1): Remove the no longer necessary input partitioning.
(do_sincos_2): Likewise.
(__sin): Likewise.
(__cos): Likewise.
(slow1): Likewise.
(slow2): Likewise.
(sloww1): Likewise.
(sloww2): Likewise.
(bsloww1): Likewise.
(bsloww2): Likewise.
(cslow2): Likewise.
---
sysdeps/ieee754/dbl-64/s_sin.c | 191 ++++++++++++++++++-----------------------
1 file changed, 82 insertions(+), 109 deletions(-)
diff --git a/sysdeps/ieee754/dbl-64/s_sin.c b/sysdeps/ieee754/dbl-64/s_sin.c
index 7f6cd09..e03c75a 100644
--- a/sysdeps/ieee754/dbl-64/s_sin.c
+++ b/sysdeps/ieee754/dbl-64/s_sin.c
@@ -141,14 +141,21 @@ static double bsloww2 (double x, double dx, double orig, int n);
int __branred (double x, double *a, double *aa);
static double cslow2 (double x);
-/* Given a number partitioned into U and X such that U is an index into the
- sin/cos table, this macro computes the cosine 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 result in RES and a
- correction value in COR. */
+/* 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 Taylor series) with the values looked up from the sin/cos table to
+ get the result in RES and a correction value in COR. */
static double
-do_cos (mynumber u, double x, double *corp)
+do_cos (double x, double dx, double *corp)
{
+ mynumber u;
+
+ if (x < 0)
+ dx = -dx;
+
+ u.x = big + fabs (x);
+ x = fabs (x) - (u.x - big) + dx;
+
double xx, s, sn, ssn, c, cs, ccs, res, cor;
xx = x * x;
s = x + x * xx * (sn3 + xx * sn5);
@@ -161,11 +168,19 @@ do_cos (mynumber u, double x, double *corp)
return res;
}
-/* A more precise variant of DO_COS where the number is partitioned into U, X
- and DX. EPS is the adjustment to the correction COR. */
+/* A more precise variant of DO_COS. EPS is the adjustment to the correction
+ COR. */
static double
-do_cos_slow (mynumber u, double x, double dx, double eps, double *corp)
+do_cos_slow (double x, double dx, double eps, double *corp)
{
+ mynumber u;
+
+ if (x <= 0)
+ dx = -dx;
+
+ u.x = big + fabs (x);
+ x = fabs (x) - (u.x - big);
+
double xx, y, x1, x2, e1, e2, res, cor;
double s, sn, ssn, c, cs, ccs;
xx = x * x;
@@ -186,14 +201,20 @@ do_cos_slow (mynumber u, double x, double dx, double eps, double *corp)
return res;
}
-/* Given a number partitioned into U and X and DX such that U is an index into
- the sin/cos table, this macro computes the sine 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 result in RES and a
- correction value in COR. */
+/* 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 Taylor series) with the values looked up from the sin/cos table to get
+ the result in RES and a correction value in COR. */
static double
-do_sin (mynumber u, double x, double dx, double *corp)
+do_sin (double x, double dx, double *corp)
{
+ mynumber u;
+
+ if (x <= 0)
+ dx = -dx;
+ u.x = big + fabs (x);
+ x = fabs (x) - (u.x - big);
+
double xx, s, sn, ssn, c, cs, ccs, cor, res;
xx = x * x;
s = x + (dx + x * xx * (sn3 + xx * sn5));
@@ -206,11 +227,18 @@ do_sin (mynumber u, double x, double dx, double *corp)
return res;
}
-/* A more precise variant of res = do_sin where the number is partitioned into U, X
- and DX. EPS is the adjustment to the correction COR. */
+/* A more precise variant of DO_SIN. EPS is the adjustment to the correction
+ COR. */
static double
-do_sin_slow (mynumber u, double x, double dx, double eps, double *corp)
+do_sin_slow (double x, double dx, double eps, double *corp)
{
+ mynumber u;
+
+ if (x <= 0)
+ dx = -dx;
+ u.x = big + fabs (x);
+ x = fabs (x) - (u.x - big);
+
double xx, y, x1, x2, c1, c2, res, cor;
double s, sn, ssn, c, cs, ccs;
xx = x * x;
@@ -288,8 +316,7 @@ static double
__always_inline
do_sincos_1 (double a, double da, double x, int4 n, int4 k)
{
- double xx, retval, res, cor, y;
- mynumber u;
+ double xx, retval, res, cor;
double eps = fabs (x) * 1.2e-30;
int k1 = (n + k) & 3;
@@ -309,10 +336,7 @@ do_sincos_1 (double a, double da, double x, int4 n, int4 k)
}
else
{
- double db = (a > 0 ? da : -da);
- u.x = big + fabs (a);
- y = fabs (a) - (u.x - big);
- res = do_sin (u, y, db, &cor);
+ res = do_sin (a, da, &cor);
cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps;
retval = ((res == res + cor) ? ((a > 0) ? res : -res)
: sloww1 (a, da, x, k));
@@ -321,16 +345,11 @@ do_sincos_1 (double a, double da, double x, int4 n, int4 k)
case 1:
case 3:
- {
- double db = (a > 0 ? da : -da);
- u.x = big + fabs (a);
- y = fabs (a) - (u.x - big) + db;
- res = do_cos (u, y, &cor);
- cor = (cor > 0) ? 1.025 * cor + eps : 1.025 * cor - eps;
- retval = ((res == res + cor) ? ((k1 & 2) ? -res : res)
- : sloww2 (a, da, x, n));
- break;
- }
+ res = do_cos (a, da, &cor);
+ cor = (cor > 0) ? 1.025 * cor + eps : 1.025 * cor - eps;
+ retval = ((res == res + cor) ? ((k1 & 2) ? -res : res)
+ : sloww2 (a, da, x, n));
+ break;
}
return retval;
@@ -369,7 +388,6 @@ __always_inline
do_sincos_2 (double a, double da, double x, int4 n, int4 k)
{
double res, retval, cor, xx;
- mynumber u;
double eps = 1.0e-24;
@@ -392,10 +410,7 @@ do_sincos_2 (double a, double da, double x, int4 n, int4 k)
}
else
{
- double db = (a > 0 ? da : -da);
- u.x = big + fabs (a);
- double y = fabs (a) - (u.x - big);
- res = do_sin (u, y, db, &cor);
+ res = do_sin (a, da, &cor);
cor = (cor > 0) ? 1.035 * cor + eps : 1.035 * cor - eps;
retval = ((res == res + cor) ? ((a > 0) ? res : -res)
: bsloww1 (a, da, x, n));
@@ -404,16 +419,11 @@ do_sincos_2 (double a, double da, double x, int4 n, int4 k)
case 1:
case 3:
- {
- double db = (a > 0 ? da : -da);
- u.x = big + fabs (a);
- double y = fabs (a) - (u.x - big) + db;
- res = do_cos (u, y, &cor);
- cor = (cor > 0) ? 1.025 * cor + eps : 1.025 * cor - eps;
- retval = ((res == res + cor) ? ((n & 2) ? -res : res)
- : bsloww2 (a, da, x, n));
- break;
- }
+ res = do_cos (a, da, &cor);
+ cor = (cor > 0) ? 1.025 * cor + eps : 1.025 * cor - eps;
+ retval = ((res == res + cor) ? ((n & 2) ? -res : res)
+ : bsloww2 (a, da, x, n));
+ break;
}
return retval;
@@ -485,11 +495,7 @@ __sin (double x)
{
t = hp0 - fabs (x);
- u.x = big + fabs (t);
- y = fabs (t) - (u.x - big);
- y = ((t >= 0) ? hp1 : -hp1) + y;
-
- res = do_cos (u, y, &cor);
+ res = do_cos (t, hp1, &cor);
retval = (res == res + 1.020 * cor) ? ((m > 0) ? res : -res) : slow2 (x);
} /* else if (k < 0x400368fd) */
@@ -561,10 +567,7 @@ __cos (double x)
else if (k < 0x3feb6000)
{ /* 2^-27 < |x| < 0.855469 */
- y = fabs (x);
- u.x = big + y;
- y = y - (u.x - big);
- res = do_cos (u, y, &cor);
+ res = do_cos (x, 0, &cor);
retval = (res == res + 1.020 * cor) ? res : cslow2 (x);
} /* else if (k < 0x3feb6000) */
@@ -582,10 +585,7 @@ __cos (double x)
}
else
{
- double db = (a > 0 ? da : -da);
- u.x = big + fabs (a);
- y = fabs (a) - (u.x - big);
- res = do_sin (u, y, db, &cor);
+ res = do_sin (a, da, &cor);
cor = (cor > 0) ? 1.035 * cor + 1.0e-31 : 1.035 * cor - 1.0e-31;
retval = ((res == res + cor) ? ((a > 0) ? res : -res)
: sloww1 (a, da, x, 1));
@@ -655,12 +655,9 @@ static double
SECTION
slow1 (double x)
{
- mynumber u;
- double w[2], y, cor, res;
- y = fabs (x);
- u.x = big + y;
- y = y - (u.x - big);
- res = do_sin_slow (u, y, 0, 0, &cor);
+ double w[2], cor, res;
+
+ res = do_sin_slow (x, 0, 0, &cor);
if (res == res + cor)
return (x > 0) ? res : -res;
@@ -679,15 +676,10 @@ static double
SECTION
slow2 (double x)
{
- mynumber u;
- double w[2], y, y1, y2, cor, res, del;
+ double w[2], y, y1, y2, cor, res;
double t = hp0 - fabs (x);
- u.x = big + fabs (t);
- y = fabs (t) - (u.x - big);
- del = (t >= 0) ? hp1 : -hp1;
-
- res = do_cos_slow (u, y, del, 0, &cor);
+ res = do_cos_slow (t, hp1, 0, &cor);
if (res == res + cor)
return (x > 0) ? res : -res;
@@ -774,17 +766,14 @@ static double
SECTION
sloww1 (double x, double dx, double orig, int k)
{
- mynumber u;
- double w[2], y, cor, res;
+ double w[2], cor, res;
- u.x = big + fabs (x);
- y = fabs (x) - (u.x - big);
- dx = (x > 0 ? dx : -dx);
- res = do_sin_slow (u, y, dx, 3.1e-30 * fabs (orig), &cor);
+ res = do_sin_slow (x, dx, 3.1e-30 * fabs (orig), &cor);
if (res == res + cor)
return (x > 0) ? res : -res;
+ dx = (x > 0 ? dx : -dx);
__dubsin (fabs (x), dx, w);
double eps = 1.1e-30 * fabs (orig);
@@ -807,17 +796,14 @@ static double
SECTION
sloww2 (double x, double dx, double orig, int n)
{
- mynumber u;
- double w[2], y, cor, res;
+ double w[2], cor, res;
- u.x = big + fabs (x);
- y = fabs (x) - (u.x - big);
- dx = (x > 0 ? dx : -dx);
- res = do_cos_slow (u, y, dx, 3.1e-30 * fabs (orig), &cor);
+ res = do_cos_slow (x, dx, 3.1e-30 * fabs (orig), &cor);
if (res == res + cor)
return (n & 2) ? -res : res;
+ dx = x > 0 ? dx : -dx;
__docos (fabs (x), dx, w);
double eps = 1.1e-30 * fabs (orig);
@@ -870,17 +856,13 @@ static double
SECTION
bsloww1 (double x, double dx, double orig, int n)
{
- mynumber u;
- double w[2], y, cor, res;
+ double w[2], cor, res;
- y = fabs (x);
- u.x = big + y;
- y = y - (u.x - big);
- dx = (x > 0) ? dx : -dx;
- res = do_sin_slow (u, y, dx, 1.1e-24, &cor);
+ res = do_sin_slow (x, dx, 1.1e-24, &cor);
if (res == res + cor)
return (x > 0) ? res : -res;
+ dx = (x > 0) ? dx : -dx;
__dubsin (fabs (x), dx, w);
cor = 1.000000005 * w[1] + ((w[1] > 0) ? 1.1e-24 : -1.1e-24);
@@ -902,17 +884,13 @@ static double
SECTION
bsloww2 (double x, double dx, double orig, int n)
{
- mynumber u;
- double w[2], y, cor, res;
+ double w[2], cor, res;
- y = fabs (x);
- u.x = big + y;
- y = y - (u.x - big);
- dx = (x > 0) ? dx : -dx;
- res = do_cos_slow (u, y, dx, 1.1e-24, &cor);
+ res = do_cos_slow (x, dx, 1.1e-24, &cor);
if (res == res + cor)
return (n & 2) ? -res : res;
+ dx = (x > 0) ? dx : -dx;
__docos (fabs (x), dx, w);
cor = 1.000000005 * w[1] + ((w[1] > 0) ? 1.1e-24 : -1.1e-24);
@@ -932,18 +910,13 @@ static double
SECTION
cslow2 (double x)
{
- mynumber u;
- double w[2], y, cor, res;
+ double w[2], cor, res;
- y = fabs (x);
- u.x = big + y;
- y = y - (u.x - big);
- res = do_cos_slow (u, y, 0, 0, &cor);
+ res = do_cos_slow (x, 0, 0, &cor);
if (res == res + cor)
return res;
- y = fabs (x);
- __docos (y, 0, w);
+ __docos (fabs (x), 0, w);
if (w[0] == w[0] + 1.000000005 * w[1])
return w[0];
--
2.7.4