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

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.

Reviewers: abrachet

Differential Revision: https://reviews.llvm.org/D79278

diff --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp
index 75ee2ad..980557a 100644 (file)
--- a/libc/utils/MPFRWrapper/MPFRUtils.cpp
+++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp
@@ -8,15 +8,55 @@
 
 #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.
@@ -45,6 +85,38 @@ public:
   }
 
   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);
@@ -65,20 +137,9 @@ public:
 
   ~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);
@@ -112,10 +173,14 @@ 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 &);
@@ -124,9 +189,10 @@ 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 &);