1 //===-- Utils which wrap MPFR ---------------------------------------------===//
2 //
3 // Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
4 // See https://llvm.org/LICENSE.txt for license information.
5 // SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
6 //
7 //===----------------------------------------------------------------------===//
8 
9 #include "MPFRUtils.h"
10 
11 #include "utils/FPUtil/FPBits.h"
12 #include "utils/FPUtil/TestHelpers.h"
13 
14 #include "llvm/ADT/StringExtras.h"
15 #include "llvm/ADT/StringRef.h"
16 
17 #include <mpfr.h>
18 #include <stdint.h>
19 #include <string>
20 
21 template <typename T> using FPBits = __llvm_libc::fputil::FPBits<T>;
22 
23 namespace __llvm_libc {
24 namespace testing {
25 namespace mpfr {
26 
27 class MPFRNumber {
28   // A precision value which allows sufficiently large additional
29   // precision even compared to quad-precision floating point values.
30   static constexpr unsigned int mpfrPrecision = 128;
31 
32   mpfr_t value;
33 
34 public:
35   MPFRNumber() { mpfr_init2(value, mpfrPrecision); }
36 
37   // We use explicit EnableIf specializations to disallow implicit
38   // conversions. Implicit conversions can potentially lead to loss of
39   // precision.
40   template <typename XType,
41             cpp::EnableIfType<cpp::IsSame<float, XType>::Value, int> = 0>
42   explicit MPFRNumber(XType x) {
43     mpfr_init2(value, mpfrPrecision);
44     mpfr_set_flt(value, x, MPFR_RNDN);
45   }
46 
47   template <typename XType,
48             cpp::EnableIfType<cpp::IsSame<double, XType>::Value, int> = 0>
49   explicit MPFRNumber(XType x) {
50     mpfr_init2(value, mpfrPrecision);
51     mpfr_set_d(value, x, MPFR_RNDN);
52   }
53 
54   template <typename XType,
55             cpp::EnableIfType<cpp::IsSame<long double, XType>::Value, int> = 0>
56   explicit MPFRNumber(XType x) {
57     mpfr_init2(value, mpfrPrecision);
58     mpfr_set_ld(value, x, MPFR_RNDN);
59   }
60 
61   template <typename XType,
62             cpp::EnableIfType<cpp::IsIntegral<XType>::Value, int> = 0>
63   explicit MPFRNumber(XType x) {
64     mpfr_init2(value, mpfrPrecision);
65     mpfr_set_sj(value, x, MPFR_RNDN);
66   }
67 
68   template <typename XType,
69             cpp::EnableIfType<cpp::IsFloatingPointType<XType>::Value, int> = 0>
70   MPFRNumber(Operation op, XType rawValue) {
71     mpfr_init2(value, mpfrPrecision);
72     MPFRNumber mpfrInput(rawValue);
73     switch (op) {
74     case Operation::Abs:
75       mpfr_abs(value, mpfrInput.value, MPFR_RNDN);
76       break;
77     case Operation::Ceil:
78       mpfr_ceil(value, mpfrInput.value);
79       break;
80     case Operation::Cos:
81       mpfr_cos(value, mpfrInput.value, MPFR_RNDN);
82       break;
83     case Operation::Exp:
84       mpfr_exp(value, mpfrInput.value, MPFR_RNDN);
85       break;
86     case Operation::Exp2:
87       mpfr_exp2(value, mpfrInput.value, MPFR_RNDN);
88       break;
89     case Operation::Floor:
90       mpfr_floor(value, mpfrInput.value);
91       break;
92     case Operation::Round:
93       mpfr_round(value, mpfrInput.value);
94       break;
95     case Operation::Sin:
96       mpfr_sin(value, mpfrInput.value, MPFR_RNDN);
97       break;
98     case Operation::Sqrt:
99       mpfr_sqrt(value, mpfrInput.value, MPFR_RNDN);
100       break;
101     case Operation::Trunc:
102       mpfr_trunc(value, mpfrInput.value);
103       break;
104     }
105   }
106 
107   MPFRNumber(const MPFRNumber &other) {
108     mpfr_set(value, other.value, MPFR_RNDN);
109   }
110 
111   ~MPFRNumber() { mpfr_clear(value); }
112 
113   std::string str() const {
114     // 200 bytes should be more than sufficient to hold a 100-digit number
115     // plus additional bytes for the decimal point, '-' sign etc.
116     constexpr size_t printBufSize = 200;
117     char buffer[printBufSize];
118     mpfr_snprintf(buffer, printBufSize, "%100.50Rf", value);
119     llvm::StringRef ref(buffer);
120     ref = ref.trim();
121     return ref.str();
122   }
123 
124   // These functions are useful for debugging.
125   template <typename T> T as() const;
126 
127   template <> float as<float>() const { return mpfr_get_flt(value, MPFR_RNDN); }
128   template <> double as<double>() const { return mpfr_get_d(value, MPFR_RNDN); }
129   template <> long double as<long double>() const {
130     return mpfr_get_ld(value, MPFR_RNDN);
131   }
132 
133   void dump(const char *msg) const { mpfr_printf("%s%.128Rf\n", msg, value); }
134 
135   // Return the ULP (units-in-the-last-place) difference between the
136   // stored MPFR and a floating point number.
137   //
138   // We define:
139   //   ULP(mpfr_value, value) = abs(mpfr_value - value) / eps(value)
140   //
141   // Remarks:
142   // 1. ULP < 0.5 will imply that the value is correctly rounded.
143   // 2. We expect that this value and the value to be compared (the [input]
144   //    argument) are reasonable close, and we will provide an upper bound
145   //    of ULP value for testing.  Morever, most of the fractional parts of
146   //    ULP value do not matter much, so using double as the return type
147   //    should be good enough.
148   template <typename T>
149   cpp::EnableIfType<cpp::IsFloatingPointType<T>::Value, double> ulp(T input) {
150     fputil::FPBits<T> bits(input);
151     MPFRNumber mpfrInput(input);
152 
153     // abs(value - input)
154     mpfr_sub(mpfrInput.value, value, mpfrInput.value, MPFR_RNDN);
155     mpfr_abs(mpfrInput.value, mpfrInput.value, MPFR_RNDN);
156 
157     // get eps(input)
158     int epsExponent = bits.exponent - fputil::FPBits<T>::exponentBias -
159                       fputil::MantissaWidth<T>::value;
160     if (bits.exponent == 0) {
161       // correcting denormal exponent
162       ++epsExponent;
163     } else if ((bits.mantissa == 0) && (bits.exponent > 1) &&
164                mpfr_less_p(value, mpfrInput.value)) {
165       // when the input is exactly 2^n, distance (epsilon) between the input
166       // and the next floating point number is different from the distance to
167       // the previous floating point number.  So in that case, if the correct
168       // value from MPFR is smaller than the input, we use the smaller epsilon
169       --epsExponent;
170     }
171 
172     // Since eps(value) is of the form 2^e, instead of dividing such number,
173     // we multiply by its inverse 2^{-e}.
174     mpfr_mul_2si(mpfrInput.value, mpfrInput.value, -epsExponent, MPFR_RNDN);
175 
176     return mpfrInput.as<double>();
177   }
178 };
179 
180 namespace internal {
181 
182 template <typename T>
183 void MPFRMatcher<T>::explainError(testutils::StreamWrapper &OS) {
184   MPFRNumber mpfrResult(operation, input);
185   MPFRNumber mpfrInput(input);
186   MPFRNumber mpfrMatchValue(matchValue);
187   FPBits<T> inputBits(input);
188   FPBits<T> matchBits(matchValue);
189   FPBits<T> mpfrResultBits(mpfrResult.as<T>());
190   OS << "Match value not within tolerance value of MPFR result:\n"
191      << "  Input decimal: " << mpfrInput.str() << '\n';
192   __llvm_libc::fputil::testing::describeValue("     Input bits: ", input, OS);
193   OS << '\n' << "  Match decimal: " << mpfrMatchValue.str() << '\n';
194   __llvm_libc::fputil::testing::describeValue("     Match bits: ", matchValue,
195                                               OS);
196   OS << '\n' << "    MPFR result: " << mpfrResult.str() << '\n';
197   __llvm_libc::fputil::testing::describeValue(
198       "   MPFR rounded: ", mpfrResult.as<T>(), OS);
199   OS << '\n';
200   OS << "      ULP error: " << std::to_string(mpfrResult.ulp(matchValue))
201      << '\n';
202 }
203 
204 template void MPFRMatcher<float>::explainError(testutils::StreamWrapper &);
205 template void MPFRMatcher<double>::explainError(testutils::StreamWrapper &);
206 template void
207 MPFRMatcher<long double>::explainError(testutils::StreamWrapper &);
208 
209 template <typename T>
210 bool compare(Operation op, T input, T libcResult, double ulpError) {
211   // If the ulp error is exactly 0.5 (i.e a tie), we would check that the result
212   // is rounded to the nearest even.
213   MPFRNumber mpfrResult(op, input);
214   double ulp = mpfrResult.ulp(libcResult);
215   bool bitsAreEven = ((FPBits<T>(libcResult).bitsAsUInt() & 1) == 0);
216   return (ulp < ulpError) ||
217          ((ulp == ulpError) && ((ulp != 0.5) || bitsAreEven));
218 }
219 
220 template bool compare<float>(Operation, float, float, double);
221 template bool compare<double>(Operation, double, double, double);
222 template bool compare<long double>(Operation, long double, long double, double);
223 
224 } // namespace internal
225 
226 } // namespace mpfr
227 } // namespace testing
228 } // namespace __llvm_libc
229