This is the mail archive of the glibc-cvs@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]

GNU C Library master sources branch master updated. glibc-2.26.9000-993-g2e77dee


This is an automated email from the git hooks/post-receive script. It was
generated because a ref change was pushed to the repository containing
the project "GNU C Library master sources".

The branch, master has been updated
       via  2e77deef676677a7bed97d87d2519679654a2772 (commit)
       via  0b9bef6d2253c9af4d49a8b2d59a76f721f4cea5 (commit)
       via  984ae9967b49830173490a33ae6130880f3f70d9 (commit)
      from  93930ea9351c0c4a239e3dcb83f1398cce4e4d43 (commit)

Those revisions listed above that are new to this repository have
not appeared on any other notification email; so we list those
revisions in full, below.

- Log -----------------------------------------------------------------
http://sourceware.org/git/gitweb.cgi?p=glibc.git;a=commitdiff;h=2e77deef676677a7bed97d87d2519679654a2772

commit 2e77deef676677a7bed97d87d2519679654a2772
Author: Rajalakshmi Srinivasaraghavan <raji@linux.vnet.ibm.com>
Date:   Sat Dec 16 14:11:56 2017 +0530

    s390: Update ulps
    
    	* sysdeps/s390/fpu/libm-test-ulps: Update.

diff --git a/ChangeLog b/ChangeLog
index 02d9cdd..67d11c9 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,9 @@
 2017-12-16  Rajalakshmi Srinivasaraghavan  <raji@linux.vnet.ibm.com>
 
+	* sysdeps/s390/fpu/libm-test-ulps: Update.
+
+2017-12-16  Rajalakshmi Srinivasaraghavan  <raji@linux.vnet.ibm.com>
+
 	* sysdeps/powerpc/fpu/libm-test-ulps: Update.
 
 2017-12-16  Rajalakshmi Srinivasaraghavan  <raji@linux.vnet.ibm.com>
diff --git a/sysdeps/s390/fpu/libm-test-ulps b/sysdeps/s390/fpu/libm-test-ulps
index 5741bc7..781f46f 100644
--- a/sysdeps/s390/fpu/libm-test-ulps
+++ b/sysdeps/s390/fpu/libm-test-ulps
@@ -1325,9 +1325,9 @@ ldouble: 3
 
 Function: Imaginary part of "ctan":
 double: 2
-float: 1
+float: 2
 idouble: 2
-ifloat: 1
+ifloat: 2
 ildouble: 3
 ldouble: 3
 
@@ -1341,9 +1341,9 @@ ldouble: 4
 
 Function: Imaginary part of "ctan_downward":
 double: 2
-float: 1
+float: 2
 idouble: 2
-ifloat: 1
+ifloat: 2
 ildouble: 5
 ldouble: 5
 
@@ -1365,9 +1365,9 @@ ldouble: 5
 
 Function: Real part of "ctan_upward":
 double: 2
-float: 3
+float: 4
 idouble: 2
-ifloat: 3
+ifloat: 4
 ildouble: 5
 ldouble: 5
 
@@ -1397,9 +1397,9 @@ ldouble: 3
 
 Function: Real part of "ctanh_downward":
 double: 4
-float: 1
+float: 2
 idouble: 4
-ifloat: 1
+ifloat: 2
 ildouble: 5
 ldouble: 5
 
@@ -1691,9 +1691,9 @@ ldouble: 2
 
 Function: "j0_downward":
 double: 2
-float: 3
+float: 4
 idouble: 2
-ifloat: 3
+ifloat: 4
 ildouble: 4
 ldouble: 4
 
@@ -2157,9 +2157,9 @@ ldouble: 3
 
 Function: "y0_downward":
 double: 3
-float: 2
+float: 4
 idouble: 3
-ifloat: 2
+ifloat: 4
 ildouble: 4
 ldouble: 4
 
@@ -2173,9 +2173,9 @@ ldouble: 3
 
 Function: "y0_upward":
 double: 2
-float: 3
+float: 5
 idouble: 2
-ifloat: 3
+ifloat: 5
 ildouble: 3
 ldouble: 3
 
@@ -2213,17 +2213,17 @@ ldouble: 5
 
 Function: "yn":
 double: 3
-float: 2
+float: 3
 idouble: 3
-ifloat: 2
+ifloat: 3
 ildouble: 5
 ldouble: 5
 
 Function: "yn_downward":
 double: 3
-float: 2
+float: 4
 idouble: 3
-ifloat: 2
+ifloat: 4
 ildouble: 5
 ldouble: 5
 
@@ -2237,9 +2237,9 @@ ldouble: 5
 
 Function: "yn_upward":
 double: 4
-float: 3
+float: 5
 idouble: 4
-ifloat: 3
+ifloat: 5
 ildouble: 5
 ldouble: 5
 

http://sourceware.org/git/gitweb.cgi?p=glibc.git;a=commitdiff;h=0b9bef6d2253c9af4d49a8b2d59a76f721f4cea5

commit 0b9bef6d2253c9af4d49a8b2d59a76f721f4cea5
Author: Rajalakshmi Srinivasaraghavan <raji@linux.vnet.ibm.com>
Date:   Sat Dec 16 14:04:14 2017 +0530

    powerpc: Update ulps
    
    	* sysdeps/powerpc/fpu/libm-test-ulps: Update.

diff --git a/ChangeLog b/ChangeLog
index 82b6048..02d9cdd 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,5 +1,9 @@
 2017-12-16  Rajalakshmi Srinivasaraghavan  <raji@linux.vnet.ibm.com>
 
+	* sysdeps/powerpc/fpu/libm-test-ulps: Update.
+
+2017-12-16  Rajalakshmi Srinivasaraghavan  <raji@linux.vnet.ibm.com>
+
 	* sysdeps/ieee754/flt-32/s_cosf.c: Move reduced() and
 	constants to s_sincosf.h file.
 	* sysdeps/ieee754/flt-32/s_sinf.c: Likewise.
diff --git a/sysdeps/powerpc/fpu/libm-test-ulps b/sysdeps/powerpc/fpu/libm-test-ulps
index a685cab..5309a4a 100644
--- a/sysdeps/powerpc/fpu/libm-test-ulps
+++ b/sysdeps/powerpc/fpu/libm-test-ulps
@@ -1669,10 +1669,10 @@ ldouble: 3
 
 Function: Imaginary part of "ctan":
 double: 2
-float: 1
+float: 2
 float128: 3
 idouble: 2
-ifloat: 1
+ifloat: 2
 ifloat128: 3
 ildouble: 2
 ldouble: 2
@@ -1689,10 +1689,10 @@ ldouble: 6
 
 Function: Imaginary part of "ctan_downward":
 double: 2
-float: 1
+float: 2
 float128: 5
 idouble: 2
-ifloat: 1
+ifloat: 2
 ifloat128: 5
 ildouble: 9
 ldouble: 9
@@ -1719,10 +1719,10 @@ ldouble: 13
 
 Function: Real part of "ctan_upward":
 double: 2
-float: 3
+float: 4
 float128: 5
 idouble: 2
-ifloat: 3
+ifloat: 4
 ifloat128: 5
 ildouble: 7
 ldouble: 7
@@ -1759,10 +1759,10 @@ ldouble: 3
 
 Function: Real part of "ctanh_downward":
 double: 4
-float: 1
+float: 2
 float128: 5
 idouble: 4
-ifloat: 1
+ifloat: 2
 ifloat128: 5
 ildouble: 9
 ldouble: 9
@@ -2149,10 +2149,10 @@ ldouble: 2
 
 Function: "j0_downward":
 double: 2
-float: 3
+float: 4
 float128: 4
 idouble: 2
-ifloat: 3
+ifloat: 4
 ifloat128: 4
 ildouble: 11
 ldouble: 11
@@ -2759,10 +2759,10 @@ ldouble: 1
 
 Function: "y0_downward":
 double: 3
-float: 2
+float: 4
 float128: 4
 idouble: 3
-ifloat: 2
+ifloat: 4
 ifloat128: 4
 ildouble: 10
 ldouble: 10
@@ -2779,10 +2779,10 @@ ldouble: 8
 
 Function: "y0_upward":
 double: 2
-float: 3
+float: 5
 float128: 3
 idouble: 2
-ifloat: 3
+ifloat: 5
 ifloat128: 3
 ildouble: 9
 ldouble: 9
@@ -2829,20 +2829,20 @@ ldouble: 9
 
 Function: "yn":
 double: 3
-float: 2
+float: 3
 float128: 5
 idouble: 3
-ifloat: 2
+ifloat: 3
 ifloat128: 5
 ildouble: 2
 ldouble: 2
 
 Function: "yn_downward":
 double: 3
-float: 2
+float: 4
 float128: 5
 idouble: 3
-ifloat: 2
+ifloat: 4
 ifloat128: 5
 ildouble: 10
 ldouble: 10
@@ -2859,10 +2859,10 @@ ldouble: 8
 
 Function: "yn_upward":
 double: 4
-float: 3
+float: 5
 float128: 5
 idouble: 4
-ifloat: 3
+ifloat: 5
 ifloat128: 5
 ildouble: 9
 ldouble: 9

http://sourceware.org/git/gitweb.cgi?p=glibc.git;a=commitdiff;h=984ae9967b49830173490a33ae6130880f3f70d9

commit 984ae9967b49830173490a33ae6130880f3f70d9
Author: Rajalakshmi Srinivasaraghavan <raji@linux.vnet.ibm.com>
Date:   Sat Dec 16 14:01:37 2017 +0530

    New generic sincosf
    
    This implementation is based on generic s_sinf.c and s_cosf.c.
    Tested on s390x, powerpc64le and powerpc32.

diff --git a/ChangeLog b/ChangeLog
index 0f2936c..82b6048 100644
--- a/ChangeLog
+++ b/ChangeLog
@@ -1,3 +1,13 @@
+2017-12-16  Rajalakshmi Srinivasaraghavan  <raji@linux.vnet.ibm.com>
+
+	* sysdeps/ieee754/flt-32/s_cosf.c: Move reduced() and
+	constants to s_sincosf.h file.
+	* sysdeps/ieee754/flt-32/s_sinf.c: Likewise.
+	* sysdeps/ieee754/flt-32/s_sincosf.c: New
+	implementation.
+	* sysdeps/ieee754/flt-32/s_sincosf.h:
+	New file.
+
 2017-12-12  Carlos O'Donell <carlos@redhat.com>
 
 	[BZ #14681]
diff --git a/sysdeps/ieee754/flt-32/s_cosf.c b/sysdeps/ieee754/flt-32/s_cosf.c
index ac6d044..f0ebee2 100644
--- a/sysdeps/ieee754/flt-32/s_cosf.c
+++ b/sysdeps/ieee754/flt-32/s_cosf.c
@@ -20,6 +20,7 @@
 #include <math.h>
 #include <math_private.h>
 #include <libm-alias-float.h>
+#include "s_sincosf.h"
 
 #ifndef COSF
 # define COSF_FUNC __cosf
@@ -27,95 +28,6 @@
 # define COSF_FUNC COSF
 #endif
 
-/* Chebyshev constants for cos, range -PI/4 - PI/4.  */
-static const double C0 = -0x1.ffffffffe98aep-2;
-static const double C1 =  0x1.55555545c50c7p-5;
-static const double C2 = -0x1.6c16b348b6874p-10;
-static const double C3 =  0x1.a00eb9ac43ccp-16;
-static const double C4 = -0x1.23c97dd8844d7p-22;
-
-/* Chebyshev constants for sin, range -PI/4 - PI/4.  */
-static const double S0 = -0x1.5555555551cd9p-3;
-static const double S1 =  0x1.1111110c2688bp-7;
-static const double S2 = -0x1.a019f8b4bd1f9p-13;
-static const double S3 =  0x1.71d7264e6b5b4p-19;
-static const double S4 = -0x1.a947e1674b58ap-26;
-
-/* Chebyshev constants for cos, range 2^-27 - 2^-5.  */
-static const double CC0 = -0x1.fffffff5cc6fdp-2;
-static const double CC1 =  0x1.55514b178dac5p-5;
-
-/* PI/2 with 98 bits of accuracy.  */
-static const double PI_2_hi = 0x1.921fb544p+0;
-static const double PI_2_lo = 0x1.0b4611a626332p-34;
-
-static const double inv_PI_4 = 0x1.45f306dc9c883p+0; /* 4/PI.  */
-
-#define FLOAT_EXPONENT_SHIFT 23
-#define FLOAT_EXPONENT_BIAS 127
-
-static const double pio2_table[] = {
-  0 * M_PI_2,
-  1 * M_PI_2,
-  2 * M_PI_2,
-  3 * M_PI_2,
-  4 * M_PI_2,
-  5 * M_PI_2
-};
-
-static const double invpio4_table[] = {
-  0x0p+0,
-  0x1.45f306cp+0,
-  0x1.c9c882ap-28,
-  0x1.4fe13a8p-58,
-  0x1.f47d4dp-85,
-  0x1.bb81b6cp-112,
-  0x1.4acc9ep-142,
-  0x1.0e4107cp-169
-};
-
-static const double ones[] = { 1.0, -1.0 };
-
-/* Compute the cosine value using Chebyshev polynomials where
-   THETA is the range reduced absolute value of the input
-   and it is less than Pi/4,
-   N is calculated as trunc(|x|/(Pi/4)) + 1 and it is used to decide
-   whether a sine or cosine approximation is more accurate and
-   the sign of the result.  */
-static inline float
-reduced (double theta, unsigned int n)
-{
-  double sign, cx;
-  const double theta2 = theta * theta;
-
-  /* Determine positive or negative primary interval.  */
-  n += 2;
-  sign = ones[(n >> 2) & 1];
-
-  /* Are we in the primary interval of sin or cos?  */
-  if ((n & 2) == 0)
-    {
-      /* Here cosf() is calculated using sin Chebyshev polynomial:
-	x+x^3*(S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4)))).  */
-      cx = S3 + theta2 * S4;
-      cx = S2 + theta2 * cx;
-      cx = S1 + theta2 * cx;
-      cx = S0 + theta2 * cx;
-      cx = theta + theta * theta2 * cx;
-    }
-  else
-    {
-     /* Here cosf() is calculated using cos Chebyshev polynomial:
-	1.0+x^2*(C0+x^2*(C1+x^2*(C2+x^2*(C3+x^2*C4)))).  */
-      cx = C3 + theta2 * C4;
-      cx = C2 + theta2 * cx;
-      cx = C1 + theta2 * cx;
-      cx = C0 + theta2 * cx;
-      cx = 1. + theta2 * cx;
-    }
-  return sign * cx;
-}
-
 float
 COSF_FUNC (float x)
 {
@@ -161,7 +73,7 @@ COSF_FUNC (float x)
 	     pio2_table must go to 5 (9 / 2 + 1).  */
 	  unsigned int n = (abstheta * inv_PI_4) + 1;
 	  theta = abstheta - pio2_table[n / 2];
-	  return reduced (theta, n);
+	  return reduced_cos (theta, n);
 	}
       else if (isless (abstheta, INFINITY))
 	{
@@ -171,7 +83,7 @@ COSF_FUNC (float x)
 	      double x = n / 2;
 	      theta = (abstheta - x * PI_2_hi) - x * PI_2_lo;
 	      /* Argument reduction needed.  */
-	      return reduced (theta, n);
+	      return reduced_cos (theta, n);
 	    }
 	  else /* |theta| >= 2^23.  */
 	    {
@@ -199,7 +111,7 @@ COSF_FUNC (float x)
 		  e += c;
 		  e += d;
 		  e *= M_PI_4;
-		  return reduced (e, l + 1);
+		  return reduced_cos (e, l + 1);
 		}
 	      else
 		{
@@ -209,14 +121,14 @@ COSF_FUNC (float x)
 		  if (e <= 1.0)
 		    {
 		      e *= M_PI_4;
-		      return reduced (e, l + 1);
+		      return reduced_cos (e, l + 1);
 		    }
 		  else
 		    {
 		      l++;
 		      e -= 2.0;
 		      e *= M_PI_4;
-		      return reduced (e, l + 1);
+		      return reduced_cos (e, l + 1);
 		    }
 		}
 	    }
diff --git a/sysdeps/ieee754/flt-32/s_sincosf.c b/sysdeps/ieee754/flt-32/s_sincosf.c
index 4946a6e..c376d20 100644
--- a/sysdeps/ieee754/flt-32/s_sincosf.c
+++ b/sysdeps/ieee754/flt-32/s_sincosf.c
@@ -1,7 +1,6 @@
 /* Compute sine and cosine of argument.
-   Copyright (C) 1997-2017 Free Software Foundation, Inc.
+   Copyright (C) 2017 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
-   Contributed by Ulrich Drepper <drepper@cygnus.com>, 1997.
 
    The GNU C Library is free software; you can redistribute it and/or
    modify it under the terms of the GNU Lesser General Public
@@ -19,9 +18,9 @@
 
 #include <errno.h>
 #include <math.h>
-
 #include <math_private.h>
 #include <libm-alias-float.h>
+#include "s_sincosf.h"
 
 #ifndef SINCOSF
 # define SINCOSF_FUNC __sincosf
@@ -32,50 +31,137 @@
 void
 SINCOSF_FUNC (float x, float *sinx, float *cosx)
 {
-  int32_t ix;
-
-  /* High word of x. */
-  GET_FLOAT_WORD (ix, x);
-
-  /* |x| ~< pi/4 */
-  ix &= 0x7fffffff;
-  if (ix <= 0x3f490fd8)
-    {
-      *sinx = __kernel_sinf (x, 0.0, 0);
-      *cosx = __kernel_cosf (x, 0.0);
-    }
-  else if (ix>=0x7f800000)
+  double cx;
+  double theta = x;
+  double abstheta = fabs (theta);
+  /* If |x|< Pi/4.  */
+  if (isless (abstheta, M_PI_4))
     {
-      /* sin(Inf or NaN) is NaN */
-      *sinx = *cosx = x - x;
-      if (ix == 0x7f800000)
-	__set_errno (EDOM);
+      if (abstheta >= 0x1p-5) /* |x| >= 2^-5.  */
+	{
+	  const double theta2 = theta * theta;
+	  /* Chebyshev polynomial of the form for sin and cos.  */
+	  cx = C3 + theta2 * C4;
+	  cx = C2 + theta2 * cx;
+	  cx = C1 + theta2 * cx;
+	  cx = C0 + theta2 * cx;
+	  cx = 1.0 + theta2 * cx;
+	  *cosx = cx;
+	  cx = S3 + theta2 * S4;
+	  cx = S2 + theta2 * cx;
+	  cx = S1 + theta2 * cx;
+	  cx = S0 + theta2 * cx;
+	  cx = theta + theta * theta2 * cx;
+	  *sinx = cx;
+	}
+      else if (abstheta >= 0x1p-27)     /* |x| >= 2^-27.  */
+	{
+	  /* A simpler Chebyshev approximation is close enough for this range:
+	     for sin: x+x^3*(SS0+x^2*SS1)
+	     for cos: 1.0+x^2*(CC0+x^3*CC1).  */
+	  const double theta2 = theta * theta;
+	  cx = CC0 + theta * theta2 * CC1;
+	  cx = 1.0 + theta2 * cx;
+	  *cosx = cx;
+	  cx = SS0 + theta2 * SS1;
+	  cx = theta + theta * theta2 * cx;
+	  *sinx = cx;
+	}
+      else
+	{
+	  /* Handle some special cases.  */
+	  if (theta)
+	    *sinx = theta - (theta * SMALL);
+	  else
+	    *sinx = theta;
+	  *cosx = 1.0 - abstheta;
+	}
     }
-  else
+  else                          /* |x| >= Pi/4.  */
     {
-      /* Argument reduction needed.  */
-      float y[2];
-      int n;
-
-      n = __ieee754_rem_pio2f (x, y);
-      switch (n & 3)
+      unsigned int signbit = isless (x, 0);
+      if (isless (abstheta, 9 * M_PI_4))        /* |x| < 9*Pi/4.  */
+	{
+	  /* There are cases where FE_UPWARD rounding mode can
+	     produce a result of abstheta * inv_PI_4 == 9,
+	     where abstheta < 9pi/4, so the domain for
+	     pio2_table must go to 5 (9 / 2 + 1).  */
+	  unsigned int n = (abstheta * inv_PI_4) + 1;
+	  theta = abstheta - pio2_table[n / 2];
+	  *sinx = reduced_sin (theta, n, signbit);
+	  *cosx = reduced_cos (theta, n);
+	}
+      else if (isless (abstheta, INFINITY))
+	{
+	  if (abstheta < 0x1p+23)     /* |x| < 2^23.  */
+	    {
+	      unsigned int n = ((unsigned int) (abstheta * inv_PI_4)) + 1;
+	      double x = n / 2;
+	      theta = (abstheta - x * PI_2_hi) - x * PI_2_lo;
+	      /* Argument reduction needed.  */
+	      *sinx = reduced_sin (theta, n, signbit);
+	      *cosx = reduced_cos (theta, n);
+	    }
+	  else                  /* |x| >= 2^23.  */
+	    {
+	      x = fabsf (x);
+	      int exponent;
+	      GET_FLOAT_WORD (exponent, x);
+	      exponent
+	        = (exponent >> FLOAT_EXPONENT_SHIFT) - FLOAT_EXPONENT_BIAS;
+	      exponent += 3;
+	      exponent /= 28;
+	      double a = invpio4_table[exponent] * x;
+	      double b = invpio4_table[exponent + 1] * x;
+	      double c = invpio4_table[exponent + 2] * x;
+	      double d = invpio4_table[exponent + 3] * x;
+	      uint64_t l = a;
+	      l &= ~0x7;
+	      a -= l;
+	      double e = a + b;
+	      l = e;
+	      e = a - l;
+	      if (l & 1)
+	        {
+	          e -= 1.0;
+	          e += b;
+	          e += c;
+	          e += d;
+	          e *= M_PI_4;
+		  *sinx = reduced_sin (e, l + 1, signbit);
+		  *cosx = reduced_cos (e, l + 1);
+	        }
+	      else
+		{
+		  e += b;
+		  e += c;
+		  e += d;
+		  if (e <= 1.0)
+		    {
+		      e *= M_PI_4;
+		      *sinx = reduced_sin (e, l + 1, signbit);
+		      *cosx = reduced_cos (e, l + 1);
+		    }
+		  else
+		    {
+		      l++;
+		      e -= 2.0;
+		      e *= M_PI_4;
+		      *sinx = reduced_sin (e, l + 1, signbit);
+		      *cosx = reduced_cos (e, l + 1);
+		    }
+		}
+	    }
+	}
+      else
 	{
-	case 0:
-	  *sinx = __kernel_sinf (y[0], y[1], 1);
-	  *cosx = __kernel_cosf (y[0], y[1]);
-	  break;
-	case 1:
-	  *sinx = __kernel_cosf (y[0], y[1]);
-	  *cosx = -__kernel_sinf (y[0], y[1], 1);
-	  break;
-	case 2:
-	  *sinx = -__kernel_sinf (y[0], y[1], 1);
-	  *cosx = -__kernel_cosf (y[0], y[1]);
-	  break;
-	default:
-	  *sinx = -__kernel_cosf (y[0], y[1]);
-	  *cosx = __kernel_sinf (y[0], y[1], 1);
-	  break;
+	  int32_t ix;
+	  /* High word of x.  */
+	  GET_FLOAT_WORD (ix, abstheta);
+	  /* sin/cos(Inf or NaN) is NaN.  */
+	  *sinx = *cosx = x - x;
+	  if (ix == 0x7f800000)
+	    __set_errno (EDOM);
 	}
     }
 }
diff --git a/sysdeps/ieee754/flt-32/s_cosf.c b/sysdeps/ieee754/flt-32/s_sincosf.h
similarity index 53%
copy from sysdeps/ieee754/flt-32/s_cosf.c
copy to sysdeps/ieee754/flt-32/s_sincosf.h
index ac6d044..b0110fc 100644
--- a/sysdeps/ieee754/flt-32/s_cosf.c
+++ b/sysdeps/ieee754/flt-32/s_sincosf.h
@@ -1,4 +1,4 @@
-/* Compute cosine of argument.
+/* Used by sinf, cosf and sincosf functions.
    Copyright (C) 2017 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
 
@@ -16,17 +16,6 @@
    License along with the GNU C Library; if not, see
    <http://www.gnu.org/licenses/>.  */
 
-#include <errno.h>
-#include <math.h>
-#include <math_private.h>
-#include <libm-alias-float.h>
-
-#ifndef COSF
-# define COSF_FUNC __cosf
-#else
-# define COSF_FUNC COSF
-#endif
-
 /* Chebyshev constants for cos, range -PI/4 - PI/4.  */
 static const double C0 = -0x1.ffffffffe98aep-2;
 static const double C1 =  0x1.55555545c50c7p-5;
@@ -41,6 +30,10 @@ static const double S2 = -0x1.a019f8b4bd1f9p-13;
 static const double S3 =  0x1.71d7264e6b5b4p-19;
 static const double S4 = -0x1.a947e1674b58ap-26;
 
+/* Chebyshev constants for sin, range 2^-27 - 2^-5.  */
+static const double SS0 = -0x1.555555543d49dp-3;
+static const double SS1 =  0x1.110f475cec8c5p-7;
+
 /* Chebyshev constants for cos, range 2^-27 - 2^-5.  */
 static const double CC0 = -0x1.fffffff5cc6fdp-2;
 static const double CC1 =  0x1.55514b178dac5p-5;
@@ -49,6 +42,7 @@ static const double CC1 =  0x1.55514b178dac5p-5;
 static const double PI_2_hi = 0x1.921fb544p+0;
 static const double PI_2_lo = 0x1.0b4611a626332p-34;
 
+static const double SMALL = 0x1p-50; /* 2^-50.  */
 static const double inv_PI_4 = 0x1.45f306dc9c883p+0; /* 4/PI.  */
 
 #define FLOAT_EXPONENT_SHIFT 23
@@ -76,6 +70,50 @@ static const double invpio4_table[] = {
 
 static const double ones[] = { 1.0, -1.0 };
 
+/* Compute the sine value using Chebyshev polynomials where
+   THETA is the range reduced absolute value of the input
+   and it is less than Pi/4,
+   N is calculated as trunc(|x|/(Pi/4)) + 1 and it is used to decide
+   whether a sine or cosine approximation is more accurate and
+   SIGNBIT is used to add the correct sign after the Chebyshev
+   polynomial is computed.  */
+static inline float
+reduced_sin (const double theta, const unsigned int n,
+	 const unsigned int signbit)
+{
+  double sx;
+  const double theta2 = theta * theta;
+  /* We are operating on |x|, so we need to add back the original
+     signbit for sinf.  */
+  double sign;
+  /* Determine positive or negative primary interval.  */
+  sign = ones[((n >> 2) & 1) ^ signbit];
+  /* Are we in the primary interval of sin or cos?  */
+  if ((n & 2) == 0)
+    {
+      /* Here sinf() is calculated using sin Chebyshev polynomial:
+	x+x^3*(S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4)))).  */
+      sx = S3 + theta2 * S4;     /* S3+x^2*S4.  */
+      sx = S2 + theta2 * sx;     /* S2+x^2*(S3+x^2*S4).  */
+      sx = S1 + theta2 * sx;     /* S1+x^2*(S2+x^2*(S3+x^2*S4)).  */
+      sx = S0 + theta2 * sx;     /* S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4))).  */
+      sx = theta + theta * theta2 * sx;
+    }
+  else
+    {
+     /* Here sinf() is calculated using cos Chebyshev polynomial:
+	1.0+x^2*(C0+x^2*(C1+x^2*(C2+x^2*(C3+x^2*C4)))).  */
+      sx = C3 + theta2 * C4;     /* C3+x^2*C4.  */
+      sx = C2 + theta2 * sx;     /* C2+x^2*(C3+x^2*C4).  */
+      sx = C1 + theta2 * sx;     /* C1+x^2*(C2+x^2*(C3+x^2*C4)).  */
+      sx = C0 + theta2 * sx;     /* C0+x^2*(C1+x^2*(C2+x^2*(C3+x^2*C4))).  */
+      sx = 1.0 + theta2 * sx;
+    }
+
+  /* Add in the signbit and assign the result.  */
+  return sign * sx;
+}
+
 /* Compute the cosine value using Chebyshev polynomials where
    THETA is the range reduced absolute value of the input
    and it is less than Pi/4,
@@ -83,7 +121,7 @@ static const double ones[] = { 1.0, -1.0 };
    whether a sine or cosine approximation is more accurate and
    the sign of the result.  */
 static inline float
-reduced (double theta, unsigned int n)
+reduced_cos (double theta, unsigned int n)
 {
   double sign, cx;
   const double theta2 = theta * theta;
@@ -115,124 +153,3 @@ reduced (double theta, unsigned int n)
     }
   return sign * cx;
 }
-
-float
-COSF_FUNC (float x)
-{
-  double theta = x;
-  double abstheta = fabs (theta);
-  if (isless (abstheta, M_PI_4))
-    {
-      double cx;
-      if (abstheta >= 0x1p-5)
-	{
-	  const double theta2 = theta * theta;
-	  /* Chebyshev polynomial of the form for cos:
-	   * 1 + x^2 (C0 + x^2 (C1 + x^2 (C2 + x^2 (C3 + x^2 * C4)))).  */
-	  cx = C3 + theta2 * C4;
-	  cx = C2 + theta2 * cx;
-	  cx = C1 + theta2 * cx;
-	  cx = C0 + theta2 * cx;
-	  cx = 1. + theta2 * cx;
-	  return cx;
-	}
-      else if (abstheta >= 0x1p-27)
-	{
-	  /* A simpler Chebyshev approximation is close enough for this range:
-	   * 1 + x^2 (CC0 + x^3 * CC1).  */
-	  const double theta2 = theta * theta;
-	  cx = CC0 + theta * theta2 * CC1;
-	  cx = 1.0 + theta2 * cx;
-	  return cx;
-	}
-      else
-	{
-	  /* For small enough |theta|, this is close enough.  */
-	  return 1.0 - abstheta;
-	}
-    }
-  else /* |theta| >= Pi/4.  */
-    {
-      if (isless (abstheta, 9 * M_PI_4))
-	{
-	  /* There are cases where FE_UPWARD rounding mode can
-	     produce a result of abstheta * inv_PI_4 == 9,
-	     where abstheta < 9pi/4, so the domain for
-	     pio2_table must go to 5 (9 / 2 + 1).  */
-	  unsigned int n = (abstheta * inv_PI_4) + 1;
-	  theta = abstheta - pio2_table[n / 2];
-	  return reduced (theta, n);
-	}
-      else if (isless (abstheta, INFINITY))
-	{
-	  if (abstheta < 0x1p+23)
-	    {
-	      unsigned int n = ((unsigned int) (abstheta * inv_PI_4)) + 1;
-	      double x = n / 2;
-	      theta = (abstheta - x * PI_2_hi) - x * PI_2_lo;
-	      /* Argument reduction needed.  */
-	      return reduced (theta, n);
-	    }
-	  else /* |theta| >= 2^23.  */
-	    {
-	      x = fabsf (x);
-	      int exponent;
-	      GET_FLOAT_WORD (exponent, x);
-	      exponent = (exponent >> FLOAT_EXPONENT_SHIFT)
-			 - FLOAT_EXPONENT_BIAS;
-	      exponent += 3;
-	      exponent /= 28;
-	      double a = invpio4_table[exponent] * x;
-	      double b = invpio4_table[exponent + 1] * x;
-	      double c = invpio4_table[exponent + 2] * x;
-	      double d = invpio4_table[exponent + 3] * x;
-	      uint64_t l = a;
-	      l &= ~0x7;
-	      a -= l;
-	      double e = a + b;
-	      l = e;
-	      e = a - l;
-	      if (l & 1)
-		{
-		  e -= 1.0;
-		  e += b;
-		  e += c;
-		  e += d;
-		  e *= M_PI_4;
-		  return reduced (e, l + 1);
-		}
-	      else
-		{
-		  e += b;
-		  e += c;
-		  e += d;
-		  if (e <= 1.0)
-		    {
-		      e *= M_PI_4;
-		      return reduced (e, l + 1);
-		    }
-		  else
-		    {
-		      l++;
-		      e -= 2.0;
-		      e *= M_PI_4;
-		      return reduced (e, l + 1);
-		    }
-		}
-	    }
-	}
-      else
-	{
-	  int32_t ix;
-	  GET_FLOAT_WORD (ix, abstheta);
-	  /* cos(Inf or NaN) is NaN.  */
-	  if (ix == 0x7f800000) /* Inf.  */
-	    __set_errno (EDOM);
-	  return x - x;
-	}
-    }
-}
-
-#ifndef COSF
-libm_alias_float (__cos, cos)
-#endif
diff --git a/sysdeps/ieee754/flt-32/s_sinf.c b/sysdeps/ieee754/flt-32/s_sinf.c
index 418d448..1fd1fd1 100644
--- a/sysdeps/ieee754/flt-32/s_sinf.c
+++ b/sysdeps/ieee754/flt-32/s_sinf.c
@@ -20,6 +20,7 @@
 #include <math.h>
 #include <math_private.h>
 #include <libm-alias-float.h>
+#include "s_sincosf.h"
 
 #ifndef SINF
 # define SINF_FUNC __sinf
@@ -27,100 +28,6 @@
 # define SINF_FUNC SINF
 #endif
 
-/* Chebyshev constants for cos, range -PI/4 - PI/4.  */
-static const double C0 = -0x1.ffffffffe98aep-2;
-static const double C1 =  0x1.55555545c50c7p-5;
-static const double C2 = -0x1.6c16b348b6874p-10;
-static const double C3 =  0x1.a00eb9ac43ccp-16;
-static const double C4 = -0x1.23c97dd8844d7p-22;
-
-/* Chebyshev constants for sin, range -PI/4 - PI/4.  */
-static const double S0 = -0x1.5555555551cd9p-3;
-static const double S1 =  0x1.1111110c2688bp-7;
-static const double S2 = -0x1.a019f8b4bd1f9p-13;
-static const double S3 =  0x1.71d7264e6b5b4p-19;
-static const double S4 = -0x1.a947e1674b58ap-26;
-
-/* Chebyshev constants for sin, range 2^-27 - 2^-5.  */
-static const double SS0 = -0x1.555555543d49dp-3;
-static const double SS1 =  0x1.110f475cec8c5p-7;
-
-/* PI/2 with 98 bits of accuracy.  */
-static const double PI_2_hi = -0x1.921fb544p+0;
-static const double PI_2_lo = -0x1.0b4611a626332p-34;
-
-static const double SMALL = 0x1p-50; /* 2^-50.  */
-static const double inv_PI_4 = 0x1.45f306dc9c883p+0; /* 4/PI.  */
-
-#define FLOAT_EXPONENT_SHIFT 23
-#define FLOAT_EXPONENT_BIAS 127
-
-static const double pio2_table[] = {
-  0 * M_PI_2,
-  1 * M_PI_2,
-  2 * M_PI_2,
-  3 * M_PI_2,
-  4 * M_PI_2,
-  5 * M_PI_2
-};
-
-static const double invpio4_table[] = {
-  0x0p+0,
-  0x1.45f306cp+0,
-  0x1.c9c882ap-28,
-  0x1.4fe13a8p-58,
-  0x1.f47d4dp-85,
-  0x1.bb81b6cp-112,
-  0x1.4acc9ep-142,
-  0x1.0e4107cp-169
-};
-
-static const double ones[] = { 1.0, -1.0 };
-
-/* Compute the sine value using Chebyshev polynomials where
-   THETA is the range reduced absolute value of the input
-   and it is less than Pi/4,
-   N is calculated as trunc(|x|/(Pi/4)) + 1 and it is used to decide
-   whether a sine or cosine approximation is more accurate and
-   SIGNBIT is used to add the correct sign after the Chebyshev
-   polynomial is computed.  */
-static inline float
-reduced (const double theta, const unsigned int n,
-	 const unsigned int signbit)
-{
-  double sx;
-  const double theta2 = theta * theta;
-  /* We are operating on |x|, so we need to add back the original
-     signbit for sinf.  */
-  double sign;
-  /* Determine positive or negative primary interval.  */
-  sign = ones[((n >> 2) & 1) ^ signbit];
-  /* Are we in the primary interval of sin or cos?  */
-  if ((n & 2) == 0)
-    {
-      /* Here sinf() is calculated using sin Chebyshev polynomial:
-	x+x^3*(S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4)))).  */
-      sx = S3 + theta2 * S4;     /* S3+x^2*S4.  */
-      sx = S2 + theta2 * sx;     /* S2+x^2*(S3+x^2*S4).  */
-      sx = S1 + theta2 * sx;     /* S1+x^2*(S2+x^2*(S3+x^2*S4)).  */
-      sx = S0 + theta2 * sx;     /* S0+x^2*(S1+x^2*(S2+x^2*(S3+x^2*S4))).  */
-      sx = theta + theta * theta2 * sx;
-    }
-  else
-    {
-     /* Here sinf() is calculated using cos Chebyshev polynomial:
-	1.0+x^2*(C0+x^2*(C1+x^2*(C2+x^2*(C3+x^2*C4)))).  */
-      sx = C3 + theta2 * C4;     /* C3+x^2*C4.  */
-      sx = C2 + theta2 * sx;     /* C2+x^2*(C3+x^2*C4).  */
-      sx = C1 + theta2 * sx;     /* C1+x^2*(C2+x^2*(C3+x^2*C4)).  */
-      sx = C0 + theta2 * sx;     /* C0+x^2*(C1+x^2*(C2+x^2*(C3+x^2*C4))).  */
-      sx = 1.0 + theta2 * sx;
-    }
-
-  /* Add in the signbit and assign the result.  */
-  return sign * sx;
-}
-
 float
 SINF_FUNC (float x)
 {
@@ -171,7 +78,7 @@ SINF_FUNC (float x)
 	     pio2_table must go to 5 (9 / 2 + 1).  */
 	  unsigned int n = (abstheta * inv_PI_4) + 1;
 	  theta = abstheta - pio2_table[n / 2];
-	  return reduced (theta, n, signbit);
+	  return reduced_sin (theta, n, signbit);
 	}
       else if (isless (abstheta, INFINITY))
 	{
@@ -179,9 +86,9 @@ SINF_FUNC (float x)
 	    {
 	      unsigned int n = ((unsigned int) (abstheta * inv_PI_4)) + 1;
 	      double x = n / 2;
-	      theta = x * PI_2_lo + (x * PI_2_hi + abstheta);
+	      theta = (abstheta - x * PI_2_hi) - x * PI_2_lo;
 	      /* Argument reduction needed.  */
-	      return reduced (theta, n, signbit);
+	      return reduced_sin (theta, n, signbit);
 	    }
 	  else                  /* |x| >= 2^23.  */
 	    {
@@ -209,7 +116,7 @@ SINF_FUNC (float x)
 	          e += c;
 	          e += d;
 	          e *= M_PI_4;
-	          return reduced (e, l + 1, signbit);
+	          return reduced_sin (e, l + 1, signbit);
 	        }
 	      else
 		{
@@ -219,14 +126,14 @@ SINF_FUNC (float x)
 		  if (e <= 1.0)
 		    {
 		      e *= M_PI_4;
-		      return reduced (e, l + 1, signbit);
+		      return reduced_sin (e, l + 1, signbit);
 		    }
 		  else
 		    {
 		      l++;
 		      e -= 2.0;
 		      e *= M_PI_4;
-		      return reduced (e, l + 1, signbit);
+		      return reduced_sin (e, l + 1, signbit);
 		    }
 		}
 	    }

-----------------------------------------------------------------------

Summary of changes:
 ChangeLog                                        |   18 ++
 sysdeps/ieee754/flt-32/s_cosf.c                  |  100 +-----------
 sysdeps/ieee754/flt-32/s_sincosf.c               |  172 +++++++++++++++-----
 sysdeps/ieee754/flt-32/{s_cosf.c => s_sincosf.h} |  185 ++++++----------------
 sysdeps/ieee754/flt-32/s_sinf.c                  |  107 +------------
 sysdeps/powerpc/fpu/libm-test-ulps               |   40 +++---
 sysdeps/s390/fpu/libm-test-ulps                  |   40 +++---
 7 files changed, 251 insertions(+), 411 deletions(-)
 copy sysdeps/ieee754/flt-32/{s_cosf.c => s_sincosf.h} (53%)


hooks/post-receive
-- 
GNU C Library master sources


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