about summary refs log tree commit diff
path: root/sysdeps/ieee754/flt-32/s_sinf.c
diff options
context:
space:
mode:
Diffstat (limited to 'sysdeps/ieee754/flt-32/s_sinf.c')
-rw-r--r--sysdeps/ieee754/flt-32/s_sinf.c107
1 files changed, 7 insertions, 100 deletions
diff --git a/sysdeps/ieee754/flt-32/s_sinf.c b/sysdeps/ieee754/flt-32/s_sinf.c
index 418d4487c5..1fd1fd17d9 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);
 		    }
 		}
 	    }