Consolidate multiple precision sin/cos functions

This commit is contained in:
Siddhesh Poyarekar 2013-10-08 11:50:17 +05:30
parent 7602d070dc
commit 09544cbcd6
4 changed files with 140 additions and 129 deletions

View File

@ -1,3 +1,30 @@
2013-10-08 Siddhesh Poyarekar <siddhesh@redhat.com>
* sysdeps/generic/math_private.h (__mpsin1): Remove
declaration.
(__mpcos1): Likewise.
(__mpsin): New argument __range_reduce.
(__mpcos): Likewise.
* sysdeps/ieee754/dbl-64/s_sin.c: Likewise.
(slow): Use __mpsin and __mpcos.
(slow1): Likewise.
(slow2): Likewise.
(sloww): Likewise.
(sloww1): Likewise.
(sloww2): Likewise.
(bsloww): Likewise.
(bsloww1): Likewise.
(bsloww2): Likewise.
(cslow2): Likewise.
(csloww): Likewise.
(csloww1): Likewise.
(csloww2): Likewise.
* sysdeps/ieee754/dbl-64/sincos32.c (__mpsin): Add argument
range_reduce. Merge in __mpsin1.
(__mpcos): Likewise.
(__mpsin1): Remove.
(__mpcos1): Likewise.
2013-10-07 Joseph Myers <joseph@codesourcery.com>
* locale/loadlocale.c (_nl_intern_locale_data): Use

View File

@ -356,10 +356,8 @@ extern void __dubcos (double __x, double __dx, double __v[]);
extern double __halfulp (double __x, double __y);
extern double __sin32 (double __x, double __res, double __res1);
extern double __cos32 (double __x, double __res, double __res1);
extern double __mpsin (double __x, double __dx);
extern double __mpcos (double __x, double __dx);
extern double __mpsin1 (double __x);
extern double __mpcos1 (double __x);
extern double __mpsin (double __x, double __dx, bool __range_reduce);
extern double __mpcos (double __x, double __dx, bool __range_reduce);
extern double __slowexp (double __x);
extern double __slowpow (double __x, double __y, double __z);
extern void __docos (double __x, double __dx, double __v[]);

View File

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

View File

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