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 <memory>
18 #include <stdint.h>
19 #include <string>
20 
21 #ifdef CUSTOM_MPFR_INCLUDER
22 // Some downstream repos are monoliths carrying MPFR sources in their third
23 // party directory. In such repos, including the MPFR header as
24 // `#include <mpfr.h>` is either disallowed or not possible. If that is the
25 // case, a file named `CustomMPFRIncluder.h` should be added through which the
26 // MPFR header can be included in manner allowed in that repo.
27 #include "CustomMPFRIncluder.h"
28 #else
29 #include <mpfr.h>
30 #endif
31 
32 template <typename T> using FPBits = __llvm_libc::fputil::FPBits<T>;
33 
34 namespace __llvm_libc {
35 namespace testing {
36 namespace mpfr {
37 
38 class MPFRNumber {
39   // A precision value which allows sufficiently large additional
40   // precision even compared to quad-precision floating point values.
41   static constexpr unsigned int mpfrPrecision = 128;
42 
43   mpfr_t value;
44 
45 public:
46   MPFRNumber() { mpfr_init2(value, mpfrPrecision); }
47 
48   // We use explicit EnableIf specializations to disallow implicit
49   // conversions. Implicit conversions can potentially lead to loss of
50   // precision.
51   template <typename XType,
52             cpp::EnableIfType<cpp::IsSame<float, XType>::Value, int> = 0>
53   explicit MPFRNumber(XType x) {
54     mpfr_init2(value, mpfrPrecision);
55     mpfr_set_flt(value, x, MPFR_RNDN);
56   }
57 
58   template <typename XType,
59             cpp::EnableIfType<cpp::IsSame<double, XType>::Value, int> = 0>
60   explicit MPFRNumber(XType x) {
61     mpfr_init2(value, mpfrPrecision);
62     mpfr_set_d(value, x, MPFR_RNDN);
63   }
64 
65   template <typename XType,
66             cpp::EnableIfType<cpp::IsSame<long double, XType>::Value, int> = 0>
67   explicit MPFRNumber(XType x) {
68     mpfr_init2(value, mpfrPrecision);
69     mpfr_set_ld(value, x, MPFR_RNDN);
70   }
71 
72   template <typename XType,
73             cpp::EnableIfType<cpp::IsIntegral<XType>::Value, int> = 0>
74   explicit MPFRNumber(XType x) {
75     mpfr_init2(value, mpfrPrecision);
76     mpfr_set_sj(value, x, MPFR_RNDN);
77   }
78 
79   MPFRNumber(const MPFRNumber &other) {
80     mpfr_set(value, other.value, MPFR_RNDN);
81   }
82 
83   ~MPFRNumber() {
84     mpfr_clear(value);
85   }
86 
87   MPFRNumber &operator=(const MPFRNumber &rhs) {
88     mpfr_set(value, rhs.value, MPFR_RNDN);
89     return *this;
90   }
91 
92   MPFRNumber abs() const {
93     MPFRNumber result;
94     mpfr_abs(result.value, value, MPFR_RNDN);
95     return result;
96   }
97 
98   MPFRNumber ceil() const {
99     MPFRNumber result;
100     mpfr_ceil(result.value, value);
101     return result;
102   }
103 
104   MPFRNumber cos() const {
105     MPFRNumber result;
106     mpfr_cos(result.value, value, MPFR_RNDN);
107     return result;
108   }
109 
110   MPFRNumber exp() const {
111     MPFRNumber result;
112     mpfr_exp(result.value, value, MPFR_RNDN);
113     return result;
114   }
115 
116   MPFRNumber exp2() const {
117     MPFRNumber result;
118     mpfr_exp2(result.value, value, MPFR_RNDN);
119     return result;
120   }
121 
122   MPFRNumber floor() const {
123     MPFRNumber result;
124     mpfr_floor(result.value, value);
125     return result;
126   }
127 
128   MPFRNumber frexp(int &exp) {
129     MPFRNumber result;
130     mpfr_exp_t resultExp;
131     mpfr_frexp(&resultExp, result.value, value, MPFR_RNDN);
132     exp = resultExp;
133     return result;
134   }
135 
136   MPFRNumber remquo(const MPFRNumber &divisor, int &quotient) {
137     MPFRNumber remainder;
138     long q;
139     mpfr_remquo(remainder.value, &q, value, divisor.value, MPFR_RNDN);
140     quotient = q;
141     return remainder;
142   }
143 
144   MPFRNumber round() const {
145     MPFRNumber result;
146     mpfr_round(result.value, value);
147     return result;
148   }
149 
150   MPFRNumber sin() const {
151     MPFRNumber result;
152     mpfr_sin(result.value, value, MPFR_RNDN);
153     return result;
154   }
155 
156   MPFRNumber sqrt() const {
157     MPFRNumber result;
158     mpfr_sqrt(result.value, value, MPFR_RNDN);
159     return result;
160   }
161 
162   MPFRNumber trunc() const {
163     MPFRNumber result;
164     mpfr_trunc(result.value, value);
165     return result;
166   }
167 
168   std::string str() const {
169     // 200 bytes should be more than sufficient to hold a 100-digit number
170     // plus additional bytes for the decimal point, '-' sign etc.
171     constexpr size_t printBufSize = 200;
172     char buffer[printBufSize];
173     mpfr_snprintf(buffer, printBufSize, "%100.50Rf", value);
174     llvm::StringRef ref(buffer);
175     ref = ref.trim();
176     return ref.str();
177   }
178 
179   // These functions are useful for debugging.
180   template <typename T> T as() const;
181 
182   template <> float as<float>() const { return mpfr_get_flt(value, MPFR_RNDN); }
183   template <> double as<double>() const { return mpfr_get_d(value, MPFR_RNDN); }
184   template <> long double as<long double>() const {
185     return mpfr_get_ld(value, MPFR_RNDN);
186   }
187 
188   void dump(const char *msg) const { mpfr_printf("%s%.128Rf\n", msg, value); }
189 
190   // Return the ULP (units-in-the-last-place) difference between the
191   // stored MPFR and a floating point number.
192   //
193   // We define:
194   //   ULP(mpfr_value, value) = abs(mpfr_value - value) / eps(value)
195   //
196   // Remarks:
197   // 1. ULP < 0.5 will imply that the value is correctly rounded.
198   // 2. We expect that this value and the value to be compared (the [input]
199   //    argument) are reasonable close, and we will provide an upper bound
200   //    of ULP value for testing.  Morever, most of the fractional parts of
201   //    ULP value do not matter much, so using double as the return type
202   //    should be good enough.
203   template <typename T>
204   cpp::EnableIfType<cpp::IsFloatingPointType<T>::Value, double> ulp(T input) {
205     fputil::FPBits<T> bits(input);
206     MPFRNumber mpfrInput(input);
207 
208     // abs(value - input)
209     mpfr_sub(mpfrInput.value, value, mpfrInput.value, MPFR_RNDN);
210     mpfr_abs(mpfrInput.value, mpfrInput.value, MPFR_RNDN);
211 
212     // get eps(input)
213     int epsExponent = bits.exponent - fputil::FPBits<T>::exponentBias -
214                       fputil::MantissaWidth<T>::value;
215     if (bits.exponent == 0) {
216       // correcting denormal exponent
217       ++epsExponent;
218     } else if ((bits.mantissa == 0) && (bits.exponent > 1) &&
219                mpfr_less_p(value, mpfrInput.value)) {
220       // when the input is exactly 2^n, distance (epsilon) between the input
221       // and the next floating point number is different from the distance to
222       // the previous floating point number.  So in that case, if the correct
223       // value from MPFR is smaller than the input, we use the smaller epsilon
224       --epsExponent;
225     }
226 
227     // Since eps(value) is of the form 2^e, instead of dividing such number,
228     // we multiply by its inverse 2^{-e}.
229     mpfr_mul_2si(mpfrInput.value, mpfrInput.value, -epsExponent, MPFR_RNDN);
230 
231     return mpfrInput.as<double>();
232   }
233 };
234 
235 namespace internal {
236 
237 template <typename InputType>
238 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
239 unaryOperation(Operation op, InputType input) {
240   MPFRNumber mpfrInput(input);
241   switch (op) {
242   case Operation::Abs:
243     return mpfrInput.abs();
244   case Operation::Ceil:
245     return mpfrInput.ceil();
246   case Operation::Cos:
247     return mpfrInput.cos();
248   case Operation::Exp:
249     return mpfrInput.exp();
250   case Operation::Exp2:
251     return mpfrInput.exp2();
252   case Operation::Floor:
253     return mpfrInput.floor();
254   case Operation::Round:
255     return mpfrInput.round();
256   case Operation::Sin:
257     return mpfrInput.sin();
258   case Operation::Sqrt:
259     return mpfrInput.sqrt();
260   case Operation::Trunc:
261     return mpfrInput.trunc();
262   default:
263     __builtin_unreachable();
264   }
265 }
266 
267 template <typename InputType>
268 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
269 unaryOperationTwoOutputs(Operation op, InputType input, int &output) {
270   MPFRNumber mpfrInput(input);
271   switch (op) {
272   case Operation::Frexp:
273     return mpfrInput.frexp(output);
274   default:
275     __builtin_unreachable();
276   }
277 }
278 
279 template <typename InputType>
280 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
281 binaryOperationTwoOutputs(Operation op, InputType x, InputType y, int &output) {
282   MPFRNumber inputX(x), inputY(y);
283   switch (op) {
284   case Operation::RemQuo:
285     return inputX.remquo(inputY, output);
286   default:
287     __builtin_unreachable();
288   }
289 }
290 
291 template <typename T>
292 void explainUnaryOperationSingleOutputError(Operation op, T input, T matchValue,
293                                             testutils::StreamWrapper &OS) {
294   MPFRNumber mpfrInput(input);
295   MPFRNumber mpfrResult = unaryOperation(op, input);
296   MPFRNumber mpfrMatchValue(matchValue);
297   FPBits<T> inputBits(input);
298   FPBits<T> matchBits(matchValue);
299   FPBits<T> mpfrResultBits(mpfrResult.as<T>());
300   OS << "Match value not within tolerance value of MPFR result:\n"
301      << "  Input decimal: " << mpfrInput.str() << '\n';
302   __llvm_libc::fputil::testing::describeValue("     Input bits: ", input, OS);
303   OS << '\n' << "  Match decimal: " << mpfrMatchValue.str() << '\n';
304   __llvm_libc::fputil::testing::describeValue("     Match bits: ", matchValue,
305                                               OS);
306   OS << '\n' << "    MPFR result: " << mpfrResult.str() << '\n';
307   __llvm_libc::fputil::testing::describeValue(
308       "   MPFR rounded: ", mpfrResult.as<T>(), OS);
309   OS << '\n';
310   OS << "      ULP error: " << std::to_string(mpfrResult.ulp(matchValue))
311      << '\n';
312 }
313 
314 template void
315 explainUnaryOperationSingleOutputError<float>(Operation op, float, float,
316                                               testutils::StreamWrapper &);
317 template void
318 explainUnaryOperationSingleOutputError<double>(Operation op, double, double,
319                                                testutils::StreamWrapper &);
320 template void explainUnaryOperationSingleOutputError<long double>(
321     Operation op, long double, long double, testutils::StreamWrapper &);
322 
323 template <typename T>
324 void explainUnaryOperationTwoOutputsError(Operation op, T input,
325                                           const BinaryOutput<T> &libcResult,
326                                           testutils::StreamWrapper &OS) {
327   MPFRNumber mpfrInput(input);
328   FPBits<T> inputBits(input);
329   int mpfrIntResult;
330   MPFRNumber mpfrResult = unaryOperationTwoOutputs(op, input, mpfrIntResult);
331 
332   if (mpfrIntResult != libcResult.i) {
333     OS << "MPFR integral result: " << mpfrIntResult << '\n'
334        << "Libc integral result: " << libcResult.i << '\n';
335   } else {
336     OS << "Integral result from libc matches integral result from MPFR.\n";
337   }
338 
339   MPFRNumber mpfrMatchValue(libcResult.f);
340   OS << "Libc floating point result is not within tolerance value of the MPFR "
341      << "result.\n\n";
342 
343   OS << "            Input decimal: " << mpfrInput.str() << "\n\n";
344 
345   OS << "Libc floating point value: " << mpfrMatchValue.str() << '\n';
346   __llvm_libc::fputil::testing::describeValue(
347       " Libc floating point bits: ", libcResult.f, OS);
348   OS << "\n\n";
349 
350   OS << "              MPFR result: " << mpfrResult.str() << '\n';
351   __llvm_libc::fputil::testing::describeValue(
352       "             MPFR rounded: ", mpfrResult.as<T>(), OS);
353   OS << '\n'
354      << "                ULP error: "
355      << std::to_string(mpfrResult.ulp(libcResult.f)) << '\n';
356 }
357 
358 template void explainUnaryOperationTwoOutputsError<float>(
359     Operation, float, const BinaryOutput<float> &, testutils::StreamWrapper &);
360 template void
361 explainUnaryOperationTwoOutputsError<double>(Operation, double,
362                                              const BinaryOutput<double> &,
363                                              testutils::StreamWrapper &);
364 template void explainUnaryOperationTwoOutputsError<long double>(
365     Operation, long double, const BinaryOutput<long double> &,
366     testutils::StreamWrapper &);
367 
368 template <typename T>
369 void explainBinaryOperationTwoOutputsError(Operation op,
370                                            const BinaryInput<T> &input,
371                                            const BinaryOutput<T> &libcResult,
372                                            testutils::StreamWrapper &OS) {
373   MPFRNumber mpfrX(input.x);
374   MPFRNumber mpfrY(input.y);
375   FPBits<T> xbits(input.x);
376   FPBits<T> ybits(input.y);
377   int mpfrIntResult;
378   MPFRNumber mpfrResult =
379       binaryOperationTwoOutputs(op, input.x, input.y, mpfrIntResult);
380   MPFRNumber mpfrMatchValue(libcResult.f);
381 
382   OS << "Input decimal: x: " << mpfrX.str() << " y: " << mpfrY.str() << '\n'
383      << "MPFR integral result: " << mpfrIntResult << '\n'
384      << "Libc integral result: " << libcResult.i << '\n'
385      << "Libc floating point result: " << mpfrMatchValue.str() << '\n'
386      << "               MPFR result: " << mpfrResult.str() << '\n';
387   __llvm_libc::fputil::testing::describeValue(
388       "Libc floating point result bits: ", libcResult.f, OS);
389   __llvm_libc::fputil::testing::describeValue(
390       "              MPFR rounded bits: ", mpfrResult.as<T>(), OS);
391   OS << "ULP error: " << std::to_string(mpfrResult.ulp(libcResult.f)) << '\n';
392 }
393 
394 template void explainBinaryOperationTwoOutputsError<float>(
395     Operation, const BinaryInput<float> &, const BinaryOutput<float> &,
396     testutils::StreamWrapper &);
397 template void explainBinaryOperationTwoOutputsError<double>(
398     Operation, const BinaryInput<double> &, const BinaryOutput<double> &,
399     testutils::StreamWrapper &);
400 template void explainBinaryOperationTwoOutputsError<long double>(
401     Operation, const BinaryInput<long double> &,
402     const BinaryOutput<long double> &, testutils::StreamWrapper &);
403 
404 template <typename T>
405 bool compareUnaryOperationSingleOutput(Operation op, T input, T libcResult,
406                                        double ulpError) {
407   // If the ulp error is exactly 0.5 (i.e a tie), we would check that the result
408   // is rounded to the nearest even.
409   MPFRNumber mpfrResult = unaryOperation(op, input);
410   double ulp = mpfrResult.ulp(libcResult);
411   bool bitsAreEven = ((FPBits<T>(libcResult).bitsAsUInt() & 1) == 0);
412   return (ulp < ulpError) ||
413          ((ulp == ulpError) && ((ulp != 0.5) || bitsAreEven));
414 }
415 
416 template bool compareUnaryOperationSingleOutput<float>(Operation, float, float,
417                                                        double);
418 template bool compareUnaryOperationSingleOutput<double>(Operation, double,
419                                                         double, double);
420 template bool compareUnaryOperationSingleOutput<long double>(Operation,
421                                                              long double,
422                                                              long double,
423                                                              double);
424 
425 template <typename T>
426 bool compareUnaryOperationTwoOutputs(Operation op, T input,
427                                      const BinaryOutput<T> &libcResult,
428                                      double ulpError) {
429   int mpfrIntResult;
430   MPFRNumber mpfrResult = unaryOperationTwoOutputs(op, input, mpfrIntResult);
431   double ulp = mpfrResult.ulp(libcResult.f);
432 
433   if (mpfrIntResult != libcResult.i)
434     return false;
435 
436   bool bitsAreEven = ((FPBits<T>(libcResult.f).bitsAsUInt() & 1) == 0);
437   return (ulp < ulpError) ||
438          ((ulp == ulpError) && ((ulp != 0.5) || bitsAreEven));
439 }
440 
441 template bool
442 compareUnaryOperationTwoOutputs<float>(Operation, float,
443                                        const BinaryOutput<float> &, double);
444 template bool
445 compareUnaryOperationTwoOutputs<double>(Operation, double,
446                                         const BinaryOutput<double> &, double);
447 template bool compareUnaryOperationTwoOutputs<long double>(
448     Operation, long double, const BinaryOutput<long double> &, double);
449 
450 template <typename T>
451 bool compareBinaryOperationTwoOutputs(Operation op, const BinaryInput<T> &input,
452                                       const BinaryOutput<T> &libcResult,
453                                       double ulpError) {
454   int mpfrIntResult;
455   MPFRNumber mpfrResult =
456       binaryOperationTwoOutputs(op, input.x, input.y, mpfrIntResult);
457   double ulp = mpfrResult.ulp(libcResult.f);
458 
459   if (mpfrIntResult != libcResult.i) {
460     if (op == Operation::RemQuo) {
461       if ((0x7 & mpfrIntResult) != (0x7 & libcResult.i))
462         return false;
463     } else {
464       return false;
465     }
466   }
467 
468   bool bitsAreEven = ((FPBits<T>(libcResult.f).bitsAsUInt() & 1) == 0);
469   return (ulp < ulpError) ||
470          ((ulp == ulpError) && ((ulp != 0.5) || bitsAreEven));
471 }
472 
473 template bool
474 compareBinaryOperationTwoOutputs<float>(Operation, const BinaryInput<float> &,
475                                         const BinaryOutput<float> &, double);
476 template bool
477 compareBinaryOperationTwoOutputs<double>(Operation, const BinaryInput<double> &,
478                                          const BinaryOutput<double> &, double);
479 template bool compareBinaryOperationTwoOutputs<long double>(
480     Operation, const BinaryInput<long double> &,
481     const BinaryOutput<long double> &, double);
482 
483 } // namespace internal
484 
485 } // namespace mpfr
486 } // namespace testing
487 } // namespace __llvm_libc
488