From 55976fad351dbdfce34b87a7054ef96859247e18 Mon Sep 17 00:00:00 2001
From: Ewan Crawford <ewan@codeplay.com>
Date: Wed, 28 Oct 2020 10:13:40 +0000
Subject: [PATCH] Permit half overflow within allowable ULP (#600)

* Permit half overflow within allowable ULP

Modify the algorithm for calculating half precision ULP error so
that it duplicates the behaviour of the single precision ULP algorithm,
in regards to allowing overflow within the defined ULP error.

In the case where the test value is infinity, but the reference is
finite, pretend the test value is 63336.0 and calculate the ULP error
against that.

Encountered this while testing half precision `hypot()` in PR !529,
for inputs `hypot(-48864.0, 43648.0)` which has reference
`65519.755799`. With RTE rounding this only just rounds to `65504` as half,
and returning INF is currently infinite ULP error. Using the leniency
introduced by this change however the error is `~0.5` within the `2` ULP
bounds defined by the spec.

* Run clang-format over changes

Code now conforms to style guidelines and allows `check-format.sh` to pass.
---
 test_common/harness/errorHelpers.cpp        | 27 +++++++-----
 test_conformance/math_brute_force/Utility.h |  1 -
 test_conformance/math_brute_force/main.cpp  | 47 ---------------------
 3 files changed, 17 insertions(+), 58 deletions(-)

diff --git a/test_common/harness/errorHelpers.cpp b/test_common/harness/errorHelpers.cpp
index 5fee9a2c..fd5a75e5 100644
--- a/test_common/harness/errorHelpers.cpp
+++ b/test_common/harness/errorHelpers.cpp
@@ -307,8 +307,6 @@ const char *GetQueuePropertyName(cl_command_queue_properties property)
 #define scalbnl(_a, _i )    ldexpl( _a, _i )
 #endif
 
-static float Ulp_Error_Half_Float( float test, double reference );
-
 // taken from math tests
 #define HALF_MIN_EXP    -13
 #define HALF_MANT_DIG    11
@@ -327,6 +325,23 @@ static float Ulp_Error_Half_Float( float test, double reference )
     // results.
 
     double testVal = test;
+
+    if (isinf(reference))
+    {
+        if (testVal == reference) return 0.0f;
+
+        return (float)(testVal - reference);
+    }
+
+    if (isinf(testVal))
+    {
+        // Allow overflow within the limit of the allowed ulp error. Towards
+        // that end we pretend the test value is actually 2**16, the next value
+        // that would appear in the number line if half had sufficient range.
+        testVal = copysign(65536.0, testVal);
+    }
+
+
     if( u.u & 0x000fffffffffffffULL )
     { // Non-power of two and NaN
         if( isnan( reference ) && isnan( test ) )
@@ -339,14 +354,6 @@ static float Ulp_Error_Half_Float( float test, double reference )
         return (float) scalbn( testVal - reference, ulp_exp );
     }
 
-    if( isinf( reference ) )
-    {
-        if( (double) test == reference )
-            return 0.0f;
-
-        return (float) (testVal - reference );
-    }
-
     // reference is a normal power of two or a zero
     int ulp_exp =  HALF_MANT_DIG - 1 - MAX( ilogb( reference) - 1, HALF_MIN_EXP-1 );
 
diff --git a/test_conformance/math_brute_force/Utility.h b/test_conformance/math_brute_force/Utility.h
index 7a4cac6c..31256358 100644
--- a/test_conformance/math_brute_force/Utility.h
+++ b/test_conformance/math_brute_force/Utility.h
@@ -95,7 +95,6 @@ extern cl_device_fp_config gDoubleCapabilities;
 
 float Abs_Error( float test, double reference );
 float Ulp_Error( float test, double reference );
-//float Ulp_Error_Half( float test, double reference );
 float Bruteforce_Ulp_Error_Double( double test, long double reference );
 
 uint64_t GetTime( void );
diff --git a/test_conformance/math_brute_force/main.cpp b/test_conformance/math_brute_force/main.cpp
index 142ac64c..8f2e0a0c 100644
--- a/test_conformance/math_brute_force/main.cpp
+++ b/test_conformance/math_brute_force/main.cpp
@@ -1739,53 +1739,6 @@ float Abs_Error( float test, double reference )
   return fabs((float)(reference-(double)test));
 }
 
-/*
-#define HALF_MIN_EXP    -13
-#define HALF_MANT_DIG    11
-float Ulp_Error_Half( float test, double reference )
-{
-    union{ double d; uint64_t u; }u;     u.d = reference;
-
-  // Note: This function presumes that someone has already tested whether the result is correctly,
-  // rounded before calling this function.  That test:
-  //
-  //    if( (float) reference == test )
-  //        return 0.0f;
-  //
-  // would ensure that cases like fabs(reference) > FLT_MAX are weeded out before we get here.
-  // Otherwise, we'll return inf ulp error here, for what are otherwise correctly rounded
-  // results.
-
-    double testVal = test;
-    if( u.u & 0x000fffffffffffffULL )
-    { // Non-power of two and NaN
-        if( isnan( reference ) && isnan( test ) )
-            return 0.0f;    // if we are expecting a NaN, any NaN is fine
-
-        // The unbiased exponent of the ulp unit place
-        int ulp_exp = HALF_MANT_DIG - 1 - MAX( ilogb( reference), HALF_MIN_EXP-1 );
-
-        // Scale the exponent of the error
-        return (float) scalbn( testVal - reference, ulp_exp );
-    }
-
-    if( isinf( reference ) )
-    {
-        if( (double) test == reference )
-            return 0.0f;
-
-        return (float) (testVal - reference );
-    }
-
-    // reference is a normal power of two or a zero
-    int ulp_exp =  HALF_MANT_DIG - 1 - MAX( ilogb( reference) - 1, HALF_MIN_EXP-1 );
-
-    // Scale the exponent of the error
-    return (float) scalbn( testVal - reference, ulp_exp );
-}
-*/
-
-
 #if defined( __APPLE__ )
     #include <mach/mach_time.h>
 #endif