Update copyright dates with scripts/update-copyrights.

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

diff --git a/sysdeps/ieee754/ldbl-128ibm/mpn2ldbl.c b/sysdeps/ieee754/ldbl-128ibm/mpn2ldbl.c
index 21d1e62dabd64f9fd2781d17399cbae762c693a5..42f5e6a02d4f58341ec91dc7f08ee8718cda68c9 100644 (file)
--- a/sysdeps/ieee754/ldbl-128ibm/mpn2ldbl.c
+++ b/sysdeps/ieee754/ldbl-128ibm/mpn2ldbl.c
@@ -1,5 +1,4 @@
-/* Copyright (C) 1995, 1996, 1997, 1998, 1999, 2002, 2003, 2004, 2006, 2007
-       Free Software Foundation, Inc.
+/* Copyright (C) 1995-2016 Free Software Foundation, Inc.
    This file is part of the GNU C Library.
 
    The GNU C Library is free software; you can redistribute it and/or
@@ -13,16 +12,20 @@
    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 "gmp.h"
-#include "gmp-impl.h"
 #include <ieee754.h>
+#include <errno.h>
 #include <float.h>
 #include <math.h>
 
+/* Need to set this when including gmp headers after system headers.  */
+#define HAVE_ALLOCA 1
+
+#include "gmp.h"
+#include "gmp-impl.h"
+
 /* Convert a multi-precision integer of the needed number of bits (106
    for long double) and an integral power of two to a `long double' in
    IBM extended format.  */
@@ -35,11 +38,11 @@ __mpn_construct_long_double (mp_srcptr frac_ptr, int expt, int sign)
   unsigned long long hi, lo;
   int exponent2;
 
-  u.ieee.negative = sign;
-  u.ieee.negative2 = sign;
-  u.ieee.exponent = expt + IBM_EXTENDED_LONG_DOUBLE_BIAS;
-  u.ieee.exponent2 = 0;
-  exponent2 = expt - 53 + IBM_EXTENDED_LONG_DOUBLE_BIAS;
+  u.d[0].ieee.negative = sign;
+  u.d[1].ieee.negative = sign;
+  u.d[0].ieee.exponent = expt + IEEE754_DOUBLE_BIAS;
+  u.d[1].ieee.exponent = 0;
+  exponent2 = expt - 53 + IEEE754_DOUBLE_BIAS;
 
 #if BITS_PER_MP_LIMB == 32
   /* The low order 53 bits (52 + hidden) go into the lower double */
@@ -71,19 +74,19 @@ __mpn_construct_long_double (mp_srcptr frac_ptr, int expt, int sign)
       else
        lzcount = __builtin_clzl ((long) val) + 32;
       if (hi)
-       lzcount = lzcount - 11;
+       lzcount = lzcount - (64 - 53);
       else
-       lzcount = lzcount + 42;
+       lzcount = lzcount + 53 - (64 - 53);
 
-      if (lzcount > u.ieee.exponent)
+      if (lzcount > u.d[0].ieee.exponent)
        {
-         lzcount = u.ieee.exponent;
-         u.ieee.exponent = 0;
+         lzcount = u.d[0].ieee.exponent;
+         u.d[0].ieee.exponent = 0;
          exponent2 -= lzcount;
        }
       else
        {
-         u.ieee.exponent -= (lzcount - 1);
+         u.d[0].ieee.exponent -= (lzcount - 1);
          exponent2 -= (lzcount - 1);
        }
 
@@ -99,29 +102,35 @@ __mpn_construct_long_double (mp_srcptr frac_ptr, int expt, int sign)
        }
     }
 
-  if (lo != 0L)
+  if (lo != 0)
     {
-      /* hidden2 bit of low double controls rounding of the high double.
-        If hidden2 is '1' and either the explicit mantissa is non-zero
+      /* hidden bit of low double controls rounding of the high double.
+        If hidden is '1' and either the explicit mantissa is non-zero
         or hi is odd, then round up hi and adjust lo (2nd mantissa)
         plus change the sign of the low double to compensate.  */
       if ((lo & (1LL << 52)) != 0
-         && ((hi & 1) != 0 || (lo & ((1LL << 52) - 1))))
+         && ((hi & 1) != 0 || (lo & ((1LL << 52) - 1)) != 0))
        {
          hi++;
-         if ((hi & ((1LL << 52) - 1)) == 0)
+         if ((hi & (1LL << 53)) != 0)
            {
-             if ((hi & (1LL << 53)) != 0)
-               hi -= 1LL << 52;
-             u.ieee.exponent++;
+             hi >>= 1;
+             u.d[0].ieee.exponent++;
+             if (u.d[0].ieee.exponent == IEEE754_DOUBLE_BIAS + DBL_MAX_EXP)
+               {
+                 /* Overflow.  The appropriate overflowed result must
+                    be produced (if an infinity, that means the low
+                    part must be zero).  */
+                 __set_errno (ERANGE);
+                 return (sign ? -LDBL_MAX : LDBL_MAX) * LDBL_MAX;
+               }
            }
-         u.ieee.negative2 = !sign;
+         u.d[1].ieee.negative = !sign;
          lo = (1LL << 53) - lo;
        }
 
-      /* The hidden bit of the lo mantissa is zero so we need to normalize
-        it for the low double.  Shift it left until the hidden bit is '1'
-        then adjust the 2nd exponent accordingly.  */
+      /* Normalize the low double.  Shift the mantissa left until
+        the hidden bit is '1' and adjust the exponent accordingly.  */
 
       if (sizeof (lo) == sizeof (long))
        lzcount = __builtin_clzl (lo);
@@ -129,24 +138,24 @@ __mpn_construct_long_double (mp_srcptr frac_ptr, int expt, int sign)
        lzcount = __builtin_clzl ((long) (lo >> 32));
       else
        lzcount = __builtin_clzl ((long) lo) + 32;
-      lzcount = lzcount - 11;
-      if (lzcount > 0)
-       {
-         lo = lo << lzcount;
-         exponent2 = exponent2 - lzcount;
-       }
+      lzcount = lzcount - (64 - 53);
+      lo <<= lzcount;
+      exponent2 -= lzcount;
+
       if (exponent2 > 0)
-       u.ieee.exponent2 = exponent2;
-      else
+       u.d[1].ieee.exponent = exponent2;
+      else if (exponent2 > -53)
        lo >>= 1 - exponent2;
+      else
+       lo = 0;
     }
   else
-    u.ieee.negative2 = 0;
+    u.d[1].ieee.negative = 0;
 
-  u.ieee.mantissa3 = lo & 0xffffffffLL;
-  u.ieee.mantissa2 = (lo >> 32) & 0xfffff;
-  u.ieee.mantissa1 = hi & 0xffffffffLL;
-  u.ieee.mantissa0 = (hi >> 32) & ((1LL << (LDBL_MANT_DIG - 86)) - 1);
+  u.d[1].ieee.mantissa1 = lo;
+  u.d[1].ieee.mantissa0 = lo >> 32;
+  u.d[0].ieee.mantissa1 = hi;
+  u.d[0].ieee.mantissa0 = hi >> 32;
 
-  return u.d;
+  return u.ld;
 }