From 9e329b544ac04b389e12984362e06b92e1506399 Mon Sep 17 00:00:00 2001
From: Jesse Huang via Newlib <newlib@sourceware.org>
Date: Thu, 27 Jul 2023 01:49:41 -0700
Subject: Fix rounding results in lrint() & llrint() when close to 0

soft-fp should round floating pointer numbers according to the current
rounding mode. However, in the current code of lrint() and llrint(),
there are if statements before the actual rounding computation

      if(j0 < -1)
        return 0;

Where j0 is the exponent of the floating point number.

It means any number having a exponent less than -1
(i.e. interval (-0.5, 0.5)) will be rounded to 0 regardeless of the
rounding mode.

The bug already fixed in glibc in 2006 by moving the check afterwards
the rounding computation, but still persists in newlib.

This patch fixed it in a similar way to glibc
Ref Commit in glibc: 6624dbc07b5a9fb316ed188ef01f65b8eea8b47c
---
 newlib/libm/common/s_llrint.c  | 31 +++++++++++--------------------
 newlib/libm/common/s_lrint.c   | 31 +++++++++++--------------------
 newlib/libm/common/sf_llrint.c |  6 ++----
 newlib/libm/common/sf_lrint.c  |  6 ++----
 4 files changed, 26 insertions(+), 48 deletions(-)

(limited to 'newlib/libm')

diff --git a/newlib/libm/common/s_llrint.c b/newlib/libm/common/s_llrint.c
index 72452db..f2c4876 100644
--- a/newlib/libm/common/s_llrint.c
+++ b/newlib/libm/common/s_llrint.c
@@ -66,26 +66,17 @@ long long int
   if(j0 < 20)
     {
       /* j0 in [-1023,19] */
-      if(j0 < -1)
-        return 0;
-      else
-        {
-          /* j0 in [0,19] */
-	  /* shift amt in [0,19] */
-          w = TWO52[sx] + x;
-          t = w - TWO52[sx];
-          GET_HIGH_WORD(i0, t);
-          /* Detect the all-zeros representation of plus and
-             minus zero, which fails the calculation below. */
-          if ((i0 & ~((__int32_t)1 << 31)) == 0)
-              return 0;
-          /* After round:  j0 in [0,20] */
-          j0 = ((i0 & 0x7ff00000) >> 20) - 1023;
-          i0 &= 0x000fffff;
-          i0 |= 0x00100000;
-	  /* shift amt in [20,0] */
-          result = i0 >> (20 - j0);
-        }
+      w = TWO52[sx] + x;
+      t = w - TWO52[sx];
+      GET_HIGH_WORD(i0, t);
+      /* Detect the all-zeros representation of plus and
+         minus zero, which fails the calculation below. */
+      if ((i0 & ~((__int32_t)1 << 31)) == 0)
+          return 0;
+      j0 = ((i0 & 0x7ff00000) >> 20) - 1023;
+      i0 &= 0x000fffff;
+      i0 |= 0x00100000;
+      result = (j0 < 0 ? 0 : i0 >> (20 - j0));
     }
   else if (j0 < (int)(8 * sizeof (long long int)) - 1)
     {
diff --git a/newlib/libm/common/s_lrint.c b/newlib/libm/common/s_lrint.c
index b37f50f..b37b93a 100644
--- a/newlib/libm/common/s_lrint.c
+++ b/newlib/libm/common/s_lrint.c
@@ -103,26 +103,17 @@ TWO52[2]={
   if(j0 < 20)
     {
       /* j0 in [-1023,19] */
-      if(j0 < -1)
-        return 0;
-      else
-        {
-          /* j0 in [0,19] */
-	  /* shift amt in [0,19] */
-          w = TWO52[sx] + x;
-          t = w - TWO52[sx];
-          GET_HIGH_WORD(i0, t);
-          /* Detect the all-zeros representation of plus and
-             minus zero, which fails the calculation below. */
-          if ((i0 & ~(1L << 31)) == 0)
-              return 0;
-          /* After round:  j0 in [0,20] */
-          j0 = ((i0 & 0x7ff00000) >> 20) - 1023;
-          i0 &= 0x000fffff;
-          i0 |= 0x00100000;
-	  /* shift amt in [20,0] */
-          result = i0 >> (20 - j0);
-        }
+      w = TWO52[sx] + x;
+      t = w - TWO52[sx];
+      GET_HIGH_WORD(i0, t);
+      /* Detect the all-zeros representation of plus and
+         minus zero, which fails the calculation below. */
+      if ((i0 & ~(1L << 31)) == 0)
+          return 0;
+      j0 = ((i0 & 0x7ff00000) >> 20) - 1023;
+      i0 &= 0x000fffff;
+      i0 |= 0x00100000;
+      result = (j0 < 0 ? 0 : i0 >> (20 - j0));
     }
   else if (j0 < (int)(8 * sizeof (long int)) - 1)
     {
diff --git a/newlib/libm/common/sf_llrint.c b/newlib/libm/common/sf_llrint.c
index 7558e89..905a5b2 100644
--- a/newlib/libm/common/sf_llrint.c
+++ b/newlib/libm/common/sf_llrint.c
@@ -60,9 +60,7 @@ TWO23[2]={
   
   if (j0 < (int)(sizeof (long long int) * 8) - 1)
     {
-      if (j0 < -1)
-        return 0;
-      else if (j0 >= 23)
+      if (j0 >= 23)
         result = (long long int) ((i0 & 0x7fffff) | 0x800000) << (j0 - 23);
       else
         {
@@ -76,7 +74,7 @@ TWO23[2]={
           j0 = ((i0 >> 23) & 0xff) - 0x7f;
           i0 &= 0x7fffff;
           i0 |= 0x800000;
-          result = i0 >> (23 - j0);
+          result = (j0 < 0 ? 0 : i0 >> (23 - j0));
         }
     }
   else
diff --git a/newlib/libm/common/sf_lrint.c b/newlib/libm/common/sf_lrint.c
index 3c58c5d..7fe47ae 100644
--- a/newlib/libm/common/sf_lrint.c
+++ b/newlib/libm/common/sf_lrint.c
@@ -60,9 +60,7 @@ TWO23[2]={
   
   if (j0 < (int)(sizeof (long int) * 8) - 1)
     {
-      if (j0 < -1)
-        return 0;
-      else if (j0 >= 23)
+      if (j0 >= 23)
         result = (long int) ((i0 & 0x7fffff) | 0x800000) << (j0 - 23);
       else
         {
@@ -76,7 +74,7 @@ TWO23[2]={
           j0 = ((i0 >> 23) & 0xff) - 0x7f;
           i0 &= 0x7fffff;
           i0 |= 0x800000;
-          result = i0 >> (23 - j0);
+          result = (j0 < 0 ? 0 : i0 >> (23 - j0));
         }
     }
   else
-- 
cgit v1.1