[newlib-cygwin] libm: Fixing overflow handling issue for scalbnf and scalbn

[Corinna Vinschen <corinna@sourceware.org>]

https://sourceware.org/git/gitweb.cgi?p=newlib-cygwin.git;h=ca7b4bd23676a9cb6d82974c2ed452e08abb7ecf

commit ca7b4bd23676a9cb6d82974c2ed452e08abb7ecf
Author: Kito Cheng <kito.cheng@sifive.com>
Date:   Mon Jul 19 16:50:22 2021 +0800

    libm: Fixing overflow handling issue for scalbnf and scalbn
    
    cc Aldy Hernandez <aldyh@redhat.com> and Andrew MacLeod <amacleod@redhat.com>,
    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`.

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