Update copyright dates with scripts/update-copyrights.

[thirdparty/glibc.git] / sysdeps / ieee754 / ldbl-128ibm / k_sinl.c

diff --git a/sysdeps/ieee754/ldbl-128ibm/k_sinl.c b/sysdeps/ieee754/ldbl-128ibm/k_sinl.c
index c4dce26f9a86964ba864dd78df8f66e024f9fe4e..44da02b0f3f0317903e50fd04e56d44e83d64fbc 100644 (file)
--- a/sysdeps/ieee754/ldbl-128ibm/k_sinl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/k_sinl.c
@@ -1,5 +1,5 @@
 /* Quad-precision floating point sine on <-pi/4,pi/4>.
-   Copyright (C) 1999,2004,2006 Free Software Foundation, Inc.
+   Copyright (C) 1999-2016 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
    Contributed by Jakub Jelinek <jj@ultra.linux.cz>
 
@@ -14,12 +14,12 @@
    Lesser General Public License for more details.
 
    You should have received a copy of the GNU Lesser General Public
-   License along with the GNU C Library; if not, write to the Free
-   Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA
-   02111-1307 USA.  */
+   License along with the GNU C Library; if not, see
+   <http://www.gnu.org/licenses/>.  */
 
-#include "math.h"
-#include "math_private.h"
+#include <float.h>
+#include <math.h>
+#include <math_private.h>
 
 static const long double c[] = {
 #define ONE c[0]
@@ -83,7 +83,10 @@ __kernel_sinl(long double x, long double y, int iy)
   long double h, l, z, sin_l, cos_l_m1;
   int64_t ix;
   u_int32_t tix, hix, index;
-  GET_LDOUBLE_MSW64 (ix, x);
+  double xhi, hhi;
+
+  xhi = ldbl_high (x);
+  EXTRACT_WORDS64 (ix, xhi);
   tix = ((u_int64_t)ix) >> 32;
   tix &= ~0x80000000;                  /* tix = |x|'s high 32 bits */
   if (tix < 0x3fc30000)                        /* |x| < 0.1484375 */
@@ -91,7 +94,10 @@ __kernel_sinl(long double x, long double y, int iy)
       /* Argument is small enough to approximate it by a Chebyshev
         polynomial of degree 17.  */
       if (tix < 0x3c600000)            /* |x| < 2^-57 */
-       if (!((int)x)) return x;        /* generate inexact */
+       {
+         math_check_force_underflow (x);
+         if (!((int)x)) return x;      /* generate inexact */
+       }
       z = x * x;
       return x + (x * (z*(SIN1+z*(SIN2+z*(SIN3+z*(SIN4+
                       z*(SIN5+z*(SIN6+z*(SIN7+z*SIN8)))))))));
@@ -104,6 +110,24 @@ __kernel_sinl(long double x, long double y, int iy)
         pre-computed tables,  compute cosl(l) and sinl(l) using a
         Chebyshev polynomial of degree 10(11) and compute
         sinl(h+l) = sinl(h)cosl(l) + cosl(h)sinl(l).  */
+      int six = tix;
+      tix = ((six - 0x3ff00000) >> 4) + 0x3fff0000;
+      index = 0x3ffe - (tix >> 16);
+      hix = (tix + (0x200 << index)) & (0xfffffc00 << index);
+      x = fabsl (x);
+      switch (index)
+       {
+       case 0: index = ((45 << 10) + hix - 0x3ffe0000) >> 8; break;
+       case 1: index = ((13 << 11) + hix - 0x3ffd0000) >> 9; break;
+       default:
+       case 2: index = (hix - 0x3ffc3000) >> 10; break;
+       }
+      hix = (hix << 4) & 0x3fffffff;
+/*
+    The following should work for double but generates the wrong index.
+    For now the code above converts double to ieee extended to compute
+    the index back to double for the h value.
+
       index = 0x3fe - (tix >> 20);
       hix = (tix + (0x2000 << index)) & (0xffffc000 << index);
       x = fabsl (x);
@@ -114,10 +138,11 @@ __kernel_sinl(long double x, long double y, int iy)
        default:
        case 2: index = (hix - 0x3fc30000) >> 14; break;
        }
-
-      SET_LDOUBLE_WORDS64(h, ((u_int64_t)hix) << 32, 0);
+*/
+      INSERT_WORDS64 (hhi, ((uint64_t)hix) << 32);
+      h = hhi;
       if (iy)
-       l = y - (h - x);
+       l = (ix < 0 ? -y : y) - (h - x);
       else
        l = x - h;
       z = l * l;