Patchwork RFC: Proposal for implementing libmvec on aarch64

login
register
mail settings
Submitter Steve Ellcey
Date March 6, 2018, 9:33 p.m.
Message ID <1520371988.6774.82.camel@cavium.com>
Download mbox | patch
Permalink /patch/26221/
State New
Headers show

Comments

Steve Ellcey - March 6, 2018, 9:33 p.m.
This is an RFC for an approach to implementing libmvec vector math routines
for aarch64.  There are three patches attached to this email, one is a GCC
patch that allows aarch64 to generate calls to vector functions that are
marked with __attribute__ ((simd)), the second is the infrastructure part
of glibc changes that mark sin and cos with this attribute in the header files,
and the third is an implementation of vector sin and cos by adding ifdefs to
the existing scalar sin/cos functions in sysdeps/ieee754/dbl-64/s_sin.c.

I have several questions that I would like to get some feedback on:

Does the approach of modifying the existing scalar math routines to also
work on vector types, as opposed to having completely seperate vector
routines sound like a good approach?  I like it because it reduced the
amount of duplicated code between scalar and vector functions and should
keep the accuracy of scalar and vector routines in sync.  And since it is
written in C it can also be shared among multiple platforms.  The main
downside right now is the slow paths are still mostly handled serially
and so I am only seeing about a 5% improvement in a test that does sin
calls on an array of doubles.  I think more changes could be done to improve
the amount of parallelism in the vector versions but this is a start.
Perhaps we could remove some of the slow paths in the vector versions like
we do with the sincos function.

The current implementation does not build because when linking libmvec
it gets undefined externals for __branred, __docos, __dubsin, __mpcos,
__mpsin, and __sincostab.  How should this be handled?  Should they be
exported from libm so that libmvec can call them or should they be 
copied into libmvec so that they don't have to be externally visible?
Are there any issues with making libmvec depend on libm?

The vector ABI for aarch64 is not yet defined, I believe ARM is working
on that.  As you can see from the GCC patch I am just using the same
approach as x86 and just changed the vecsize_mangle variable to 'n' for
128 bit vectors on Aarch64.

Comments?

Steve Ellcey
sellcey@cavium.com

Patch

diff --git a/sysdeps/aarch64/fpu/libmvec_d_sincos2.c b/sysdeps/aarch64/fpu/libmvec_d_sincos2.c
index e69de29..f92291f 100644
--- a/sysdeps/aarch64/fpu/libmvec_d_sincos2.c
+++ b/sysdeps/aarch64/fpu/libmvec_d_sincos2.c
@@ -0,0 +1,21 @@ 
+
+/* We need the includes before the defines so that the define of __sin
+   and __cos do not mess up the bits/mathcalls.h declaration of __sin
+   and __cos.  */
+
+#include <float.h>
+#include "endian.h"
+#include "mydefs.h"
+#include "usncs.h"
+#include "MathLib.h"
+#include <math.h>
+#include <math_private.h>
+#include <libm-alias-double.h>
+#include <fenv.h>
+
+#define VEC_SIZE 2
+#define OP_TYPE __Float64x2_t
+#define __sin _ZGVnN2v_sin
+#define __cos _ZGVnN2v_cos
+
+#include "sysdeps/ieee754/dbl-64/s_sin.c"
diff --git a/sysdeps/ieee754/dbl-64/s_sin.c b/sysdeps/ieee754/dbl-64/s_sin.c
index 8c589cb..e1db599 100644
--- a/sysdeps/ieee754/dbl-64/s_sin.c
+++ b/sysdeps/ieee754/dbl-64/s_sin.c
@@ -55,6 +55,18 @@ 
 #include <libm-alias-double.h>
 #include <fenv.h>
 
+#ifndef VEC_SIZE
+# define VEC_SIZE 1
+#endif
+
+#if VEC_SIZE != 1 && VEC_SIZE != 2
+# error "This file must be updated to support other vector lengths."
+#endif
+
+#ifndef OP_TYPE
+# define OP_TYPE double
+#endif
+
 /* Helper macros to compute sin of the input values.  */
 #define POLYNOMIAL2(xx) ((((s5 * (xx) + s4) * (xx) + s3) * (xx) + s2) * (xx))
 
@@ -67,13 +79,25 @@ 
 
    The constants s1, s2, s3, etc. are pre-computed values of 1/3!, 1/5! and so
    on.  The result is returned to LHS and correction in COR.  */
-#define TAYLOR_SIN(xx, a, da, cor) \
-({									      \
-  double t = ((POLYNOMIAL (xx)  * (a) - 0.5 * (da))  * (xx) + (da));	      \
-  double res = (a) + t;							      \
-  (cor) = ((a) - res) + t;						      \
-  res;									      \
-})
+static inline double
+__always_inline
+taylor_sin_scalar (double xx, double a, double da, double *cor)
+{
+  double t = (POLYNOMIAL (xx) * a - 0.5 * da)  * xx + da;
+  double res = a + t;
+  *cor = (a - res) + t;
+  return res;
+}
+
+static inline OP_TYPE
+__always_inline
+taylor_sin_vector (OP_TYPE xx, OP_TYPE a, OP_TYPE da, OP_TYPE *cor)
+{
+  OP_TYPE t = (POLYNOMIAL (xx) * a - 0.5 * da)  * xx + da;
+  OP_TYPE res = a + t;
+  *cor = (a - res) + t;
+  return res;
+}
 
 /* This is again a variation of the Taylor series expansion with the term
    x^3/3! expanded into the following for better accuracy:
@@ -126,10 +150,11 @@  static const double
 
 static const double t22 = 0x1.8p22;
 
-void __dubsin (double x, double dx, double w[]);
-void __docos (double x, double dx, double w[]);
-double __mpsin (double x, double dx, bool reduce_range);
-double __mpcos (double x, double dx, bool reduce_range);
+extern void __dubsin (double x, double dx, double w[]);
+extern void __docos (double x, double dx, double w[]);
+extern double __mpsin (double x, double dx, bool reduce_range);
+extern double __mpcos (double x, double dx, bool reduce_range);
+
 static double slow (double x);
 static double slow1 (double x);
 static double slow2 (double x);
@@ -148,7 +173,7 @@  static double cslow2 (double x);
    get the result in RES and a correction value in COR.  */
 static inline double
 __always_inline
-do_cos (double x, double dx, double *corp)
+do_cos_scalar (double x, double dx, double *corp)
 {
   mynumber u;
 
@@ -170,6 +195,45 @@  do_cos (double x, double dx, double *corp)
   return res;
 }
 
+static inline OP_TYPE
+__always_inline
+do_cos_vector (OP_TYPE x, OP_TYPE dx, OP_TYPE *corp)
+{
+#if VEC_SIZE == 1
+  mynumber u;
+  if (x < 0)
+    dx = -dx;
+  u.x = big + fabs (x);
+  x = fabs (x) - (u.x - big) + dx;
+#else
+  mynumber u0, u1;
+  if (x[0] < 0)
+    dx[0] = -dx[0];
+  if (x[1] < 0)
+    dx[1] = -dx[1];
+  u0.x = big + fabs (x[0]);
+  u1.x = big + fabs (x[1]);
+  x[0] = fabs (x[0]) - (u0.x - big) + dx[0];
+  x[1] = fabs (x[1]) - (u1.x - big) + dx[1];
+#endif
+
+  OP_TYPE xx, s, sn, ssn, c, cs, ccs, res, cor;
+  xx = x * x;
+  s = x + x * xx * (sn3 + xx * sn5);
+  c = xx * (cs2 + xx * (cs4 + xx * cs6));
+#if VEC_SIZE == 1
+  SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
+#else
+  SINCOS_TABLE_LOOKUP (u0, sn[0], ssn[0], cs[0], ccs[0]);
+  SINCOS_TABLE_LOOKUP (u1, sn[1], ssn[1], cs[1], ccs[1]);
+#endif
+  cor = (ccs - s * ssn - cs * c) - sn * s;
+  res = cs + cor;
+  cor = (cs - res) + cor;
+  *corp = cor;
+  return res;
+}
+
 /* A more precise variant of DO_COS.  EPS is the adjustment to the correction
    COR.  */
 static inline double
@@ -210,7 +274,7 @@  do_cos_slow (double x, double dx, double eps, double *corp)
    the result in RES and a correction value in COR.  */
 static inline double
 __always_inline
-do_sin (double x, double dx, double *corp)
+do_sin_scalar (double x, double dx, double *corp)
 {
   mynumber u;
 
@@ -231,6 +295,45 @@  do_sin (double x, double dx, double *corp)
   return res;
 }
 
+static inline OP_TYPE
+__always_inline
+do_sin_vector (OP_TYPE x, OP_TYPE dx, OP_TYPE *corp)
+{
+#if VEC_SIZE == 1
+  mynumber u;
+  if (x <= 0)
+    dx = -dx;
+  u.x = big + fabs (x);
+  x = fabs (x) - (u.x - big);
+#else
+  mynumber u0, u1;
+  if (x[0] <= 0)
+    dx[0] = -dx[0];
+  if (x[1] <= 0)
+    dx[1] = -dx[1];
+  u0.x = big + fabs (x[0]);
+  u1.x = big + fabs (x[1]);
+  x[0] = fabs (x[0]) - (u0.x - big);
+  x[1] = fabs (x[1]) - (u1.x - big);
+#endif
+
+  OP_TYPE xx, s, sn, ssn, c, cs, ccs, cor, res;
+  xx = x * x;
+  s = x + (dx + x * xx * (sn3 + xx * sn5));
+  c = x * dx + xx * (cs2 + xx * (cs4 + xx * cs6));
+#if VEC_SIZE == 1
+  SINCOS_TABLE_LOOKUP (u, sn, ssn, cs, ccs);
+#else
+  SINCOS_TABLE_LOOKUP (u0, sn[0], ssn[0], cs[0], ccs[0]);
+  SINCOS_TABLE_LOOKUP (u1, sn[1], ssn[1], cs[1], ccs[1]);
+#endif
+  cor = (ssn + s * ccs - sn * c) + cs * s;
+  res = sn + cor;
+  cor = (sn - res) + cor;
+  *corp = cor;
+  return res;
+}
+
 /* A more precise variant of DO_SIN.  EPS is the adjustment to the correction
    COR.  */
 static inline double
@@ -337,13 +440,13 @@  do_sincos_1 (double a, double da, double x, int4 n, bool shift_quadrant)
       if (xx < 0.01588)
 	{
 	  /* Taylor series.  */
-	  res = TAYLOR_SIN (xx, a, da, cor);
+	  res = taylor_sin_scalar (xx, a, da, &cor);
 	  cor = 1.02 * cor + __copysign (eps, cor);
 	  retval = (res == res + cor) ? res : sloww (a, da, x, shift_quadrant);
 	}
       else
 	{
-	  res = do_sin (a, da, &cor);
+	  res = do_sin_scalar (a, da, &cor);
 	  cor = 1.035 * cor + __copysign (eps, cor);
 	  retval = ((res == res + cor) ? __copysign (res, a)
 		    : sloww1 (a, da, x, shift_quadrant));
@@ -352,7 +455,7 @@  do_sincos_1 (double a, double da, double x, int4 n, bool shift_quadrant)
 
     case 1:
     case 3:
-      res = do_cos (a, da, &cor);
+      res = do_cos_scalar (a, da, &cor);
       cor = 1.025 * cor + __copysign (eps, cor);
       retval = ((res == res + cor) ? ((n & 2) ? -res : res)
 		: sloww2 (a, da, x, n));
@@ -411,13 +514,13 @@  do_sincos_2 (double a, double da, double x, int4 n, bool shift_quadrant)
       if (xx < 0.01588)
 	{
 	  /* Taylor series.  */
-	  res = TAYLOR_SIN (xx, a, da, cor);
+	  res = taylor_sin_scalar (xx, a, da, &cor);
 	  cor = 1.02 * cor + __copysign (eps, cor);
 	  retval = (res == res + cor) ? res : bsloww (a, da, x, n);
 	}
       else
 	{
-	  res = do_sin (a, da, &cor);
+	  res = do_sin_scalar (a, da, &cor);
 	  cor = 1.035 * cor + __copysign (eps, cor);
 	  retval = ((res == res + cor) ? __copysign (res, a)
 		    : bsloww1 (a, da, x, n));
@@ -426,7 +529,7 @@  do_sincos_2 (double a, double da, double x, int4 n, bool shift_quadrant)
 
     case 1:
     case 3:
-      res = do_cos (a, da, &cor);
+      res = do_cos_scalar (a, da, &cor);
       cor = 1.025 * cor + __copysign (eps, cor);
       retval = ((res == res + cor) ? ((n & 2) ? -res : res)
 		: bsloww2 (a, da, x, n));
@@ -436,33 +539,84 @@  do_sincos_2 (double a, double da, double x, int4 n, bool shift_quadrant)
   return retval;
 }
 
+static double
+__always_inline
+do_sincos_tail (double x, int4 k, int4 kl, bool shift_quadrant)
+{
+  int4 n;
+  double retval;
+  if (k < 0x419921FB)
+    {
+      double a, da;
+      n = reduce_sincos_1 (x, &a, &da);
+      retval = do_sincos_1 (a, da, x, n, shift_quadrant);
+    }
+  else if (k < 0x42F00000)
+    {
+      double a, da;
+      n = reduce_sincos_2 (x, &a, &da);
+      retval = do_sincos_2 (a, da, x, n, shift_quadrant);
+    }
+  else if (k < 0x7ff00000)
+    {
+      retval = reduce_and_compute (x, shift_quadrant);
+    }
+  else
+    {
+      if (k == 0x7ff00000 && kl == 0)
+	__set_errno (EDOM);
+      retval = x / x;
+    }
+  return retval;
+}
+
 /*******************************************************************/
 /* An ultimate sin routine. Given an IEEE double machine number x   */
 /* it computes the correctly rounded (to nearest) value of sin(x)  */
 /*******************************************************************/
 #ifdef IN_SINCOS
-static double
+static OP_TYPE
 #else
-double
+OP_TYPE
 SECTION
 #endif
-__sin (double x)
+__sin (OP_TYPE x)
 {
-  double xx, res, t, cor;
-  mynumber u;
-  int4 k, m;
-  double retval = 0;
+  OP_TYPE xx, res, res2, t, cor;
+  OP_TYPE retval;
 
 #ifndef IN_SINCOS
   SET_RESTORE_ROUND_53BIT (FE_TONEAREST);
 #endif
 
+#if VEC_SIZE == 1
+  mynumber u;
+  int4 m, k;
+  OP_TYPE zero = 0.0;
   u.x = x;
   m = u.i[HIGH_HALF];
   k = 0x7fffffff & m;		/* no sign           */
+#else /* VEC_SIZE == 2 */
+  mynumber u0, u1;
+  int4 m0, m1, k0, k1, k;
+  OP_TYPE zero = {0.0, 0.0};
+  u0.x = x[0];
+  u1.x = x[1];
+  m0 = u0.i[HIGH_HALF];
+  m1 = u1.i[HIGH_HALF];
+  k0 = 0x7fffffff & m0;
+  k1 = 0x7fffffff & m1;
+  k = k0 > k1 ? k0 : k1;
+#endif
+
   if (k < 0x3e500000)		/* if x->0 =>sin(x)=x */
     {
+#if VEC_SIZE == 1
       math_check_force_underflow (x);
+#else
+      math_check_force_underflow (x[0]);
+      math_check_force_underflow (x[1]);
+#endif
       retval = x;
     }
  /*---------------------------- 2^-26 < |x|< 0.25 ----------------------*/
@@ -473,54 +627,64 @@  __sin (double x)
       t = POLYNOMIAL (xx) * (xx * x);
       res = x + t;
       cor = (x - res) + t;
-      retval = (res == res + 1.07 * cor) ? res : slow (x);
+      res2 = res + 1.07 * cor;
+#if VEC_SIZE == 1
+      retval = (res == res2) ? res : slow (x);
+#else
+      retval[0] = (res[0] == res2[0]) ? res[0] : slow (x[0]);
+      retval[1] = (res[1] == res2[1]) ? res[1] : slow (x[1]);
+#endif
     }				/*  else  if (k < 0x3fd00000)    */
+
 /*---------------------------- 0.25<|x|< 0.855469---------------------- */
   else if (k < 0x3feb6000)
     {
-      res = do_sin (x, 0, &cor);
-      retval = (res == res + 1.096 * cor) ? res : slow1 (x);
+      res = do_sin_vector (x, zero, &cor);
+      res2 = res + 1.096 * cor;
+#if VEC_SIZE == 1
+      retval = (res == res2) ? res : slow1 (x);
       retval = __copysign (retval, x);
+#else
+      retval[0] = (res[0] == res2[0]) ? res[0] : slow1 (x[0]);
+      retval[1] = (res[1] == res2[1]) ? res[1] : slow1 (x[1]);
+      retval[0] = __copysign (retval[0], x[0]);
+      retval[1] = __copysign (retval[1], x[1]);
+#endif
     }				/*   else  if (k < 0x3feb6000)    */
 
 /*----------------------- 0.855469  <|x|<2.426265  ----------------------*/
   else if (k < 0x400368fd)
     {
-
+#if VEC_SIZE == 1
       t = hp0 - fabs (x);
-      res = do_cos (t, hp1, &cor);
-      retval = (res == res + 1.020 * cor) ? res : slow2 (x);
+      res = do_cos_vector (t, hp1, &cor);
+      res2 = res + 1.020 * cor;
+      retval = (res == res2) ? res : slow2 (x);
       retval = __copysign (retval, x);
+#else
+      OP_TYPE hp1_vec;
+      t[0] = hp0 - fabs (x[0]);
+      t[1] = hp0 - fabs (x[1]);
+      hp1_vec[0] = hp1_vec[1] = hp1;
+      res = do_cos_vector (t, hp1_vec, &cor);
+      res2 = res + 1.020 * cor;
+      retval[0] = (res[0] == res2[0]) ? res[0] : slow2 (x[0]);
+      retval[1] = (res[1] == res2[1]) ? res[1] : slow2 (x[1]);
+      retval[0] = __copysign (retval[0], x[0]);
+      retval[1] = __copysign (retval[1], x[1]);
+#endif
     }				/*   else  if (k < 0x400368fd)    */
 
 #ifndef IN_SINCOS
 /*-------------------------- 2.426265<|x|< 105414350 ----------------------*/
-  else if (k < 0x419921FB)
-    {
-      double a, da;
-      int4 n = reduce_sincos_1 (x, &a, &da);
-      retval = do_sincos_1 (a, da, x, n, false);
-    }				/*   else  if (k <  0x419921FB )    */
-
-/*---------------------105414350 <|x|< 281474976710656 --------------------*/
-  else if (k < 0x42F00000)
-    {
-      double a, da;
-
-      int4 n = reduce_sincos_2 (x, &a, &da);
-      retval = do_sincos_2 (a, da, x, n, false);
-    }				/*   else  if (k <  0x42F00000 )   */
-
-/* -----------------281474976710656 <|x| <2^1024----------------------------*/
-  else if (k < 0x7ff00000)
-    retval = reduce_and_compute (x, false);
-
-/*--------------------- |x| > 2^1024 ----------------------------------*/
   else
     {
-      if (k == 0x7ff00000 && u.i[LOW_HALF] == 0)
-	__set_errno (EDOM);
-      retval = x / x;
+#if VEC_SIZE == 1
+      retval = do_sincos_tail (x, k, u.i[LOW_HALF], false);
+#else
+      retval[0] = do_sincos_tail (x[0], k0, u0.i[LOW_HALF], false);
+      retval[1] = do_sincos_tail (x[1], k1, u1.i[LOW_HALF], false);
+#endif
     }
 #endif
 
@@ -534,85 +698,122 @@  __sin (double x)
 /*******************************************************************/
 
 #ifdef IN_SINCOS
-static double
+static OP_TYPE
 #else
-double
+OP_TYPE
 SECTION
 #endif
-__cos (double x)
+__cos (OP_TYPE x)
 {
-  double y, xx, res, cor, a, da;
-  mynumber u;
-  int4 k, m;
-
-  double retval = 0;
+  OP_TYPE y, xx, res, res2, a, da, cor;
+  OP_TYPE retval;
 
 #ifndef IN_SINCOS
   SET_RESTORE_ROUND_53BIT (FE_TONEAREST);
 #endif
 
+#if VEC_SIZE == 1
+  mynumber u;
+  int4 m, k;
+  OP_TYPE zero = 0.0;
+  OP_TYPE one = 1.0;
   u.x = x;
   m = u.i[HIGH_HALF];
-  k = 0x7fffffff & m;
+  k = 0x7fffffff & m;           /* no sign           */
+#else /* VEC_SIZE == 2 */
+  mynumber u0, u1;
+  int4 m0, m1, k0, k1, k;
+  OP_TYPE zero = {0.0, 0.0};
+  OP_TYPE one = {1.0, 1.0};
+  u0.x = x[0];
+  u1.x = x[1];
+  m0 = u0.i[HIGH_HALF];
+  m1 = u1.i[HIGH_HALF];
+  k0 = 0x7fffffff & m0;
+  k1 = 0x7fffffff & m1;
+  k = k0 > k1 ? k0 : k1;
+#endif
 
   /* |x|<2^-27 => cos(x)=1 */
   if (k < 0x3e400000)
-    retval = 1.0;
+      retval = one;
 
   else if (k < 0x3feb6000)
     {				/* 2^-27 < |x| < 0.855469 */
-      res = do_cos (x, 0, &cor);
-      retval = (res == res + 1.020 * cor) ? res : cslow2 (x);
+      res = do_cos_vector(x, zero, &cor);
+      res2 = res + 1.020 * cor;
+#if VEC_SIZE == 1
+      retval = (res == res2) ? res : cslow2 (x);
+#else
+      retval[0] = (res[0] == res2[0]) ? res[0] : cslow2 (x[0]);
+      retval[1] = (res[1] == res2[1]) ? res[1] : cslow2 (x[1]);
+#endif
     }				/*   else  if (k < 0x3feb6000)    */
 
   else if (k < 0x400368fd)
     { /* 0.855469  <|x|<2.426265  */ ;
+#if VEC_SIZE == 1
       y = hp0 - fabs (x);
       a = y + hp1;
       da = (y - a) + hp1;
       xx = a * a;
       if (xx < 0.01588)
 	{
-	  res = TAYLOR_SIN (xx, a, da, cor);
+	  res = taylor_sin_scalar (xx, a, da, &cor);
 	  cor = 1.02 * cor + __copysign (1.0e-31, cor);
 	  retval = (res == res + cor) ? res : sloww (a, da, x, true);
 	}
       else
 	{
-	  res = do_sin (a, da, &cor);
+	  res = do_sin_scalar (a, da, &cor);
 	  cor = 1.035 * cor + __copysign (1.0e-31, cor);
 	  retval = ((res == res + cor) ? __copysign (res, a)
 		    : sloww1 (a, da, x, true));
 	}
-
+#else
+      OP_TYPE hp1_vec, cors;
+      hp1_vec[0] = hp1_vec[1] = hp1;
+      y[0] = hp0 - fabs (x[0]);
+      y[1] = hp0 - fabs (x[1]);
+      a = y + hp1_vec;
+      da = (y - a) + hp1_vec;
+      xx = a * a;
+      if (xx[0] < 0.01588 && x[1] < 0.01588)
+	{
+	  res = taylor_sin_vector (xx, a, da, &cor);
+          cors[0] = __copysign (1.0e-31, cor[0]);
+          cors[1] = __copysign (1.0e-31, cor[1]);
+	  cor = 1.02 * cor + cors;
+          res2 = res + cor;
+	  retval[0] = (res[0] == res2[0]) ? res[0]
+					  : sloww (a[0], da[0], x[0], true);
+	  retval[1] = (res[1] == res2[1]) ? res[1]
+					  : sloww (a[1], da[1], x[1], true);
+	}
+      else
+	{
+	  res = do_sin_vector (a, da, &cor);
+	  cors[0] = __copysign (1.0e-31, cor[0]);
+	  cors[1] = __copysign (1.0e-31, cor[1]);
+	  cor = 1.035 * cor + cors;
+	  res2 = res + cor;
+          retval[0] = __copysign (res[0], a[0]);
+          retval[1] = __copysign (res[1], a[1]);
+	  retval[0] = (res[0] == res2[0]) ? retval[0] : sloww1 (a[0], da[0], x[0], true);
+	  retval[1] = (res[1] == res2[1]) ? retval[1] : sloww1 (a[1], da[1], x[1], true);
+	}
+#endif
     }				/*   else  if (k < 0x400368fd)    */
 
-
 #ifndef IN_SINCOS
-  else if (k < 0x419921FB)
-    {				/* 2.426265<|x|< 105414350 */
-      double a, da;
-      int4 n = reduce_sincos_1 (x, &a, &da);
-      retval = do_sincos_1 (a, da, x, n, true);
-    }				/*   else  if (k <  0x419921FB )    */
-
-  else if (k < 0x42F00000)
-    {
-      double a, da;
-
-      int4 n = reduce_sincos_2 (x, &a, &da);
-      retval = do_sincos_2 (a, da, x, n, true);
-    }				/*   else  if (k <  0x42F00000 )    */
-
-  /* 281474976710656 <|x| <2^1024 */
-  else if (k < 0x7ff00000)
-    retval = reduce_and_compute (x, true);
-
   else
     {
-      if (k == 0x7ff00000 && u.i[LOW_HALF] == 0)
-	__set_errno (EDOM);
-      retval = x / x;		/* |x| > 2^1024 */
+#if VEC_SIZE == 1
+      retval = do_sincos_tail (x, k, u.i[LOW_HALF], true);
+#else
+      retval[0] = do_sincos_tail (x[0], k0, u0.i[LOW_HALF], true);
+      retval[1] = do_sincos_tail (x[1], k1, u1.i[LOW_HALF], true);
+#endif
     }
 #endif