Avoid some undefined behavior in test_bruteforce. (#2400)

* Ulp_Error*: ilogb(reference) - 1 may overflow if reference is zero.
* binary_i_double Test: DoubleFromUInt32's result is a cl_double and the
attempt is to store it as a cl_double, but p was defined as a pointer to
cl_ulong, resulting in an unintended implicit conversion that is not
valid for out-of-range doubles.
* exp2, tanpi: ensure early exit for NaN.
* shift_right_sticky_128: avoid out-of-range shift if shift value is
exactly 64.
* scalbn: e += n may overflow if n is large, move it after the check for
large n.

test_common/harness/errorHelpers.cpp

@@ -387,8 +387,7 @@ static float Ulp_Error_Half_Float(float test, double reference)
     }
 
     // reference is a normal power of two or a zero
-    int ulp_exp =
-        HALF_MANT_DIG - 1 - std::max(ilogb(reference) - 1, HALF_MIN_EXP - 1);
+    int ulp_exp = HALF_MANT_DIG - std::max(ilogb(reference), HALF_MIN_EXP);
 
     // Scale the exponent of the error
     return (float)scalbn(testVal - reference, ulp_exp);
@@ -469,8 +468,7 @@ float Ulp_Error(float test, double reference)
 
     // reference is a normal power of two or a zero
     // The unbiased exponent of the ulp unit place
-    int ulp_exp =
-        FLT_MANT_DIG - 1 - std::max(ilogb(reference) - 1, FLT_MIN_EXP - 1);
+    int ulp_exp = FLT_MANT_DIG - std::max(ilogb(reference), FLT_MIN_EXP);
 
     // Scale the exponent of the error
     return (float)scalbn(testVal - reference, ulp_exp);
@@ -553,8 +551,7 @@ float Ulp_Error_Double(double test, long double reference)
 
     // reference is a normal power of two or a zero
     // The unbiased exponent of the ulp unit place
-    int ulp_exp =
-        DBL_MANT_DIG - 1 - std::max(ilogbl(reference) - 1, DBL_MIN_EXP - 1);
+    int ulp_exp = DBL_MANT_DIG - std::max(ilogbl(reference), DBL_MIN_EXP);
 
     // Scale the exponent of the error
     float result = (float)scalbnl(testVal - reference, ulp_exp);

test_conformance/math_brute_force/binary_i_double.cpp

@@ -248,7 +248,7 @@ cl_int Test(cl_uint job_id, cl_uint thread_id, void *data)
     }
 
     // Init input array
-    cl_ulong *p = (cl_ulong *)gIn + thread_id * buffer_elements;
+    cl_double *p = (cl_double *)gIn + thread_id * buffer_elements;
     cl_int *p2 = (cl_int *)gIn2 + thread_id * buffer_elements;
     size_t idx = 0;
     int totalSpecialValueCount = specialValuesCount * specialValuesIntCount;
@@ -257,7 +257,6 @@ cl_int Test(cl_uint job_id, cl_uint thread_id, void *data)
     // Test edge cases
     if (job_id <= (cl_uint)lastSpecialJobIndex)
     {
-        cl_double *fp = (cl_double *)p;
         cl_int *ip2 = (cl_int *)p2;
         uint32_t x, y;
 
@@ -266,7 +265,7 @@ cl_int Test(cl_uint job_id, cl_uint thread_id, void *data)
 
         for (; idx < buffer_elements; idx++)
         {
-            fp[idx] = specialValues[x];
+            p[idx] = specialValues[x];
             ip2[idx] = specialValuesInt[y];
             if (++x >= specialValuesCount)
             {

test_conformance/math_brute_force/main.cpp

@@ -1330,8 +1330,7 @@ float Bruteforce_Ulp_Error_Double(double test, long double reference)
 
     // reference is a normal power of two or a zero
     // The unbiased exponent of the ulp unit place
-    int ulp_exp =
-        DBL_MANT_DIG - 1 - std::max(ilogbl(reference) - 1, DBL_MIN_EXP - 1);
+    int ulp_exp = DBL_MANT_DIG - std::max(ilogbl(reference), DBL_MIN_EXP);
 
     // allow correctly rounded results to pass through unmolested. (We might add
     // error to it below.) There is something of a performance optimization here

test_conformance/math_brute_force/reference_math.cpp

@@ -721,9 +721,9 @@ double reference_tanpi(double x)
     double z = reference_fabs(x);
 
     // if big and even  -- caution: only works if x only has single precision
-    if (z >= HEX_DBL(+, 1, 0, +, 24))
+    if (!(z < HEX_DBL(+, 1, 0, +, 24)))
     {
-        if (z == INFINITY) return x - x; // nan
+        if (!isfinite(z)) return x - x; // nan
 
         return reference_copysign(
             0.0, x); // tanpi ( n ) is copysign( 0.0, n)  for even integers n.
@@ -1223,6 +1223,8 @@ double reference_relaxed_exp2(double x) { return reference_exp2(x); }
 double reference_exp2(double x)
 { // Note: only suitable for verifying single precision. Doesn't have range of a
   // full double exp2 implementation.
+    if (isnan(x)) return x;
+
     if (x == 0.0) return 1.0;
 
     // separate x into fractional and integer parts
@@ -2781,7 +2783,7 @@ static inline void shift_right_sticky_128(cl_ulong *hi, cl_ulong *lo, int shift)
             sticky |= (0 != l);
             l = 0;
         }
-        else
+        else if (shift > 0)
         {
             sticky |= (0 != (l << (64 - shift)));
             l >>= shift;
@@ -3088,9 +3090,9 @@ long double reference_tanpil(long double x)
     long double z = reference_fabsl(x);
 
     // if big and even  -- caution: only works if x only has single precision
-    if (z >= HEX_LDBL(+, 1, 0, +, 53))
+    if (!(z < HEX_LDBL(+, 1, 0, +, 53)))
     {
-        if (z == INFINITY) return x - x; // nan
+        if (!isfinite(z)) return x - x; // nan
 
         return reference_copysignl(
             0.0L, x); // tanpi ( n ) is copysign( 0.0, n)  for even integers n.
@@ -5027,8 +5029,9 @@ static double reference_scalbn(double x, int n)
         u.d -= 1.0;
         e = (int)((u.l & 0x7ff0000000000000LL) >> 52) - 1022;
     }
+    if (n >= 2098) return reference_copysign(INFINITY, x);
     e += n;
-    if (e >= 2047 || n >= 2098) return reference_copysign(INFINITY, x);
+    if (e >= 2047) return reference_copysign(INFINITY, x);
     if (e < -51 || n < -2097) return reference_copysign(0.0, x);
     if (e <= 0)
     {