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> MPFRNumber(XType x, const Tolerance &t) {
69     mpfr_init2(value, mpfrPrecision);
70     mpfr_set_zero(value, 1); // Set to positive zero.
71     MPFRNumber xExponent(fputil::FPBits<XType>(x).getExponent());
72     // E = 2^E
73     mpfr_exp2(xExponent.value, xExponent.value, MPFR_RNDN);
74     uint32_t bitMask = 1 << (t.width - 1);
75     for (int n = -t.basePrecision; bitMask > 0; bitMask >>= 1) {
76       --n;
77       if (t.bits & bitMask) {
78         // delta = -n
79         MPFRNumber delta(n);
80 
81         // delta = 2^(-n)
82         mpfr_exp2(delta.value, delta.value, MPFR_RNDN);
83 
84         // delta = E * 2^(-n)
85         mpfr_mul(delta.value, delta.value, xExponent.value, MPFR_RNDN);
86 
87         // tolerance += delta
88         mpfr_add(value, value, delta.value, MPFR_RNDN);
89       }
90     }
91   }
92 
93   template <typename XType,
94             cpp::EnableIfType<cpp::IsFloatingPointType<XType>::Value, int> = 0>
95   MPFRNumber(Operation op, XType rawValue) {
96     mpfr_init2(value, mpfrPrecision);
97     MPFRNumber mpfrInput(rawValue);
98     switch (op) {
99     case Operation::Abs:
100       mpfr_abs(value, mpfrInput.value, MPFR_RNDN);
101       break;
102     case Operation::Ceil:
103       mpfr_ceil(value, mpfrInput.value);
104       break;
105     case Operation::Cos:
106       mpfr_cos(value, mpfrInput.value, MPFR_RNDN);
107       break;
108     case Operation::Exp:
109       mpfr_exp(value, mpfrInput.value, MPFR_RNDN);
110       break;
111     case Operation::Exp2:
112       mpfr_exp2(value, mpfrInput.value, MPFR_RNDN);
113       break;
114     case Operation::Floor:
115       mpfr_floor(value, mpfrInput.value);
116       break;
117     case Operation::Round:
118       mpfr_round(value, mpfrInput.value);
119       break;
120     case Operation::Sin:
121       mpfr_sin(value, mpfrInput.value, MPFR_RNDN);
122       break;
123     case Operation::Sqrt:
124       mpfr_sqrt(value, mpfrInput.value, MPFR_RNDN);
125       break;
126     case Operation::Trunc:
127       mpfr_trunc(value, mpfrInput.value);
128       break;
129     }
130   }
131 
132   MPFRNumber(const MPFRNumber &other) {
133     mpfr_set(value, other.value, MPFR_RNDN);
134   }
135 
136   ~MPFRNumber() { mpfr_clear(value); }
137 
138   // Returns true if |other| is within the |tolerance| value of this
139   // number.
140   bool isEqual(const MPFRNumber &other, const MPFRNumber &tolerance) const {
141     MPFRNumber difference;
142     if (mpfr_cmp(value, other.value) >= 0)
143       mpfr_sub(difference.value, value, other.value, MPFR_RNDN);
144     else
145       mpfr_sub(difference.value, other.value, value, MPFR_RNDN);
146 
147     return mpfr_lessequal_p(difference.value, tolerance.value);
148   }
149 
150   std::string str() const {
151     // 200 bytes should be more than sufficient to hold a 100-digit number
152     // plus additional bytes for the decimal point, '-' sign etc.
153     constexpr size_t printBufSize = 200;
154     char buffer[printBufSize];
155     mpfr_snprintf(buffer, printBufSize, "%100.50Rf", value);
156     llvm::StringRef ref(buffer);
157     ref = ref.trim();
158     return ref.str();
159   }
160 
161   // These functions are useful for debugging.
162   template <typename T> T as() const;
163 
164   template <> float as<float>() const { return mpfr_get_flt(value, MPFR_RNDN); }
165   template <> double as<double>() const { return mpfr_get_d(value, MPFR_RNDN); }
166   template <> long double as<long double>() const {
167     return mpfr_get_ld(value, MPFR_RNDN);
168   }
169 
170   void dump(const char *msg) const { mpfr_printf("%s%.128Rf\n", msg, value); }
171 
172   // Return the ULP (units-in-the-last-place) difference between the
173   // stored MPFR and a floating point number.
174   //
175   // We define:
176   //   ULP(mpfr_value, value) = abs(mpfr_value - value) / eps(value)
177   //
178   // Remarks:
179   // 1. ULP < 0.5 will imply that the value is correctly rounded.
180   // 2. We expect that this value and the value to be compared (the [input]
181   //    argument) are reasonable close, and we will provide an upper bound
182   //    of ULP value for testing.  Morever, most of the fractional parts of
183   //    ULP value do not matter much, so using double as the return type
184   //    should be good enough.
185   template <typename T>
186   cpp::EnableIfType<cpp::IsFloatingPointType<T>::Value, double> ulp(T input) {
187     fputil::FPBits<T> bits(input);
188     MPFRNumber mpfrInput(input);
189 
190     // abs(value - input)
191     mpfr_sub(mpfrInput.value, value, mpfrInput.value, MPFR_RNDN);
192     mpfr_abs(mpfrInput.value, mpfrInput.value, MPFR_RNDN);
193 
194     // get eps(input)
195     int epsExponent = bits.exponent - fputil::FPBits<T>::exponentBias -
196                       fputil::MantissaWidth<T>::value;
197     if (bits.exponent == 0) {
198       // correcting denormal exponent
199       ++epsExponent;
200     } else if ((bits.mantissa == 0) && (bits.exponent > 1) &&
201                mpfr_less_p(value, mpfrInput.value)) {
202       // when the input is exactly 2^n, distance (epsilon) between the input
203       // and the next floating point number is different from the distance to
204       // the previous floating point number.  So in that case, if the correct
205       // value from MPFR is smaller than the input, we use the smaller epsilon
206       --epsExponent;
207     }
208 
209     // Since eps(value) is of the form 2^e, instead of dividing such number,
210     // we multiply by its inverse 2^{-e}.
211     mpfr_mul_2si(mpfrInput.value, mpfrInput.value, -epsExponent, MPFR_RNDN);
212 
213     return mpfrInput.as<double>();
214   }
215 };
216 
217 namespace internal {
218 
219 template <typename T>
220 void MPFRMatcher<T>::explainError(testutils::StreamWrapper &OS) {
221   MPFRNumber mpfrResult(operation, input);
222   MPFRNumber mpfrInput(input);
223   MPFRNumber mpfrMatchValue(matchValue);
224   FPBits<T> inputBits(input);
225   FPBits<T> matchBits(matchValue);
226   FPBits<T> mpfrResultBits(mpfrResult.as<T>());
227   OS << "Match value not within tolerance value of MPFR result:\n"
228      << "  Input decimal: " << mpfrInput.str() << '\n';
229   __llvm_libc::fputil::testing::describeValue("     Input bits: ", input, OS);
230   OS << '\n' << "  Match decimal: " << mpfrMatchValue.str() << '\n';
231   __llvm_libc::fputil::testing::describeValue("     Match bits: ", matchValue,
232                                               OS);
233   OS << '\n' << "    MPFR result: " << mpfrResult.str() << '\n';
234   __llvm_libc::fputil::testing::describeValue(
235       "   MPFR rounded: ", mpfrResult.as<T>(), OS);
236   OS << '\n';
237   if (useULP) {
238     OS << "      ULP error: " << std::to_string(mpfrResult.ulp(matchValue))
239        << '\n';
240   } else {
241     MPFRNumber mpfrToleranceValue = MPFRNumber(matchValue, tolerance);
242     OS << "Tolerance value: " << mpfrToleranceValue.str() << '\n';
243   }
244 }
245 
246 template void MPFRMatcher<float>::explainError(testutils::StreamWrapper &);
247 template void MPFRMatcher<double>::explainError(testutils::StreamWrapper &);
248 template void
249 MPFRMatcher<long double>::explainError(testutils::StreamWrapper &);
250 
251 template <typename T>
252 bool compare(Operation op, T input, T libcResult, const Tolerance &t) {
253   MPFRNumber mpfrResult(op, input);
254   MPFRNumber mpfrLibcResult(libcResult);
255   MPFRNumber mpfrToleranceValue(libcResult, t);
256 
257   return mpfrResult.isEqual(mpfrLibcResult, mpfrToleranceValue);
258 };
259 
260 template bool compare<float>(Operation, float, float, const Tolerance &);
261 template bool compare<double>(Operation, double, double, const Tolerance &);
262 template bool compare<long double>(Operation, long double, long double,
263                                    const Tolerance &);
264 
265 template <typename T>
266 bool compare(Operation op, T input, T libcResult, double ulpError) {
267   // If the ulp error is exactly 0.5 (i.e a tie), we would check that the result
268   // is rounded to the nearest even.
269   MPFRNumber mpfrResult(op, input);
270   double ulp = mpfrResult.ulp(libcResult);
271   bool bitsAreEven = ((FPBits<T>(libcResult).bitsAsUInt() & 1) == 0);
272   return (ulp < ulpError) ||
273          ((ulp == ulpError) && ((ulp != 0.5) || bitsAreEven));
274 }
275 
276 template bool compare<float>(Operation, float, float, double);
277 template bool compare<double>(Operation, double, double, double);
278 template bool compare<long double>(Operation, long double, long double, double);
279 
280 } // namespace internal
281 
282 } // namespace mpfr
283 } // namespace testing
284 } // namespace __llvm_libc
285