[libc] Add ULP function to MPFRNumber class to test correctly rounded functions such as SQRT, FMA.

[lntue]

Summary

Add ULP function to MPFRNumber class to test correctly rounded functions.

[lntue]

[libc] Add ULP function to MPFRNumber class to test correctly rounded functions such as SQRT, FMA.

[lntue]

Fix explicit specialization in class scope.

[sivachandra]

The functionality LGTM. May be a few more comments would have helped me read this faster, but never mind for now. I am want to spend a little more time the layering of FPUtils, MPFRWrapper and UnitTest. Will come back soon with more comments.

[sivachandra]

We currently have UnitTest depend on FPUtil. I think that an creates unnecessary dependency for tests which do not need FPUtil for testing. So, a plan could be:

Instead of EXPECT_EQ and friends for floating point numbers, we will decouple UnitTest from FPUtil using matchers implemented in FPUtil. Say, EXPECT_FP_MATCH etc.

WDYT? If you concur, then you need not do the work. I will do it, and you can absorb those changes into this patch.

[lntue]

In D84725#2183801, @sivachandra wrote:

We currently have UnitTest depend on FPUtil. I think that an creates unnecessary dependency for tests which do not need FPUtil for testing. So, a plan could be:

Instead of EXPECT_EQ and friends for floating point numbers, we will decouple UnitTest from FPUtil using matchers implemented in FPUtil. Say, EXPECT_FP_MATCH etc.

WDYT? If you concur, then you need not do the work. I will do it, and you can absorb those changes into this patch.

That sounds good to me. Let me know when it's ready. Thanks!

[lntue]

Merge with reviews.llvm.org/D85486

[sivachandra]

libc/utils/MPFRWrapper/MPFRUtils.cpp

176

Few comments about documentation.

1. I think this function assumes that the exponent of input matches that of the MPFR number. If yes, document it here?
2. Explain why a double return value is good enough?
3. Lastly, is the name ulp of this method correct? I guess, because you defined ULP above as a ratio, the name ulp is correct and is OK as is.

219

Instead of non-member functions, can we replace the existing methods with a template method:

template <typename T>
T as() const;

We will still need specializations like below though.

235

In what circumstances would a rounded MPFR result help understand a test failure?

238

Is plural errors correct here?

[lntue]

libc/utils/MPFRWrapper/MPFRUtils.cpp

235

When the ULP error is exactly 0.5, the ULP value itself is not enough to check if the rounding is correct. This should never happen with SQRT function, but other one like HYPOT or FMA might be able to hit exactly 0.5 ULP.

[lntue]

Correct the eps when the input is a power of 2, and add a check for correct rounding to even when it is a tie.

[lntue]

Update comments

[lntue]

libc/utils/MPFRWrapper/MPFRUtils.cpp

176

I updated the comments to reflect the expectations.

235

I updated the test to reflect the case when there is a tie in rounded, i.e. ULP error == 0.5, we expect it to be rounded to even

[sivachandra]

libc/utils/MPFRWrapper/MPFRUtils.cpp

180

May be reword it as:

We expect that this value and the value to be compared (the |input| argument) are reasonably close ...

[lntue]

Update comments.

libc/utils/MPFRWrapper/MPFRUtils.h

enum class Operation : int {
   Abs,
   Ceil,
   Cos,
   Exp,
   Exp2,
   Floor,
   Round,
   Sin,
+  Sqrt,
   Trunc
 };
 
 namespace internal {
 
 template <typename T>
 bool compare(Operation op, T input, T libcOutput, const Tolerance &t);
 
+template <typename T>
+bool compare(Operation op, T input, T libcOutput, double t);
+
 template <typename T> class MPFRMatcher : public testing::Matcher<T> {
   static_assert(__llvm_libc::cpp::IsFloatingPointType<T>::Value,
                 "MPFRMatcher can only be used with floating point values.");
 
   Operation operation;
   T input;
   Tolerance tolerance;
   T matchValue;
+  double ulpTolerance;
+  bool useULP;
 
 public:
   MPFRMatcher(Operation op, T testInput, Tolerance &t)
-      : operation(op), input(testInput), tolerance(t) {}
+      : operation(op), input(testInput), tolerance(t), useULP(false) {}
+  MPFRMatcher(Operation op, T testInput, double ulpTolerance)
+      : operation(op), input(testInput), ulpTolerance(ulpTolerance),
+        useULP(true) {}
 
   bool match(T libcResult) {
     matchValue = libcResult;
-    return internal::compare(operation, input, libcResult, tolerance);
+    return (useULP
+                ? internal::compare(operation, input, libcResult, ulpTolerance)
+                : internal::compare(operation, input, libcResult, tolerance));
   }
 
   void explainError(testutils::StreamWrapper &OS) override;
 };
 
 } // namespace internal
 
-template <typename T>
+template <typename T, typename U>
 __attribute__((no_sanitize("address")))
-internal::MPFRMatcher<T> getMPFRMatcher(Operation op, T input, Tolerance t) {
+typename cpp::EnableIfType<cpp::IsSameV<U, Tolerance> ||
+                               cpp::IsSameV<U, double>,
+                           internal::MPFRMatcher<T>>
+getMPFRMatcher(Operation op, T input, U t) {
   static_assert(
       __llvm_libc::cpp::IsFloatingPointType<T>::Value,
       "getMPFRMatcher can only be used to match floating point results.");
   return internal::MPFRMatcher<T>(op, input, t);
 }
 
 } // namespace mpfr
 } // namespace testing

libc/utils/MPFRWrapper/MPFRUtils.cpp

 //===-- Utils which wrap MPFR ---------------------------------------------===//
 //
 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
 // See https://llvm.org/LICENSE.txt for license information.
 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
 //
 //===----------------------------------------------------------------------===//
 
 #include "MPFRUtils.h"
 
 #include "utils/FPUtil/FPBits.h"
+#include "utils/FPUtil/TestHelpers.h"
 
 #include "llvm/ADT/StringExtras.h"
 #include "llvm/ADT/StringRef.h"
 
 #include <mpfr.h>
 #include <stdint.h>
 #include <string>
 
case Operation::Floor:
       mpfr_floor(value, mpfrInput.value);
       break;
     case Operation::Round:
       mpfr_round(value, mpfrInput.value);
       break;
     case Operation::Sin:
       mpfr_sin(value, mpfrInput.value, MPFR_RNDN);
       break;
+    case Operation::Sqrt:
+      mpfr_sqrt(value, mpfrInput.value, MPFR_RNDN);
+      break;
     case Operation::Trunc:
       mpfr_trunc(value, mpfrInput.value);
       break;
     }
   }
 
   MPFRNumber(const MPFRNumber &other) {
     mpfr_set(value, other.value, MPFR_RNDN);
std::string str() const {
     char buffer[printBufSize];
     mpfr_snprintf(buffer, printBufSize, "%100.50Rf", value);
     llvm::StringRef ref(buffer);
     ref = ref.trim();
     return ref.str();
   }
 
   // These functions are useful for debugging.
-  float asFloat() const { return mpfr_get_flt(value, MPFR_RNDN); }
-  double asDouble() const { return mpfr_get_d(value, MPFR_RNDN); }
+  template <typename T> T as() const;
+
+  template <> float as<float>() const { return mpfr_get_flt(value, MPFR_RNDN); }
+  template <> double as<double>() const { return mpfr_get_d(value, MPFR_RNDN); }
+  template <> long double as<long double>() const {
+    return mpfr_get_ld(value, MPFR_RNDN);
+  }
+
   void dump(const char *msg) const { mpfr_printf("%s%.128Rf\n", msg, value); }
+
+  // Return the ULP (units-in-the-last-place) difference between the
+  // stored MPFR and a floating point number.
+  //
+  // We define:
+  //   ULP(mpfr_value, value) = abs(mpfr_value - value) / eps(value)
+  //
+  // Remarks:
+  // 1. ULP < 0.5 will imply that the value is correctly rounded.
+  // 2. We expect that this value and the value to be compared (the [input]
+  //    argument) are reasonable close, and we will provide an upper bound
+  //    of ULP value for testing.  Morever, most of the fractional parts of
+  //    ULP value do not matter much, so using double as the return type
+  //    should be good enough.
+  template <typename T>
+  cpp::EnableIfType<cpp::IsFloatingPointType<T>::Value, double> ulp(T input) {
+    fputil::FPBits<T> bits(input);
+    MPFRNumber mpfrInput(input);
+
+    // abs(value - input)
+    mpfr_sub(mpfrInput.value, value, mpfrInput.value, MPFR_RNDN);
+    mpfr_abs(mpfrInput.value, mpfrInput.value, MPFR_RNDN);
+
+    // get eps(input)
+    int epsExponent = bits.exponent - fputil::FPBits<T>::exponentBias -
+                      fputil::MantissaWidth<T>::value;
+    if (bits.exponent == 0) {
+      // correcting denormal exponent
+      ++epsExponent;
+    } else if ((bits.mantissa == 0) && (bits.exponent > 1) &&
+               mpfr_less_p(value, mpfrInput.value)) {
+      // when the input is exactly 2^n, distance (epsilon) between the input
+      // and the next floating point number is different from the distance to
+      // the previous floating point number.  So in that case, if the correct
+      // value from MPFR is smaller than the input, we use the smaller epsilon
+      --epsExponent;
+    }
+
+    // Since eps(value) is of the form 2^e, instead of dividing such number,
+    // we multiply by its inverse 2^{-e}.
+    mpfr_mul_2si(mpfrInput.value, mpfrInput.value, -epsExponent, MPFR_RNDN);
+
+    return mpfrInput.as<double>();
+  }
 };
 
 namespace internal {
 
 template <typename T>
 void MPFRMatcher<T>::explainError(testutils::StreamWrapper &OS) {
   MPFRNumber mpfrResult(operation, input);
   MPFRNumber mpfrInput(input);
   MPFRNumber mpfrMatchValue(matchValue);
-  MPFRNumber mpfrToleranceValue(matchValue, tolerance);
   FPBits<T> inputBits(input);
   FPBits<T> matchBits(matchValue);
-  // TODO: Call to llvm::utohexstr implicitly converts __uint128_t values to
-  // uint64_t values. This can be fixed using a custom wrapper for
-  // llvm::utohexstr to handle __uint128_t values correctly.
+  FPBits<T> mpfrResultBits(mpfrResult.as<T>());
   OS << "Match value not within tolerance value of MPFR result:\n"
-     << "  Input decimal: " << mpfrInput.str() << '\n'
-     << "     Input bits: 0x" << llvm::utohexstr(inputBits.bitsAsUInt()) << '\n'
-     << "  Match decimal: " << mpfrMatchValue.str() << '\n'
-     << "     Match bits: 0x" << llvm::utohexstr(matchBits.bitsAsUInt()) << '\n'
-     << "    MPFR result: " << mpfrResult.str() << '\n'
-     << "Tolerance value: " << mpfrToleranceValue.str() << '\n';
+     << "  Input decimal: " << mpfrInput.str() << '\n';
+  __llvm_libc::fputil::testing::describeValue("     Input bits: ", input, OS);
+  OS << '\n' << "  Match decimal: " << mpfrMatchValue.str() << '\n';
+  __llvm_libc::fputil::testing::describeValue("     Match bits: ", matchValue,
+                                              OS);
+  OS << '\n' << "    MPFR result: " << mpfrResult.str() << '\n';
+  __llvm_libc::fputil::testing::describeValue(
+      "   MPFR rounded: ", mpfrResult.as<T>(), OS);
+  OS << '\n';
+  if (useULP) {
+    OS << "      ULP error: " << std::to_string(mpfrResult.ulp(matchValue))
+       << '\n';
+  } else {
+    MPFRNumber mpfrToleranceValue = MPFRNumber(matchValue, tolerance);
+    OS << "Tolerance value: " << mpfrToleranceValue.str() << '\n';
+  }
 }
 
 template void MPFRMatcher<float>::explainError(testutils::StreamWrapper &);
 template void MPFRMatcher<double>::explainError(testutils::StreamWrapper &);
 template void
 MPFRMatcher<long double>::explainError(testutils::StreamWrapper &);
 
 template <typename T>
 bool compare(Operation op, T input, T libcResult, const Tolerance &t) {
   MPFRNumber mpfrResult(op, input);
   MPFRNumber mpfrLibcResult(libcResult);
   MPFRNumber mpfrToleranceValue(libcResult, t);
 
   return mpfrResult.isEqual(mpfrLibcResult, mpfrToleranceValue);
 };
 
 template bool compare<float>(Operation, float, float, const Tolerance &);
 template bool compare<double>(Operation, double, double, const Tolerance &);
 template bool compare<long double>(Operation, long double, long double,
                                    const Tolerance &);
 
+template <typename T>
+bool compare(Operation op, T input, T libcResult, double ulpError) {
+  // If the ulp error is exactly 0.5 (i.e a tie), we would check that the result
+  // is rounded to the nearest even.
+  MPFRNumber mpfrResult(op, input);
+  double ulp = mpfrResult.ulp(libcResult);
+  bool bitsAreEven = ((FPBits<T>(libcResult).bitsAsUInt() & 1) == 0);
+  return (ulp < ulpError) ||
+         ((ulp == ulpError) && ((ulp != 0.5) || bitsAreEven));
+}
+
+template bool compare<float>(Operation, float, float, double);
+template bool compare<double>(Operation, double, double, double);
+template bool compare<long double>(Operation, long double, long double, double);
+
 } // namespace internal
 
 } // namespace mpfr
 } // namespace testing
 } // namespace __llvm_libc