[libc] Fix a bug in MPFRUtils making ULP values off by 2^(-mantissaWidth).

[lntue]

Summary

Fix a bug in MPFRUtils making ULP values off by 2^(-mantissaWidth) and incorrect eps for denormal numbers.

[sivachandra]

Thanks for catching! My mistake.

libc/utils/MPFRWrapper/MPFRUtils.cpp

315

Add a comment saying this is for the denormal adjustment.

340

Same here.

[lntue]

[libc] Fix a bug in MPFRUtils making ULP values off by 2^(-mantissaWidth).

libc/utils/MPFRWrapper/MPFRUtils.cpp

public:
   template <typename T>
   cpp::EnableIfType<cpp::IsFloatingPointType<T>::Value, double> ulp(T input) {
     T thisAsT = as<T>();
     if (thisAsT == input)
       return T(0.0);
 
     int thisExponent = fputil::FPBits<T>(thisAsT).getExponent();
     int inputExponent = fputil::FPBits<T>(input).getExponent();
+    // Adjust the exponents for denormal numbers.
+    if (fputil::FPBits<T>(thisAsT).getUnbiasedExponent() == 0)
+      ++thisExponent;
+    if (fputil::FPBits<T>(input).getUnbiasedExponent() == 0)
+      ++inputExponent;
+
     if (thisAsT * input < 0 || thisExponent == inputExponent) {
       MPFRNumber inputMPFR(input);
       mpfr_sub(inputMPFR.value, value, inputMPFR.value, MPFR_RNDN);
       mpfr_abs(inputMPFR.value, inputMPFR.value, MPFR_RNDN);
-      mpfr_mul_2si(inputMPFR.value, inputMPFR.value, -thisExponent, MPFR_RNDN);
+      mpfr_mul_2si(inputMPFR.value, inputMPFR.value,
+                   -thisExponent + int(fputil::MantissaWidth<T>::value),
+                   MPFR_RNDN);
       return inputMPFR.as<double>();
     }
 
     // If the control reaches here, it means that this number and input are
     // of the same sign but different exponent. In such a case, ULP error is
     // calculated as sum of two parts.
     thisAsT = std::abs(thisAsT);
     input = std::abs(input);
     T min = thisAsT > input ? input : thisAsT;
     T max = thisAsT > input ? thisAsT : input;
     int minExponent = fputil::FPBits<T>(min).getExponent();
     int maxExponent = fputil::FPBits<T>(max).getExponent();
+    // Adjust the exponents for denormal numbers.
+    if (fputil::FPBits<T>(min).getUnbiasedExponent() == 0)
+      ++minExponent;
+    if (fputil::FPBits<T>(max).getUnbiasedExponent() == 0)
+      ++maxExponent;
 
     MPFRNumber minMPFR(min);
     MPFRNumber maxMPFR(max);
 
     MPFRNumber pivot(uint32_t(1));
     mpfr_mul_2si(pivot.value, pivot.value, maxExponent, MPFR_RNDN);
 
     mpfr_sub(minMPFR.value, pivot.value, minMPFR.value, MPFR_RNDN);
-    mpfr_mul_2si(minMPFR.value, minMPFR.value, -minExponent, MPFR_RNDN);
+    mpfr_mul_2si(minMPFR.value, minMPFR.value,
+                 -minExponent + int(fputil::MantissaWidth<T>::value),
+                 MPFR_RNDN);
 
     mpfr_sub(maxMPFR.value, maxMPFR.value, pivot.value, MPFR_RNDN);
-    mpfr_mul_2si(maxMPFR.value, maxMPFR.value, -maxExponent, MPFR_RNDN);
+    mpfr_mul_2si(maxMPFR.value, maxMPFR.value,
+                 -maxExponent + int(fputil::MantissaWidth<T>::value),
+                 MPFR_RNDN);
 
     mpfr_add(minMPFR.value, minMPFR.value, maxMPFR.value, MPFR_RNDN);
     return minMPFR.as<double>();
   }
 };
 
 namespace internal {