From 09544cbcd6ef9e5ea2553c8b410dd55712171c33 Mon Sep 17 00:00:00 2001 From: Siddhesh Poyarekar Date: Tue, 8 Oct 2013 11:50:17 +0530 Subject: Consolidate multiple precision sin/cos functions --- sysdeps/ieee754/dbl-64/s_sin.c | 32 +++--- sysdeps/ieee754/dbl-64/sincos32.c | 204 ++++++++++++++++++-------------------- 2 files changed, 111 insertions(+), 125 deletions(-) (limited to 'sysdeps/ieee754') diff --git a/sysdeps/ieee754/dbl-64/s_sin.c b/sysdeps/ieee754/dbl-64/s_sin.c index 2be8fe3a76..b3dfef2406 100644 --- a/sysdeps/ieee754/dbl-64/s_sin.c +++ b/sysdeps/ieee754/dbl-64/s_sin.c @@ -127,10 +127,8 @@ static const double void __dubsin (double x, double dx, double w[]); void __docos (double x, double dx, double w[]); -double __mpsin (double x, double dx); -double __mpcos (double x, double dx); -double __mpsin1 (double x); -double __mpcos1 (double x); +double __mpsin (double x, double dx, bool reduce_range); +double __mpcos (double x, double dx, bool reduce_range); static double slow (double x); static double slow1 (double x); static double slow2 (double x); @@ -722,7 +720,7 @@ slow (double x) if (w[0] == w[0] + 1.000000001 * w[1]) return (x > 0) ? w[0] : -w[0]; else - return (x > 0) ? __mpsin (x, 0) : -__mpsin (-x, 0); + return (x > 0) ? __mpsin (x, 0, false) : -__mpsin (-x, 0, false); } } @@ -762,7 +760,7 @@ slow1 (double x) if (w[0] == w[0] + 1.000000005 * w[1]) return (x > 0) ? w[0] : -w[0]; else - return (x > 0) ? __mpsin (x, 0) : -__mpsin (-x, 0); + return (x > 0) ? __mpsin (x, 0, false) : -__mpsin (-x, 0, false); } } @@ -815,7 +813,7 @@ slow2 (double x) if (w[0] == w[0] + 1.000000005 * w[1]) return (x > 0) ? w[0] : -w[0]; else - return (x > 0) ? __mpsin (x, 0) : -__mpsin (-x, 0); + return (x > 0) ? __mpsin (x, 0, false) : -__mpsin (-x, 0, false); } } @@ -882,7 +880,7 @@ sloww (double x, double dx, double orig) if (w[0] == w[0] + cor) return (a > 0) ? w[0] : -w[0]; else - return __mpsin1 (orig); + return __mpsin (orig, 0, true); } } } @@ -939,7 +937,7 @@ sloww1 (double x, double dx, double orig) if (w[0] == w[0] + cor) return (x > 0) ? w[0] : -w[0]; else - return __mpsin1 (orig); + return __mpsin (orig, 0, true); } } @@ -996,7 +994,7 @@ sloww2 (double x, double dx, double orig, int n) if (w[0] == w[0] + cor) return (n & 2) ? -w[0] : w[0]; else - return __mpsin1 (orig); + return __mpsin (orig, 0, true); } } @@ -1028,7 +1026,7 @@ bsloww (double x, double dx, double orig, int n) if (w[0] == w[0] + cor) return (x > 0) ? w[0] : -w[0]; else - return (n & 1) ? __mpcos1 (orig) : __mpsin1 (orig); + return (n & 1) ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true); } } @@ -1079,7 +1077,7 @@ bsloww1 (double x, double dx, double orig, int n) if (w[0] == w[0] + cor) return (x > 0) ? w[0] : -w[0]; else - return (n & 1) ? __mpcos1 (orig) : __mpsin1 (orig); + return (n & 1) ? __mpcos (orig, 0, true) : __mpsin (orig, 0, true); } } @@ -1131,7 +1129,7 @@ bsloww2 (double x, double dx, double orig, int n) if (w[0] == w[0] + cor) return (n & 2) ? -w[0] : w[0]; else - return (n & 1) ? __mpsin1 (orig) : __mpcos1 (orig); + return (n & 1) ? __mpsin (orig, 0, true) : __mpcos (orig, 0, true); } } @@ -1173,7 +1171,7 @@ cslow2 (double x) if (w[0] == w[0] + 1.000000005 * w[1]) return w[0]; else - return __mpcos (x, 0); + return __mpcos (x, 0, false); } } @@ -1246,7 +1244,7 @@ csloww (double x, double dx, double orig) if (w[0] == w[0] + cor) return (a > 0) ? w[0] : -w[0]; else - return __mpcos1 (orig); + return __mpcos (orig, 0, true); } } } @@ -1301,7 +1299,7 @@ csloww1 (double x, double dx, double orig) if (w[0] == w[0] + cor) return (x > 0) ? w[0] : -w[0]; else - return __mpcos1 (orig); + return __mpcos (orig, 0, true); } } @@ -1357,7 +1355,7 @@ csloww2 (double x, double dx, double orig, int n) if (w[0] == w[0] + cor) return (n) ? -w[0] : w[0]; else - return __mpcos1 (orig); + return __mpcos (orig, 0, true); } } diff --git a/sysdeps/ieee754/dbl-64/sincos32.c b/sysdeps/ieee754/dbl-64/sincos32.c index ac27fcda83..f253b8ce8b 100644 --- a/sysdeps/ieee754/dbl-64/sincos32.c +++ b/sysdeps/ieee754/dbl-64/sincos32.c @@ -187,50 +187,119 @@ __cos32 (double x, double res, double res1) return (res < res1) ? res : res1; } -/* Compute sin(x+dx) as Multi Precision number and return result as double. */ +/* Compute sin() of double-length number (X + DX) as Multi Precision number and + return result as double. If REDUCE_RANGE is true, X is assumed to be the + original input and DX is ignored. */ double SECTION -__mpsin (double x, double dx) +__mpsin (double x, double dx, bool reduce_range) { - int p; double y; - mp_no a, b, c; - p = 32; - __dbl_mp (x, &a, p); - __dbl_mp (dx, &b, p); - __add (&a, &b, &c, p); - if (x > 0.8) + mp_no a, b, c, s; + int n; + int p = 32; + + if (reduce_range) { - __sub (&hp, &c, &a, p); - __c32 (&a, &b, &c, p); + n = __mpranred (x, &a, p); /* n is 0, 1, 2 or 3. */ + __c32 (&a, &c, &s, p); } else - __c32 (&c, &a, &b, p); /* b = sin(x+dx) */ - __mp_dbl (&b, &y, p); + { + n = -1; + __dbl_mp (x, &b, p); + __dbl_mp (dx, &c, p); + __add (&b, &c, &a, p); + if (x > 0.8) + { + __sub (&hp, &a, &b, p); + __c32 (&b, &s, &c, p); + } + else + __c32 (&a, &c, &s, p); /* b = sin(x+dx) */ + } + + /* Convert result based on which quarter of unit circle y is in. */ + switch (n) + { + case 1: + __mp_dbl (&c, &y, p); + break; + + case 3: + __mp_dbl (&c, &y, p); + y = -y; + break; + + case 2: + __mp_dbl (&s, &y, p); + y = -y; + break; + + /* Quadrant not set, so the result must be sin (X + DX), which is also in + S. */ + case 0: + default: + __mp_dbl (&s, &y, p); + } return y; } -/* Compute cos() of double-length number (x+dx) as Multi Precision number and - return result as double. */ +/* Compute cos() of double-length number (X + DX) as Multi Precision number and + return result as double. If REDUCE_RANGE is true, X is assumed to be the + original input and DX is ignored. */ double SECTION -__mpcos (double x, double dx) +__mpcos (double x, double dx, bool reduce_range) { - int p; double y; - mp_no a, b, c; - p = 32; - __dbl_mp (x, &a, p); - __dbl_mp (dx, &b, p); - __add (&a, &b, &c, p); - if (x > 0.8) + mp_no a, b, c, s; + int n; + int p = 32; + + if (reduce_range) { - __sub (&hp, &c, &b, p); - __c32 (&b, &c, &a, p); + n = __mpranred (x, &a, p); /* n is 0, 1, 2 or 3. */ + __c32 (&a, &c, &s, p); } else - __c32 (&c, &a, &b, p); /* a = cos(x+dx) */ - __mp_dbl (&a, &y, p); + { + n = -1; + __dbl_mp (x, &b, p); + __dbl_mp (dx, &c, p); + __add (&b, &c, &a, p); + if (x > 0.8) + { + __sub (&hp, &a, &b, p); + __c32 (&b, &s, &c, p); + } + else + __c32 (&a, &c, &s, p); /* a = cos(x+dx) */ + } + + /* Convert result based on which quarter of unit circle y is in. */ + switch (n) + { + case 1: + __mp_dbl (&s, &y, p); + y = -y; + break; + + case 3: + __mp_dbl (&s, &y, p); + break; + + case 2: + __mp_dbl (&c, &y, p); + y = -y; + break; + + /* Quadrant not set, so the result must be cos (X + DX), which is also + stored in C. */ + case 0: + default: + __mp_dbl (&c, &y, p); + } return y; } @@ -294,84 +363,3 @@ __mpranred (double x, mp_no *y, int p) return (n & 3); } } - -/* Multi-Precision sin() function subroutine, for p = 32. It is based on the - routines mpranred() and c32(). */ -double -SECTION -__mpsin1 (double x) -{ - int p; - int n; - mp_no u, s, c; - double y; - p = 32; - n = __mpranred (x, &u, p); /* n is 0, 1, 2 or 3. */ - __c32 (&u, &c, &s, p); - /* Convert result based on which quarter of unit circle y is in. */ - switch (n) - { - case 0: - __mp_dbl (&s, &y, p); - return y; - break; - - case 2: - __mp_dbl (&s, &y, p); - return -y; - break; - - case 1: - __mp_dbl (&c, &y, p); - return y; - break; - - case 3: - __mp_dbl (&c, &y, p); - return -y; - break; - } - /* Unreachable, to make the compiler happy. */ - return 0; -} - -/* Multi-Precision cos() function subroutine, for p = 32. It is based on the - routines mpranred() and c32(). */ -double -SECTION -__mpcos1 (double x) -{ - int p; - int n; - mp_no u, s, c; - double y; - - p = 32; - n = __mpranred (x, &u, p); /* n is 0, 1, 2 or 3. */ - __c32 (&u, &c, &s, p); - /* Convert result based on which quarter of unit circle y is in. */ - switch (n) - { - case 0: - __mp_dbl (&c, &y, p); - return y; - break; - - case 2: - __mp_dbl (&c, &y, p); - return -y; - break; - - case 1: - __mp_dbl (&s, &y, p); - return -y; - break; - - case 3: - __mp_dbl (&s, &y, p); - return y; - break; - } - /* Unreachable, to make the compiler happy. */ - return 0; -} -- cgit 1.4.1