[ping][PATCH] Consolidate multiple precision sin and cos functions
Siddhesh Poyarekar
siddhesh@redhat.com
Mon Oct 7 07:00:00 GMT 2013
Ping!
On Thu, Sep 19, 2013 at 02:25:29PM +0530, Siddhesh Poyarekar wrote:
> Hi,
>
> The multiple precision functions for sin and cos have two variants,
> one for the x + dx for and another which needs range reduction. The
> following patch consolidates the two into a single function to allow
> addition of a single systemtap probe marker. Verified that it does
> not cause any noticeable performance drop on x86_64 and i686 and also
> does not regress the testsuite. OK to commit?
>
> Siddhesh
>
> * 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.
>
>
> diff --git a/sysdeps/generic/math_private.h b/sysdeps/generic/math_private.h
> index c0fc03d..9b881a3 100644
> --- a/sysdeps/generic/math_private.h
> +++ b/sysdeps/generic/math_private.h
> @@ -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[]);
> diff --git a/sysdeps/ieee754/dbl-64/s_sin.c b/sysdeps/ieee754/dbl-64/s_sin.c
> index 7f614bd..84253a9 100644
> --- a/sysdeps/ieee754/dbl-64/s_sin.c
> +++ b/sysdeps/ieee754/dbl-64/s_sin.c
> @@ -113,10 +113,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);
> @@ -708,7 +706,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);
> }
> }
>
> @@ -748,7 +746,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);
> }
> }
>
> @@ -801,7 +799,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);
> }
> }
>
> @@ -868,7 +866,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);
> }
> }
> }
> @@ -925,7 +923,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);
> }
> }
>
> @@ -982,7 +980,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);
> }
> }
>
> @@ -1014,7 +1012,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);
> }
> }
>
> @@ -1065,7 +1063,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);
> }
> }
>
> @@ -1117,7 +1115,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);
> }
> }
>
> @@ -1159,7 +1157,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);
> }
> }
>
> @@ -1232,7 +1230,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);
> }
> }
> }
> @@ -1287,7 +1285,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);
> }
> }
>
> @@ -1343,7 +1341,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 ac27fcd..f253b8c 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;
> -}
More information about the Libc-alpha
mailing list