# commit 4cf69995e26e16005d4e3843ad4d18c75cf21a04

# Author: Alan Modra <amodra@gmail.com>

# Date:   Sat Aug 17 18:19:44 2013 +0930

# 

#     Fix for [BZ #15680] IBM long double inaccuracy

#     http://sourceware.org/ml/libc-alpha/2013-06/msg00919.html

#     

#     I discovered a number of places where denormals and other corner cases

#     were being handled wrongly.

#     

#     - printf_fphex.c: Testing for the low double exponent being zero is

#     unnecessary.  If the difference in exponents is less than 53 then the

#     high double exponent must be nearing the low end of its range, and the

#     low double exponent hit rock bottom.

#     

#     - ldbl2mpn.c: A denormal (ie. exponent of zero) value is treated as

#     if the exponent was one, so shift mantissa left by one.  Code handling

#     normalisation of the low double mantissa lacked a test for shift count

#     greater than bits in type being shifted, and lacked anything to handle

#     the case where the difference in exponents is less than 53 as in

#     printf_fphex.c.

#     

#     - math_ldbl.h (ldbl_extract_mantissa): Same as above, but worse, with

#     code testing for exponent > 1 for some reason, probably a typo for >= 1.

#     

#     - math_ldbl.h (ldbl_insert_mantissa): Round the high double as per

#     mpn2ldbl.c (hi is odd or explicit mantissas non-zero) so that the

#     number we return won't change when applying ldbl_canonicalize().

#     Add missing overflow checks and normalisation of high mantissa.

#     Correct misleading comment: "The hidden bit of the lo mantissa is

#     zero" is not always true as can be seen from the code rounding the hi

#     mantissa.  Also by inspection, lzcount can never be less than zero so

#     remove that test.  Lastly, masking bitfields to their widths can be

#     left to the compiler.

#     

#     - mpn2ldbl.c: The overflow checks here on rounding of high double were

#     just plain wrong.  Incrementing the exponent must be accompanied by a

#     shift right of the mantissa to keep the value unchanged.  Above notes

#     for ldbl_insert_mantissa are also relevant.

#     

#         [BZ #15680]

#         * sysdeps/ieee754/ldbl-128ibm/e_rem_pio2l.c: Comment fix.

#         * sysdeps/ieee754/ldbl-128ibm/printf_fphex.c

#         (PRINT_FPHEX_LONG_DOUBLE): Tidy code by moving -53 into ediff

#         calculation.  Remove unnecessary test for denormal exponent.

#         * sysdeps/ieee754/ldbl-128ibm/ldbl2mpn.c (__mpn_extract_long_double):

#         Correct handling of denormals.  Avoid undefined shift behaviour.

#         Correct normalisation of low mantissa when low double is denormal.

#         * sysdeps/ieee754/ldbl-128ibm/math_ldbl.h

#         (ldbl_extract_mantissa): Likewise.  Comment.  Use uint64_t* for hi64.

#         (ldbl_insert_mantissa): Make both hi64 and lo64 parms uint64_t.

#         Correct normalisation of low mantissa.  Test for overflow of high

#         mantissa and normalise.

#         (ldbl_nearbyint): Use more readable constant for two52.

#         * sysdeps/ieee754/ldbl-128ibm/mpn2ldbl.c

#         (__mpn_construct_long_double): Fix test for overflow of high

#         mantissa and correct normalisation.  Avoid undefined shift.

# 

diff -urN glibc-2.17-c758a686/sysdeps/ieee754/ldbl-128ibm/e_rem_pio2l.c glibc-2.17-c758a686/sysdeps/ieee754/ldbl-128ibm/e_rem_pio2l.c

--- glibc-2.17-c758a686/sysdeps/ieee754/ldbl-128ibm/e_rem_pio2l.c	2014-05-27 19:13:56.000000000 -0500

+++ glibc-2.17-c758a686/sysdeps/ieee754/ldbl-128ibm/e_rem_pio2l.c	2014-05-27 19:14:45.000000000 -0500

@@ -243,7 +243,7 @@

      We split the 113 bits of the mantissa into 5 24bit integers

      stored in a double array.  */

   /* Make the IBM extended format 105 bit mantissa look like the ieee854 112

-     bit mantissa so the next operatation will give the correct result.  */

+     bit mantissa so the next operation will give the correct result.  */

   ldbl_extract_mantissa (&ixd, &lxd, &exp, x);

   exp = exp - 23;

   /* This is faster than doing this in floating point, because we

diff -urN glibc-2.17-c758a686/sysdeps/ieee754/ldbl-128ibm/ldbl2mpn.c glibc-2.17-c758a686/sysdeps/ieee754/ldbl-128ibm/ldbl2mpn.c

--- glibc-2.17-c758a686/sysdeps/ieee754/ldbl-128ibm/ldbl2mpn.c	2014-05-27 19:13:56.000000000 -0500

+++ glibc-2.17-c758a686/sysdeps/ieee754/ldbl-128ibm/ldbl2mpn.c	2014-05-27 19:14:45.000000000 -0500

@@ -36,6 +36,7 @@

   union ibm_extended_long_double u;

   unsigned long long hi, lo;

   int ediff;

+

   u.ld = value;

   *is_neg = u.d[0].ieee.negative;

@@ -43,27 +44,36 @@

   lo = ((long long) u.d[1].ieee.mantissa0 << 32) | u.d[1].ieee.mantissa1;

   hi = ((long long) u.d[0].ieee.mantissa0 << 32) | u.d[0].ieee.mantissa1;

-  /* If the lower double is not a denomal or zero then set the hidden

+

+  /* If the lower double is not a denormal or zero then set the hidden

      53rd bit.  */

-  if (u.d[1].ieee.exponent > 0)

-    {

-      lo |= 1LL << 52;

+  if (u.d[1].ieee.exponent != 0)

+    lo |= 1ULL << 52;

+  else

+    lo = lo << 1;

-      /* The lower double is normalized separately from the upper.  We may

-	 need to adjust the lower manitissa to reflect this.  */

-      ediff = u.d[0].ieee.exponent - u.d[1].ieee.exponent;

-      if (ediff > 53)

-	lo = lo >> (ediff-53);

+  /* The lower double is normalized separately from the upper.  We may

+     need to adjust the lower manitissa to reflect this.  */

+  ediff = u.d[0].ieee.exponent - u.d[1].ieee.exponent - 53;

+  if (ediff > 0)

+    {

+      if (ediff < 64)

+	lo = lo >> ediff;

+      else

+	lo = 0;

     }

+  else if (ediff < 0)

+    lo = lo << -ediff;

+

   /* The high double may be rounded and the low double reflects the

      difference between the long double and the rounded high double

      value.  This is indicated by a differnce between the signs of the

      high and low doubles.  */

-  if ((u.d[0].ieee.negative != u.d[1].ieee.negative)

-      && ((u.d[1].ieee.exponent != 0) && (lo != 0L)))

+  if (u.d[0].ieee.negative != u.d[1].ieee.negative

+      && lo != 0)

     {

       lo = (1ULL << 53) - lo;

-      if (hi == 0LL)

+      if (hi == 0)

 	{

 	  /* we have a borrow from the hidden bit, so shift left 1.  */

 	  hi = 0x0ffffffffffffeLL | (lo >> 51);

diff -urN glibc-2.17-c758a686/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h glibc-2.17-c758a686/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h

--- glibc-2.17-c758a686/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h	2014-05-27 19:13:56.000000000 -0500

+++ glibc-2.17-c758a686/sysdeps/ieee754/ldbl-128ibm/math_ldbl.h	2014-05-27 19:51:13.000000000 -0500

@@ -13,77 +13,118 @@

      the number before the decimal point and the second implicit bit

      as bit 53 of the mantissa.  */

   uint64_t hi, lo;

-  int ediff;

   union ibm_extended_long_double u;

+

   u.ld = x;

   *exp = u.d[0].ieee.exponent - IEEE754_DOUBLE_BIAS;

   lo = ((uint64_t)u.d[1].ieee.mantissa0 << 32) | u.d[1].ieee.mantissa1;

   hi = ((uint64_t)u.d[0].ieee.mantissa0 << 32) | u.d[0].ieee.mantissa1;

-  /* If the lower double is not a denomal or zero then set the hidden

-     53rd bit.  */

-  if (u.d[1].ieee.exponent > 0x001)

-    {

-      lo |= (1ULL << 52);

-      lo = lo << 7; /* pre-shift lo to match ieee854.  */

-      /* The lower double is normalized separately from the upper.  We

-	 may need to adjust the lower manitissa to reflect this.  */

-      ediff = u.d[0].ieee.exponent - u.d[1].ieee.exponent;

-      if (ediff > 53)

-	lo = lo >> (ediff-53);

-      hi |= (1ULL << 52);

-    }

-  if ((u.d[0].ieee.negative != u.d[1].ieee.negative)

-      && ((u.d[1].ieee.exponent != 0) && (lo != 0LL)))

+  if (u.d[0].ieee.exponent != 0)

     {

-      hi--;

-      lo = (1ULL << 60) - lo;

-      if (hi < (1ULL << 52))

+      int ediff;

+

+      /* If not a denormal or zero then we have an implicit 53rd bit.  */

+      hi |= (uint64_t) 1 << 52;

+

+      if (u.d[1].ieee.exponent != 0)

+	lo |= (uint64_t) 1 << 52;

+      else

+	/* A denormal is to be interpreted as having a biased exponent

+	   of 1.  */

+	lo = lo << 1;

+

+      /* We are going to shift 4 bits out of hi later, because we only

+	 want 48 bits in *hi64.  That means we want 60 bits in lo, but

+	 we currently only have 53.  Shift the value up.  */

+      lo = lo << 7;

+

+      /* The lower double is normalized separately from the upper.

+	 We may need to adjust the lower mantissa to reflect this.

+	 The difference between the exponents can be larger than 53

+	 when the low double is much less than 1ULP of the upper

+	 (in which case there are significant bits, all 0's or all

+	 1's, between the two significands).  The difference between

+	 the exponents can be less than 53 when the upper double

+	 exponent is nearing its minimum value (in which case the low

+	 double is denormal ie. has an exponent of zero).  */

+      ediff = u.d[0].ieee.exponent - u.d[1].ieee.exponent - 53;

+      if (ediff > 0)

 	{

-	  /* we have a borrow from the hidden bit, so shift left 1.  */

-	  hi = (hi << 1) | (lo >> 59);

-	  lo = 0xfffffffffffffffLL & (lo << 1);

-	  *exp = *exp - 1;

+	  if (ediff < 64)

+	    lo = lo >> ediff;

+	  else

+	    lo = 0;

+	}

+      else if (ediff < 0)

+	lo = lo << -ediff;

+

+      if (u.d[0].ieee.negative != u.d[1].ieee.negative

+	  && lo != 0)

+	{

+	  hi--;

+	  lo = ((uint64_t) 1 << 60) - lo;

+	  if (hi < (uint64_t) 1 << 52)

+	    {

+	      /* We have a borrow from the hidden bit, so shift left 1.  */

+	      hi = (hi << 1) | (lo >> 59);

+	      lo = (((uint64_t) 1 << 60) - 1) & (lo << 1);

+	      *exp = *exp - 1;

+	    }

 	}

     }

+  else

+    /* If the larger magnitude double is denormal then the smaller

+       one must be zero.  */

+    hi = hi << 1;

+

   *lo64 = (hi << 60) | lo;

   *hi64 = hi >> 4;

 }

 static inline long double

-ldbl_insert_mantissa (int sign, int exp, int64_t hi64, u_int64_t lo64)

+ldbl_insert_mantissa (int sign, int exp, int64_t hi64, uint64_t lo64)

 {

   union ibm_extended_long_double u;

-  unsigned long hidden2, lzcount;

-  unsigned long long hi, lo;

+  int expnt2;

+  uint64_t hi, lo;

   u.d[0].ieee.negative = sign;

   u.d[1].ieee.negative = sign;

   u.d[0].ieee.exponent = exp + IEEE754_DOUBLE_BIAS;

-  u.d[1].ieee.exponent = exp-53 + IEEE754_DOUBLE_BIAS;

+  u.d[1].ieee.exponent = 0;

+  expnt2 = exp - 53 + IEEE754_DOUBLE_BIAS;

+ 

   /* Expect 113 bits (112 bits + hidden) right justified in two longs.

      The low order 53 bits (52 + hidden) go into the lower double */ 

-  lo = (lo64 >> 7)& ((1ULL << 53) - 1);

-  hidden2 = (lo64 >> 59) &  1ULL;

+  lo = (lo64 >> 7) & (((uint64_t) 1 << 53) - 1);

   /* The high order 53 bits (52 + hidden) go into the upper double */

-  hi = (lo64 >> 60) & ((1ULL << 11) - 1);

-  hi |= (hi64 << 4);

+  hi = lo64 >> 60;

+  hi |= hi64 << 4;

-  if (lo != 0LL)

+  if (lo != 0)

     {

-      /* hidden2 bit of low double controls rounding of the high double.

-	 If hidden2 is '1' then round up hi and adjust lo (2nd mantissa)

+      int lzcount;

+

+      /* 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 (hidden2)

+      if ((lo & ((uint64_t) 1 << 52)) != 0

+	  && ((hi & 1) != 0 || (lo & (((uint64_t) 1 << 52) - 1)) != 0))

 	{

 	  hi++;

+	  if ((hi & ((uint64_t) 1 << 53)) != 0)

+	    {

+	      hi = hi >> 1;

+	      u.d[0].ieee.exponent++;

+	    }

 	  u.d[1].ieee.negative = !sign;

-	  lo = (1ULL << 53) - lo;

+	  lo = ((uint64_t) 1 << 53) - lo;

 	}

-      /* The hidden bit of the lo mantissa is zero so we need to

-	 normalize the 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);

@@ -91,34 +132,30 @@

 	lzcount = __builtin_clzl ((long) (lo >> 32));

       else

 	lzcount = __builtin_clzl ((long) lo) + 32;

-      lzcount = lzcount - 11;

-      if (lzcount > 0)

+      lzcount = lzcount - (64 - 53);

+      lo <<= lzcount;

+      expnt2 -= lzcount;

+

+      if (expnt2 >= 1)

+	/* Not denormal.  */

+	u.d[1].ieee.exponent = expnt2;

+      else

 	{

-	  int expnt2 = u.d[1].ieee.exponent - lzcount;

-	  if (expnt2 >= 1)

-	    {

-	      /* Not denormal.  Normalize and set low exponent.  */

-	      lo = lo << lzcount;

-	      u.d[1].ieee.exponent = expnt2;

-	    }

+	  /* Is denormal.  Note that biased exponent of 0 is treated

+	     as if it was 1, hence the extra shift.  */

+	  if (expnt2 > -53)

+	    lo >>= 1 - expnt2;

 	  else

-	    {

-	      /* Is denormal.  */

-	      lo = lo << (lzcount + expnt2);

-	      u.d[1].ieee.exponent = 0;

-	    }

+	    lo = 0;

 	}

     }

   else

-    {

-      u.d[1].ieee.negative = 0;

-      u.d[1].ieee.exponent = 0;

-    }

+    u.d[1].ieee.negative = 0;

-  u.d[1].ieee.mantissa1 = lo & ((1ULL << 32) - 1);

-  u.d[1].ieee.mantissa0 = (lo >> 32) & ((1ULL << 20) - 1);

-  u.d[0].ieee.mantissa1 = hi & ((1ULL << 32) - 1);

-  u.d[0].ieee.mantissa0 = (hi >> 32) & ((1ULL << 20) - 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.ld;

 }

@@ -133,6 +170,10 @@

   return u.ld;

 }

+/* To suit our callers we return *hi64 and *lo64 as if they came from

+   an ieee854 112 bit mantissa, that is, 48 bits in *hi64 (plus one

+   implicit bit) and 64 bits in *lo64.  */

+

 static inline void

 default_ldbl_unpack (long double l, double *a, double *aa)

 {

@@ -162,13 +203,13 @@

   *aa = xl;

 }

-/* Simple inline nearbyint (double) function .

+/* Simple inline nearbyint (double) function.

    Only works in the default rounding mode

    but is useful in long double rounding functions.  */

 static inline double

 ldbl_nearbyint (double a)

 {

-  double two52 = 0x10000000000000LL;

+  double two52 = 0x1p52;

   if (__builtin_expect ((__builtin_fabs (a) < two52), 1))

     {

diff -urN glibc-2.17-c758a686/sysdeps/ieee754/ldbl-128ibm/mpn2ldbl.c glibc-2.17-c758a686/sysdeps/ieee754/ldbl-128ibm/mpn2ldbl.c

--- glibc-2.17-c758a686/sysdeps/ieee754/ldbl-128ibm/mpn2ldbl.c	2014-05-27 19:13:56.000000000 -0500

+++ glibc-2.17-c758a686/sysdeps/ieee754/ldbl-128ibm/mpn2ldbl.c	2014-05-27 19:14:45.000000000 -0500

@@ -70,9 +70,9 @@

       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.d[0].ieee.exponent)

 	{

@@ -98,29 +98,27 @@

 	}

     }

-  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;

+	      hi >>= 1;

 	      u.d[0].ieee.exponent++;

 	    }

 	  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);

@@ -128,24 +126,24 @@

 	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.d[1].ieee.exponent = exponent2;

-      else

+      else if (exponent2 > -53)

 	lo >>= 1 - exponent2;

+      else

+	lo = 0;

     }

   else

     u.d[1].ieee.negative = 0;

-  u.d[1].ieee.mantissa1 = lo & 0xffffffffLL;

-  u.d[1].ieee.mantissa0 = (lo >> 32) & 0xfffff;

-  u.d[0].ieee.mantissa1 = hi & 0xffffffffLL;

-  u.d[0].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.ld;

 }

diff -urN glibc-2.17-c758a686/sysdeps/ieee754/ldbl-128ibm/printf_fphex.c glibc-2.17-c758a686/sysdeps/ieee754/ldbl-128ibm/printf_fphex.c

--- glibc-2.17-c758a686/sysdeps/ieee754/ldbl-128ibm/printf_fphex.c	2014-05-27 19:13:56.000000000 -0500

+++ glibc-2.17-c758a686/sysdeps/ieee754/ldbl-128ibm/printf_fphex.c	2014-05-27 19:14:45.000000000 -0500

@@ -43,15 +43,15 @@

 	lo <<= 1;							      \

       /* The lower double is normalized separately from the upper.  We	      \

 	 may need to adjust the lower manitissa to reflect this.  */	      \

-      ediff = u.d[0].ieee.exponent - u.d[1].ieee.exponent;		      \

-      if (ediff > 53 + 63)						      \

+      ediff = u.d[0].ieee.exponent - u.d[1].ieee.exponent - 53;		      \

+      if (ediff > 63)							      \

 	lo = 0;								      \

-      else if (ediff > 53)						      \

-	lo = lo >> (ediff - 53);					      \

-      else if (u.d[1].ieee.exponent == 0 && ediff < 53)			      \

-	lo = lo << (53 - ediff);					      \

+      else if (ediff > 0)						      \

+	lo = lo >> ediff;						      \

+      else if (ediff < 0)						      \

+	lo = lo << -ediff;						      \

       if (u.d[0].ieee.negative != u.d[1].ieee.negative			      \

-	  && (u.d[1].ieee.exponent != 0 || lo != 0L))			      \

+	  && lo != 0)							      \

 	{								      \

 	  lo = (1ULL << 60) - lo;					      \

 	  if (hi == 0L)							      \