[libc-commits] [libc] d5700bb - [libc] Calculate ulp error after rounding MPFR result to the result type.

[Siva Chandra Reddy via libc-commits]

Author: Siva Chandra Reddy
Date: 2021-06-23T20:29:46Z
New Revision: d5700bb694490b7e2a3c5dbd29d2a82b0843c8e9

URL: https://github.com/llvm/llvm-project/commit/d5700bb694490b7e2a3c5dbd29d2a82b0843c8e9
DIFF: https://github.com/llvm/llvm-project/commit/d5700bb694490b7e2a3c5dbd29d2a82b0843c8e9.diff

LOG: [libc] Calculate ulp error after rounding MPFR result to the result type.

Reviewed By: lntue

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

Added: 
    

Modified: 
    libc/utils/MPFRWrapper/MPFRUtils.cpp

Removed: 
    


################################################################################
diff  --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp
index e040393bc26c..959701cf9412 100644
--- a/libc/utils/MPFRWrapper/MPFRUtils.cpp
+++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp
@@ -270,9 +270,15 @@ class MPFRNumber {
   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);
+    // When calculating error, we first round this number to the floating
+    // point type T and calculate the ulp error against this rounded number.
+    // The input is always either exact or rounded. So, we if compare
+    // with this number directly, we can end up with a large ULP error.
+    MPFRNumber thisRoundedToFloat(as<T>());
+
+    // abs(thisRoundedToFloat - input)
+    mpfr_sub(mpfrInput.value, thisRoundedToFloat.value, mpfrInput.value,
+             MPFR_RNDN);
     mpfr_abs(mpfrInput.value, mpfrInput.value, MPFR_RNDN);
 
     // get eps(input)
@@ -282,7 +288,7 @@ class MPFRNumber {
       // correcting denormal exponent
       ++epsExponent;
     } else if ((bits.encoding.mantissa == 0) && (bits.encoding.exponent > 1) &&
-               mpfr_less_p(value, mpfrInput.value)) {
+               mpfr_less_p(thisRoundedToFloat.value, mpfrInput.value)) {
       // when the input is exactly 2^n, distance (epsilon) between the input
       // and the next floating point number is 
diff erent from the distance to
       // the previous floating point number.  So in that case, if the correct