Add ULP verification for fp32

Add a verifier to check two results are correct within a certain ULP
tolerance for IEEE-754 32-bit floating point values.

Add a test to check the ULP verifier is correct.

Signed-off-by: Jack Frankland <jack.frankland@arm.com>
Change-Id: Iaf43069f300999479d998e7837746b247ca5177e

4 files changed, 166 insertions, 1 deletions

diff --git a/reference_model/src/verify/verifiers.h b/reference_model/src/verify/verifiers.h
index 177eeaf..bdc8fe7 100644
--- a/reference_model/src/verify/verifiers.h
+++ b/reference_model/src/verify/verifiers.h
@@ -41,6 +41,15 @@ bool verifyDotProduct(const CTensor* ref,
 /// \return True if compliant else false
 bool verifyExact(const CTensor* referenceTensor, const CTensor* implementationTensor);
 
+/// \brief Perform ULP result verification
+///
+/// \param referenceTensor    Reference tensor
+/// \param implementationTensor    Implementation resulting tensor
+/// \param ulp    The ULP tolerence for the comparison of the two tensors
+///
+/// \return True if compliant else false
+bool verifyULP(const CTensor* referenceTensor, const CTensor* implementationTensor, uint64_t ulp);
+
 };    // namespace TosaReference
 
 #endif    // VERIFIERS_H_
diff --git a/reference_model/src/verify/verify_entry.cc b/reference_model/src/verify/verify_entry.cc
index 1c48354..1f7c680 100644
--- a/reference_model/src/verify/verify_entry.cc
+++ b/reference_model/src/verify/verify_entry.cc
@@ -34,6 +34,9 @@ bool verify(const CTensor* ref, const CTensor* refBnd, const CTensor* imp, const
         case VerifyMode::Exact: {
             return verifyExact(ref, imp);
         }
+        case VerifyMode::Ulp: {
+            return verifyULP(ref, imp, cfg.ulpInfo.ulp);
+        }
         default: {
             WARNING("unsupported verification mode.");
             break;
diff --git a/reference_model/src/verify/verify_ulp.cc b/reference_model/src/verify/verify_ulp.cc
new file mode 100644
index 0000000..622fba4
--- /dev/null
+++ b/reference_model/src/verify/verify_ulp.cc
@@ -0,0 +1,153 @@
+// Copyright (c) 2023, ARM Limited.
+//
+//    Licensed under the Apache License, Version 2.0 (the "License");
+//    you may not use this file except in compliance with the License.
+//    You may obtain a copy of the License at
+//
+//         http://www.apache.org/licenses/LICENSE-2.0
+//
+//    Unless required by applicable law or agreed to in writing, software
+//    distributed under the License is distributed on an "AS IS" BASIS,
+//    WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+//    See the License for the specific language governing permissions and
+//    limitations under the License.
+
+#include <cmath>
+#include <limits>
+#include <memory>
+#include <type_traits>
+
+#include "verifiers.h"
+
+namespace TosaReference
+{
+
+namespace
+{
+static_assert(std::numeric_limits<float>::is_iec559,
+              "TOSA Reference Model has not been built with standard IEE574 32-bit float support; ULP based "
+              "verifcation is invalid");
+static_assert(std::numeric_limits<double>::is_iec559,
+              "TOSA Reference Model has not been built with standard IEE574 64-bit float support; ULP based "
+              "verifcation is invalid");
+
+bool tosaCheckULP(float testValue, double referenceValue, int64_t ulpCount)
+{
+
+    // Start by sanitizing the input.
+
+    // The concept of ULP isn't defined for NaN's
+    if (std::isnan(referenceValue) || std::isnan(testValue))
+    {
+        return false;
+    }
+
+    // Make the sign of the reference value positive
+    // and adjust the test value appropriately.
+    if (referenceValue < 0)
+    {
+        referenceValue = -referenceValue;
+        testValue      = -testValue;
+    }
+
+    // At this point we are ready to calculate the ULP bounds for the reference value.
+    double referenceMin, referenceMax;
+
+    // If the reference is infinity e.g. the result of an overflow the test value must
+    // be infinity of an appropriate sign.
+    if (std::isinf(referenceValue))
+    {
+        // We already canonicalized the input such that the reference value is positive
+        // so no need to check again here.
+        referenceMin = std::numeric_limits<float>::infinity();
+        referenceMax = std::numeric_limits<float>::infinity();
+    }
+    else if (referenceValue == 0)
+    {
+        // For zero we require that the results match exactly with the correct sign.
+        referenceMin = 0;
+        referenceMax = 0;
+    }
+    else
+    {
+        // Find the exponent of the reference value.
+        int referenceExponent;
+        std::frexp(referenceValue, &referenceExponent);
+
+        // Work out the values magnitude - by raising 2 to the power of the
+        // exponent and taking the normalized minimum for denormal values
+        const double referencePower2 =
+            std::max(std::ldexp(1.0, referenceExponent), static_cast<double>(std::numeric_limits<float>::min()));
+        // Get the value of changing the last bit - by shifting the least significant bit to this magnitude
+        // i.e. the ULP.
+        double ulpValue = referencePower2 * std::ldexp(1.0, -23);
+
+        // It is possible that within one ULP we cross a boundary where we need to change the exponent,
+        // if this happens we will take the ULP for the larger exponent.
+        if (referenceValue + ulpValue > 2 * referencePower2)
+        {
+            ulpValue = 2 * ulpValue;
+        }
+
+        // Scale by the number of ULPs requested by the user.
+        referenceMax = referenceValue + ulpValue * ulpCount;
+        referenceMin = referenceValue - ulpValue * ulpCount;
+
+        // Handle the overflow cases.
+        if (referenceMax > std::numeric_limits<float>::max())
+        {
+            referenceMax = std::numeric_limits<float>::infinity();
+        }
+
+        if (referenceMin > std::numeric_limits<float>::max())
+        {
+            referenceMin = std::numeric_limits<float>::infinity();
+        }
+
+        // And the underflow cases.
+        if (referenceMax < std::numeric_limits<float>::min())
+        {
+            referenceMax = std::numeric_limits<float>::min();
+        }
+
+        if (referenceMin < std::numeric_limits<float>::min())
+        {
+            referenceMin = 0;
+        }
+    }
+
+    // And finally... Do the comparison.
+    return static_cast<double>(testValue) >= referenceMin && static_cast<double>(testValue) <= referenceMax;
+}
+}    // namespace
+
+bool verifyULP(const CTensor* referenceTensor, const CTensor* implementationTensor, uint64_t ulp)
+{
+    // Validate that tensors are provided
+    TOSA_REF_REQUIRE(referenceTensor != nullptr, "reference tensor is missing");
+    TOSA_REF_REQUIRE(implementationTensor != nullptr, "implementation tensor is missing");
+
+    // Get number of elements
+    const auto elementCount =
+        numElements(std::vector<int32_t>(referenceTensor->shape, referenceTensor->shape + referenceTensor->num_dims));
+    TOSA_REF_REQUIRE(elementCount > 0, "invalid shape for reference tensor");
+
+    switch (implementationTensor->data_type)
+    {
+        case tosa_datatype_fp32_t: {
+            const auto* refData = reinterpret_cast<const float*>(referenceTensor->data);
+            TOSA_REF_REQUIRE(refData != nullptr, "missing data for reference");
+            const auto* impData = reinterpret_cast<const float*>(implementationTensor->data);
+            TOSA_REF_REQUIRE(impData != nullptr, "missing data for implementation");
+            return std::equal(refData, std::next(refData, elementCount), impData, std::next(impData, elementCount),
+                              [ulp](const auto& referenceValue, const auto& implementationValue) {
+                                  return tosaCheckULP(referenceValue, implementationValue, ulp);
+                              });
+        }
+        default:
+            break;
+    }
+
+    return false;
+}
+}    // namespace TosaReference
diff --git a/reference_model/src/verify/verify_utils.h b/reference_model/src/verify/verify_utils.h
index 3a527da..510b9cb 100644
--- a/reference_model/src/verify/verify_utils.h
+++ b/reference_model/src/verify/verify_utils.h
@@ -51,7 +51,7 @@ struct UlpInfo
 {
     UlpInfo() = default;
 
-    float ulp;
+    uint64_t ulp;
 };
 
 /// \brief Dot-product verification meta-data