[libc] Fix how math results are compared with MPFR results.

[sivachandra]

Summary

Math results are compared with MPFR results by checking if they are
within a tolerance level of the MPFR result. The tolerance level is set
using additional bits of precision of the fractional part of a floating
point value. Hence, the actual value of the tolerance depends on not
only the additional bits, but also on the exponent part of the floating
point number.

Previously, the exponent part was not considered in evaluating the
tolerance value. While it was OK for small values less than 1 (hence
sinf, cosf, sincosf tests were OK), it breaks for large values which
functions like exp and friends produce. This change uses the exponent
value also to evaluate the tolerance value. LLVM libc produced results
can now be compared with MPFR produced results for large values also.

Depends on D79256.

[sivachandra]

Rebase.

[abrachet]

libc/utils/MPFRWrapper/MPFRUtils.cpp

29

Is there any way to identify the floating point standard being used? Is IEEE 754 just used everywhere?

30

Does it make semantic sense for signMask's type to be BitsType?

56

Period after comment.

191–192

I missed this before, but is this being used?

[sivachandra]

Address comments.

libc/utils/MPFRWrapper/MPFRUtils.cpp

29

Its just IEEE 754 used everywhere. If we add support for other formats, then we can specialize then may be.

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 "llvm/ADT/StringExtras.h"
 #include "llvm/ADT/StringRef.h"
 
 #include <mpfr.h>
+#include <stdint.h>
 #include <string>
 
 namespace __llvm_libc {
 namespace testing {
 namespace mpfr {
 
+template <typename T> struct FloatProperties {};
+
+template <> struct FloatProperties<float> {
+  typedef uint32_t BitsType;
+  static_assert(sizeof(BitsType) == sizeof(float),
+                "Unexpected size of 'float' type.");
+
+  static constexpr uint32_t mantissaWidth = 23;
+  static constexpr BitsType signMask = 0x7FFFFFFFU;
+  static constexpr uint32_t exponentOffset = 127;
+};
+
+template <> struct FloatProperties<double> {
+  typedef uint64_t BitsType;
+  static_assert(sizeof(BitsType) == sizeof(double),
+                "Unexpected size of 'double' type.");
+
+  static constexpr uint32_t mantissaWidth = 52;
+  static constexpr BitsType signMask = 0x7FFFFFFFFFFFFFFFULL;
+  static constexpr uint32_t exponentOffset = 1023;
+};
+
+template <typename T> typename FloatProperties<T>::BitsType getBits(T x) {
+  using BitsType = typename FloatProperties<T>::BitsType;
+  return *reinterpret_cast<BitsType *>(&x);
+}
+
+// Returns the zero adjusted exponent value of abs(x).
+template <typename T> int getExponent(T x) {
+  using Properties = FloatProperties<T>;
+  using BitsType = typename Properties::BitsType;
+  BitsType bits = *reinterpret_cast<BitsType *>(&x);
+  bits &= Properties::signMask;                // Zero the sign bit.
+  int e = (bits >> Properties::mantissaWidth); // Shift out the mantissa.
+  e -= Properties::exponentOffset;             // Zero adjust.
+  return e;
+}
+
 class MPFRNumber {
   // A precision value which allows sufficiently large additional
   // precision even compared to double precision floating point values.
   static constexpr unsigned int mpfrPrecision = 96;
 
   mpfr_t value;
 
 public:
public:
   template <typename XType,
             cpp::EnableIfType<cpp::IsSame<double, XType>::Value, int> = 0>
   explicit MPFRNumber(XType x) {
     mpfr_init2(value, mpfrPrecision);
     mpfr_set_d(value, x, MPFR_RNDN);
   }
 
   template <typename XType,
+            cpp::EnableIfType<cpp::IsIntegral<XType>::Value, int> = 0>
+  explicit MPFRNumber(XType x) {
+    mpfr_init2(value, mpfrPrecision);
+    mpfr_set_sj(value, x, MPFR_RNDN);
+  }
+
+  template <typename XType> MPFRNumber(XType x, const Tolerance &t) {
+    mpfr_init2(value, mpfrPrecision);
+    mpfr_set_zero(value, 1); // Set to positive zero.
+    MPFRNumber xExponent(getExponent(x));
+    // E = 2^E
+    mpfr_exp2(xExponent.value, xExponent.value, MPFR_RNDN);
+    uint32_t bitMask = 1 << (t.width - 1);
+    for (int n = -t.basePrecision; bitMask > 0; bitMask >>= 1) {
+      --n;
+      if (t.bits & bitMask) {
+        // delta = -n
+        MPFRNumber delta(n);
+
+        // delta = 2^(-n)
+        mpfr_exp2(delta.value, delta.value, MPFR_RNDN);
+
+        // delta = E * 2^(-n)
+        mpfr_mul(delta.value, delta.value, xExponent.value, MPFR_RNDN);
+
+        // tolerance += delta
+        mpfr_add(value, value, delta.value, MPFR_RNDN);
+      }
+    }
+  }
+
+  template <typename XType,
             cpp::EnableIfType<cpp::IsFloatingPointType<XType>::Value, int> = 0>
   MPFRNumber(Operation op, XType rawValue) {
     mpfr_init2(value, mpfrPrecision);
     MPFRNumber mpfrInput(rawValue);
     switch (op) {
     case OP_Cos:
       mpfr_cos(value, mpfrInput.value, MPFR_RNDN);
       break;
     case OP_Sin:
       mpfr_sin(value, mpfrInput.value, MPFR_RNDN);
       break;
     }
   }
 
   MPFRNumber(const MPFRNumber &other) {
     mpfr_set(value, other.value, MPFR_RNDN);
   }
 
   ~MPFRNumber() { mpfr_clear(value); }
 
-  // Returns true if |other| is within the tolerance value |t| of this
+  // Returns true if |other| is within the |tolerance| value of this
   // number.
-  bool isEqual(const MPFRNumber &other, const Tolerance &t) {
-    MPFRNumber tolerance(0.0);
-    uint32_t bitMask = 1 << (t.width - 1);
-    for (int exponent = -t.basePrecision; bitMask > 0; bitMask >>= 1) {
-      --exponent;
-      if (t.bits & bitMask) {
-        MPFRNumber delta;
-        mpfr_set_ui_2exp(delta.value, 1, exponent, MPFR_RNDN);
-        mpfr_add(tolerance.value, tolerance.value, delta.value, MPFR_RNDN);
-      }
-    }
-
+  bool isEqual(const MPFRNumber &other, const MPFRNumber &tolerance) const {
     MPFRNumber difference;
     if (mpfr_cmp(value, other.value) >= 0)
       mpfr_sub(difference.value, value, other.value, MPFR_RNDN);
     else
       mpfr_sub(difference.value, other.value, value, MPFR_RNDN);
 
     return mpfr_lessequal_p(difference.value, tolerance.value);
   }

 
 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);
   OS << "Match value not within tolerance value of MPFR result:\n"
-     << "Operation input: " << mpfrInput.str() << '\n'
-     << "    Match value: " << mpfrMatchValue.str() << '\n'
-     << "    MPFR result: " << mpfrResult.str() << '\n';
+     << "  Input decimal: " << mpfrInput.str() << '\n'
+     << "     Input bits: 0x" << llvm::utohexstr(getBits(input)) << '\n'
+     << "  Match decimal: " << mpfrMatchValue.str() << '\n'
+     << "     Match bits: 0x" << llvm::utohexstr(getBits(matchValue)) << '\n'
+     << "    MPFR result: " << mpfrResult.str() << '\n'
+     << "Tolerance value: " << mpfrToleranceValue.str() << '\n';
 }
 
 template void MPFRMatcher<float>::explainError(testutils::StreamWrapper &);
 template void MPFRMatcher<double>::explainError(testutils::StreamWrapper &);
 
 template <typename T>
 bool compare(Operation op, T input, T libcResult, const Tolerance &t) {
   MPFRNumber mpfrResult(op, input);
-  MPFRNumber mpfrInput(input);
   MPFRNumber mpfrLibcResult(libcResult);
-  return mpfrResult.isEqual(mpfrLibcResult, t);
+  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 &);
 
 } // namespace internal
 
 } // namespace mpfr
 } // namespace testing
 } // namespace __llvm_libc