This is the mail archive of the libc-alpha@sourceware.org mailing list for the glibc project.


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]
Other format: [Raw text]

[ping][PATCH] Consolidate multiple precision sin and cos functions


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


Index Nav: [Date Index] [Subject Index] [Author Index] [Thread Index]
Message Nav: [Date Prev] [Date Next] [Thread Prev] [Thread Next]