[PATCH] libm: Fixing overflow handling issue for scalbnf and scalbn

Kito Cheng <[email protected]>
Newsgroups gmane.comp.lib.newlib
Message-ID <[email protected]>
cc Aldy Hernandez <[email protected]> and Andrew MacLeod <[email protected]>,
they are author of new VRP analysis for GCC, just to make sure I didn't
mis-understanding or mis-interpreting anything on GCC site.

GCC 11 have better value range analysis, that give GCC more confidence
to perform more aggressive optimization, but it cause scalbn/scalbnf get
wrong result.

Using scalbn to demostrate what happened on GCC 11, see comments with VRP
prefix:

```c
double scalbn (double x, int n)
{
	/* VRP RESULT: n = [-INF, +INF] */
        __int32_t  k,hx,lx;
        ...
        k = (hx&0x7ff00000)>>20;
	/* VRP RESULT: k = [0, 2047] */
        if (k==0) {
	    /* VRP RESULT: k = 0 */
	    ...
	    k = ((hx&0x7ff00000)>>20) - 54;
            if (n< -50000) return tiny*x;       /*underflow*/
	    /* VRP RESULT: k = -54 */
	}
	/* VRP RESULT: k = [-54, 2047] */
        if (k==0x7ff) return x+x;               /* NaN or Inf */
	/* VRP RESULT: k = [-54, 2046] */
        k = k+n;
        if (k > 0x7fe) return huge*copysign(huge,x); /* overflow  */
	/* VRP RESULT: k = [-INF, 2046] */
	/* VRP RESULT: n = [-INF, 2100],
	   because k + n <= 0x7fe is false, so:
	   1. -INF < [-54, 2046] + n <= 0x7fe(2046) < INF
	   2. -INF < [-54, 2046] + n <= 2046 < INF
	   3. -INF < n <= 2046 - [-54, 2046] < INF
	   4. -INF < n <= [0, 2100] < INF
	   5. n = [-INF, 2100] */
        if (k > 0)                              /* normal result */
            {SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x;}
        if (k <= -54) {
	    /* VRP OPT: Evaluate n > 50000 as true...*/
            if (n > 50000)      /* in case integer overflow in n+k */
                return huge*copysign(huge,x);   /*overflow*/
            else return tiny*copysign(tiny,x);  /*underflow*/
	}
        k += 54;                                /* subnormal result */
        SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20));
        return x*twom54;
}
```

However give the input n = INT32_MAX, k = k+n will overflow, and then we
expect got `huge*copysign(huge,x)`, but new VRP optimization think
`n > 50000` is never be true, so optimize that into `tiny*copysign(tiny,x)`.

so the solution here is to moving the overflow handle logic before `k = k + n`.
---
 newlib/libm/common/s_scalbn.c  | 9 ++++-----
 newlib/libm/common/sf_scalbn.c | 9 ++++-----
 2 files changed, 8 insertions(+), 10 deletions(-)

diff --git a/newlib/libm/common/s_scalbn.c b/newlib/libm/common/s_scalbn.c
index e9dbb16ac..72da0e4a1 100644
--- a/newlib/libm/common/s_scalbn.c
+++ b/newlib/libm/common/s_scalbn.c
@@ -93,15 +93,14 @@ tiny   = 1.0e-300;
             if (n< -50000) return tiny*x; 	/*underflow*/
 	    }
         if (k==0x7ff) return x+x;		/* NaN or Inf */
+        if (n > 50000) 	/* in case integer overflow in n+k */
+	    return huge*copysign(huge,x);	/*overflow*/
         k = k+n; 
         if (k >  0x7fe) return huge*copysign(huge,x); /* overflow  */
         if (k > 0) 				/* normal result */
 	    {SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20)); return x;}
-        if (k <= -54) {
-            if (n > 50000) 	/* in case integer overflow in n+k */
-		return huge*copysign(huge,x);	/*overflow*/
-	    else return tiny*copysign(tiny,x); 	/*underflow*/
-      }
+        if (k <= -54)
+	    return tiny*copysign(tiny,x); 	/*underflow*/
         k += 54;				/* subnormal result */
 	SET_HIGH_WORD(x,(hx&0x800fffff)|(k<<20));
         return x*twom54;
diff --git a/newlib/libm/common/sf_scalbn.c b/newlib/libm/common/sf_scalbn.c
index 700060010..747a421ef 100644
--- a/newlib/libm/common/sf_scalbn.c
+++ b/newlib/libm/common/sf_scalbn.c
@@ -56,15 +56,14 @@ tiny   = 1.0e-30;
 	    k = ((ix&0x7f800000)>>23) - 25; 
             if (n< -50000) return tiny*x; 	/*underflow*/
         }
+        if (n > OVERFLOW_INT) 	/* in case integer overflow in n+k */
+	    return huge*copysignf(huge,x);	/*overflow*/
         k = k+n; 
         if (k > FLT_LARGEST_EXP) return huge*copysignf(huge,x); /* overflow  */
         if (k > 0) 				/* normal result */
 	    {SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23)); return x;}
-        if (k < FLT_SMALLEST_EXP) {
-            if (n > OVERFLOW_INT) 	/* in case integer overflow in n+k */
-		return huge*copysignf(huge,x);	/*overflow*/
-	    else return tiny*copysignf(tiny,x);	/*underflow*/
-        }
+        if (k < FLT_SMALLEST_EXP)
+	    return tiny*copysignf(tiny,x);	/*underflow*/
         k += 25;				/* subnormal result */
 	SET_FLOAT_WORD(x,(ix&0x807fffff)|(k<<23));
         return x*twom25;
-- 
2.31.1
lmpx.com only provides a reader for public news (NNTP) servers. It is not affiliated with the servers or forums shown here and is not responsible for the content of articles, which is written by their respective authors.