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

6 files changed, 223 insertions, 1 deletions

diff --git a/reference_model/CMakeLists.txt b/reference_model/CMakeLists.txt
index 5facfe1..9392221 100644
--- a/reference_model/CMakeLists.txt
+++ b/reference_model/CMakeLists.txt
@@ -74,6 +74,7 @@ set(CXX_SOURCE
     src/verify/verify_dot_product.cc
     src/verify/verify_entry.cc
     src/verify/verify_exact.cc
+    src/verify/verify_ulp.cc
     src/verify/verify_utils.cc
     src/ops/op_factory.cc
     src/ops/tensor_ops.cc
@@ -143,6 +144,7 @@ add_library(tosa_reference_verify_lib SHARED
   src/verify/verify_dot_product.cc
   src/verify/verify_entry.cc
   src/verify/verify_exact.cc
+  src/verify/verify_ulp.cc
   src/verify/verify_utils.cc
   src/verify/verify_config.cc
   src/func_debug.cc
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
diff --git a/reference_model/test/verify_tests.cpp b/reference_model/test/verify_tests.cpp
index 731a808..7482847 100644
--- a/reference_model/test/verify_tests.cpp
+++ b/reference_model/test/verify_tests.cpp
@@ -14,6 +14,8 @@
 #include "verify.h"
 
 #include <algorithm>
+#include <cmath>
+#include <cstdint>
 #include <doctest.h>
 
 #include <array>
@@ -55,6 +57,14 @@ private:
     tosa_tensor_t _tensor;
 };
 
+template <typename FP>
+std::enable_if_t<std::is_floating_point_v<FP>, FP> increment(FP input, uint64_t steps)
+{
+    for (uint64_t step = 0; step < steps; ++step)
+        input = std::nextafter(input, std::numeric_limits<FP>::infinity());
+    return input;
+}
+
 auto& getRandomGenerator()
 {
     static std::mt19937 gen(0);
@@ -227,4 +237,49 @@ TEST_CASE("positive - exact")
     }
 }
 
+TEST_CASE("positive - ulp")
+{
+    std::string json_cfg = R"({
+        "tensors" : {
+            "out1" : {
+                "mode": "ULP",
+                "ulp_info": {
+                "ulp": 5
+                }
+            }
+        }
+    })";
+
+    const auto shape        = std::vector<int32_t>{ 8, 8, 8 };
+    const auto elementCount = std::accumulate(std::begin(shape), std::end(shape), 1, std::multiplies<>());
+
+    // Generate some random floats using the full range of fp32.
+    auto data = generateRandomTensorData<float>(elementCount, false);
+    SUBCASE("same")
+    {
+        // Generate some data that meets the ULP requirements of the result.
+        auto otherData = data;
+        std::for_each(std::begin(otherData), std::end(otherData), [](auto& value) { value = increment(value, 5); });
+        const auto referenceTensor =
+            TosaTensor("out1", tosa_datatype_fp64_t, shape, reinterpret_cast<uint8_t*>(data.data()));
+        const auto implementationTensor =
+            TosaTensor("out1", tosa_datatype_fp32_t, shape, reinterpret_cast<uint8_t*>(otherData.data()));
+        REQUIRE(tvf_verify_data(referenceTensor.cTensor(), nullptr, implementationTensor.cTensor(), json_cfg.c_str()));
+    }
+
+    SUBCASE("different")
+    {
+        // Generate some data that exceeds a specified number of ULP for each value in the tensor.
+        auto otherData = std::vector<float>(elementCount);
+        std::for_each(std::begin(otherData), std::end(otherData), [](auto& value) { value = increment(value, 6); });
+
+        const auto referenceTensor =
+            TosaTensor("out1", tosa_datatype_fp64_t, shape, reinterpret_cast<uint8_t*>(data.data()));
+        const auto implementationTensor =
+            TosaTensor("out1", tosa_datatype_fp32_t, shape, reinterpret_cast<uint8_t*>(otherData.data()));
+        REQUIRE_FALSE(
+            tvf_verify_data(referenceTensor.cTensor(), nullptr, implementationTensor.cTensor(), json_cfg.c_str()));
+    }
+}
+
 TEST_SUITE_END();    // verify