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 hypot(const MPFRNumber &b) {
137     MPFRNumber result;
138     mpfr_hypot(result.value, value, b.value, MPFR_RNDN);
139     return result;
140   }
141 
142   MPFRNumber remquo(const MPFRNumber &divisor, int &quotient) {
143     MPFRNumber remainder;
144     long q;
145     mpfr_remquo(remainder.value, &q, value, divisor.value, MPFR_RNDN);
146     quotient = q;
147     return remainder;
148   }
149 
150   MPFRNumber round() const {
151     MPFRNumber result;
152     mpfr_round(result.value, value);
153     return result;
154   }
155 
156   MPFRNumber sin() const {
157     MPFRNumber result;
158     mpfr_sin(result.value, value, MPFR_RNDN);
159     return result;
160   }
161 
162   MPFRNumber sqrt() const {
163     MPFRNumber result;
164     mpfr_sqrt(result.value, value, MPFR_RNDN);
165     return result;
166   }
167 
168   MPFRNumber trunc() const {
169     MPFRNumber result;
170     mpfr_trunc(result.value, value);
171     return result;
172   }
173 
174   std::string str() const {
175     // 200 bytes should be more than sufficient to hold a 100-digit number
176     // plus additional bytes for the decimal point, '-' sign etc.
177     constexpr size_t printBufSize = 200;
178     char buffer[printBufSize];
179     mpfr_snprintf(buffer, printBufSize, "%100.50Rf", value);
180     llvm::StringRef ref(buffer);
181     ref = ref.trim();
182     return ref.str();
183   }
184 
185   // These functions are useful for debugging.
186   template <typename T> T as() const;
187 
188   template <> float as<float>() const { return mpfr_get_flt(value, MPFR_RNDN); }
189   template <> double as<double>() const { return mpfr_get_d(value, MPFR_RNDN); }
190   template <> long double as<long double>() const {
191     return mpfr_get_ld(value, MPFR_RNDN);
192   }
193 
194   void dump(const char *msg) const { mpfr_printf("%s%.128Rf\n", msg, value); }
195 
196   // Return the ULP (units-in-the-last-place) difference between the
197   // stored MPFR and a floating point number.
198   //
199   // We define:
200   //   ULP(mpfr_value, value) = abs(mpfr_value - value) / eps(value)
201   //
202   // Remarks:
203   // 1. ULP < 0.5 will imply that the value is correctly rounded.
204   // 2. We expect that this value and the value to be compared (the [input]
205   //    argument) are reasonable close, and we will provide an upper bound
206   //    of ULP value for testing.  Morever, most of the fractional parts of
207   //    ULP value do not matter much, so using double as the return type
208   //    should be good enough.
209   template <typename T>
210   cpp::EnableIfType<cpp::IsFloatingPointType<T>::Value, double> ulp(T input) {
211     fputil::FPBits<T> bits(input);
212     MPFRNumber mpfrInput(input);
213 
214     // abs(value - input)
215     mpfr_sub(mpfrInput.value, value, mpfrInput.value, MPFR_RNDN);
216     mpfr_abs(mpfrInput.value, mpfrInput.value, MPFR_RNDN);
217 
218     // get eps(input)
219     int epsExponent = bits.exponent - fputil::FPBits<T>::exponentBias -
220                       fputil::MantissaWidth<T>::value;
221     if (bits.exponent == 0) {
222       // correcting denormal exponent
223       ++epsExponent;
224     } else if ((bits.mantissa == 0) && (bits.exponent > 1) &&
225                mpfr_less_p(value, mpfrInput.value)) {
226       // when the input is exactly 2^n, distance (epsilon) between the input
227       // and the next floating point number is different from the distance to
228       // the previous floating point number.  So in that case, if the correct
229       // value from MPFR is smaller than the input, we use the smaller epsilon
230       --epsExponent;
231     }
232 
233     // Since eps(value) is of the form 2^e, instead of dividing such number,
234     // we multiply by its inverse 2^{-e}.
235     mpfr_mul_2si(mpfrInput.value, mpfrInput.value, -epsExponent, MPFR_RNDN);
236 
237     return mpfrInput.as<double>();
238   }
239 };
240 
241 namespace internal {
242 
243 template <typename InputType>
244 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
245 unaryOperation(Operation op, InputType input) {
246   MPFRNumber mpfrInput(input);
247   switch (op) {
248   case Operation::Abs:
249     return mpfrInput.abs();
250   case Operation::Ceil:
251     return mpfrInput.ceil();
252   case Operation::Cos:
253     return mpfrInput.cos();
254   case Operation::Exp:
255     return mpfrInput.exp();
256   case Operation::Exp2:
257     return mpfrInput.exp2();
258   case Operation::Floor:
259     return mpfrInput.floor();
260   case Operation::Round:
261     return mpfrInput.round();
262   case Operation::Sin:
263     return mpfrInput.sin();
264   case Operation::Sqrt:
265     return mpfrInput.sqrt();
266   case Operation::Trunc:
267     return mpfrInput.trunc();
268   default:
269     __builtin_unreachable();
270   }
271 }
272 
273 template <typename InputType>
274 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
275 unaryOperationTwoOutputs(Operation op, InputType input, int &output) {
276   MPFRNumber mpfrInput(input);
277   switch (op) {
278   case Operation::Frexp:
279     return mpfrInput.frexp(output);
280   default:
281     __builtin_unreachable();
282   }
283 }
284 
285 template <typename InputType>
286 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
287 binaryOperationOneOutput(Operation op, InputType x, InputType y) {
288   MPFRNumber inputX(x), inputY(y);
289   switch (op) {
290   case Operation::Hypot:
291     return inputX.hypot(inputY);
292   default:
293     __builtin_unreachable();
294   }
295 }
296 
297 template <typename InputType>
298 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
299 binaryOperationTwoOutputs(Operation op, InputType x, InputType y, int &output) {
300   MPFRNumber inputX(x), inputY(y);
301   switch (op) {
302   case Operation::RemQuo:
303     return inputX.remquo(inputY, output);
304   default:
305     __builtin_unreachable();
306   }
307 }
308 
309 template <typename T>
310 void explainUnaryOperationSingleOutputError(Operation op, T input, T matchValue,
311                                             testutils::StreamWrapper &OS) {
312   MPFRNumber mpfrInput(input);
313   MPFRNumber mpfrResult = unaryOperation(op, input);
314   MPFRNumber mpfrMatchValue(matchValue);
315   FPBits<T> inputBits(input);
316   FPBits<T> matchBits(matchValue);
317   FPBits<T> mpfrResultBits(mpfrResult.as<T>());
318   OS << "Match value not within tolerance value of MPFR result:\n"
319      << "  Input decimal: " << mpfrInput.str() << '\n';
320   __llvm_libc::fputil::testing::describeValue("     Input bits: ", input, OS);
321   OS << '\n' << "  Match decimal: " << mpfrMatchValue.str() << '\n';
322   __llvm_libc::fputil::testing::describeValue("     Match bits: ", matchValue,
323                                               OS);
324   OS << '\n' << "    MPFR result: " << mpfrResult.str() << '\n';
325   __llvm_libc::fputil::testing::describeValue(
326       "   MPFR rounded: ", mpfrResult.as<T>(), OS);
327   OS << '\n';
328   OS << "      ULP error: " << std::to_string(mpfrResult.ulp(matchValue))
329      << '\n';
330 }
331 
332 template void
333 explainUnaryOperationSingleOutputError<float>(Operation op, float, float,
334                                               testutils::StreamWrapper &);
335 template void
336 explainUnaryOperationSingleOutputError<double>(Operation op, double, double,
337                                                testutils::StreamWrapper &);
338 template void explainUnaryOperationSingleOutputError<long double>(
339     Operation op, long double, long double, testutils::StreamWrapper &);
340 
341 template <typename T>
342 void explainUnaryOperationTwoOutputsError(Operation op, T input,
343                                           const BinaryOutput<T> &libcResult,
344                                           testutils::StreamWrapper &OS) {
345   MPFRNumber mpfrInput(input);
346   FPBits<T> inputBits(input);
347   int mpfrIntResult;
348   MPFRNumber mpfrResult = unaryOperationTwoOutputs(op, input, mpfrIntResult);
349 
350   if (mpfrIntResult != libcResult.i) {
351     OS << "MPFR integral result: " << mpfrIntResult << '\n'
352        << "Libc integral result: " << libcResult.i << '\n';
353   } else {
354     OS << "Integral result from libc matches integral result from MPFR.\n";
355   }
356 
357   MPFRNumber mpfrMatchValue(libcResult.f);
358   OS << "Libc floating point result is not within tolerance value of the MPFR "
359      << "result.\n\n";
360 
361   OS << "            Input decimal: " << mpfrInput.str() << "\n\n";
362 
363   OS << "Libc floating point value: " << mpfrMatchValue.str() << '\n';
364   __llvm_libc::fputil::testing::describeValue(
365       " Libc floating point bits: ", libcResult.f, OS);
366   OS << "\n\n";
367 
368   OS << "              MPFR result: " << mpfrResult.str() << '\n';
369   __llvm_libc::fputil::testing::describeValue(
370       "             MPFR rounded: ", mpfrResult.as<T>(), OS);
371   OS << '\n'
372      << "                ULP error: "
373      << std::to_string(mpfrResult.ulp(libcResult.f)) << '\n';
374 }
375 
376 template void explainUnaryOperationTwoOutputsError<float>(
377     Operation, float, const BinaryOutput<float> &, testutils::StreamWrapper &);
378 template void
379 explainUnaryOperationTwoOutputsError<double>(Operation, double,
380                                              const BinaryOutput<double> &,
381                                              testutils::StreamWrapper &);
382 template void explainUnaryOperationTwoOutputsError<long double>(
383     Operation, long double, const BinaryOutput<long double> &,
384     testutils::StreamWrapper &);
385 
386 template <typename T>
387 void explainBinaryOperationTwoOutputsError(Operation op,
388                                            const BinaryInput<T> &input,
389                                            const BinaryOutput<T> &libcResult,
390                                            testutils::StreamWrapper &OS) {
391   MPFRNumber mpfrX(input.x);
392   MPFRNumber mpfrY(input.y);
393   FPBits<T> xbits(input.x);
394   FPBits<T> ybits(input.y);
395   int mpfrIntResult;
396   MPFRNumber mpfrResult =
397       binaryOperationTwoOutputs(op, input.x, input.y, mpfrIntResult);
398   MPFRNumber mpfrMatchValue(libcResult.f);
399 
400   OS << "Input decimal: x: " << mpfrX.str() << " y: " << mpfrY.str() << '\n'
401      << "MPFR integral result: " << mpfrIntResult << '\n'
402      << "Libc integral result: " << libcResult.i << '\n'
403      << "Libc floating point result: " << mpfrMatchValue.str() << '\n'
404      << "               MPFR result: " << mpfrResult.str() << '\n';
405   __llvm_libc::fputil::testing::describeValue(
406       "Libc floating point result bits: ", libcResult.f, OS);
407   __llvm_libc::fputil::testing::describeValue(
408       "              MPFR rounded bits: ", mpfrResult.as<T>(), OS);
409   OS << "ULP error: " << std::to_string(mpfrResult.ulp(libcResult.f)) << '\n';
410 }
411 
412 template void explainBinaryOperationTwoOutputsError<float>(
413     Operation, const BinaryInput<float> &, const BinaryOutput<float> &,
414     testutils::StreamWrapper &);
415 template void explainBinaryOperationTwoOutputsError<double>(
416     Operation, const BinaryInput<double> &, const BinaryOutput<double> &,
417     testutils::StreamWrapper &);
418 template void explainBinaryOperationTwoOutputsError<long double>(
419     Operation, const BinaryInput<long double> &,
420     const BinaryOutput<long double> &, testutils::StreamWrapper &);
421 
422 template <typename T>
423 void explainBinaryOperationOneOutputError(Operation op,
424                                           const BinaryInput<T> &input,
425                                           T libcResult,
426                                           testutils::StreamWrapper &OS) {
427   MPFRNumber mpfrX(input.x);
428   MPFRNumber mpfrY(input.y);
429   FPBits<T> xbits(input.x);
430   FPBits<T> ybits(input.y);
431   MPFRNumber mpfrResult = binaryOperationOneOutput(op, input.x, input.y);
432   MPFRNumber mpfrMatchValue(libcResult);
433 
434   OS << "Input decimal: x: " << mpfrX.str() << " y: " << mpfrY.str() << '\n';
435   __llvm_libc::fputil::testing::describeValue("First input bits: ", input.x,
436                                               OS);
437   __llvm_libc::fputil::testing::describeValue("Second input bits: ", input.y,
438                                               OS);
439 
440   OS << "Libc result: " << mpfrMatchValue.str() << '\n'
441      << "MPFR result: " << mpfrResult.str() << '\n';
442   __llvm_libc::fputil::testing::describeValue(
443       "Libc floating point result bits: ", libcResult, OS);
444   __llvm_libc::fputil::testing::describeValue(
445       "              MPFR rounded bits: ", mpfrResult.as<T>(), OS);
446   OS << "ULP error: " << std::to_string(mpfrResult.ulp(libcResult)) << '\n';
447 }
448 
449 template void explainBinaryOperationOneOutputError<float>(
450     Operation, const BinaryInput<float> &, float, testutils::StreamWrapper &);
451 template void explainBinaryOperationOneOutputError<double>(
452     Operation, const BinaryInput<double> &, double, testutils::StreamWrapper &);
453 template void explainBinaryOperationOneOutputError<long double>(
454     Operation, const BinaryInput<long double> &, long double,
455     testutils::StreamWrapper &);
456 
457 template <typename T>
458 bool compareUnaryOperationSingleOutput(Operation op, T input, T libcResult,
459                                        double ulpError) {
460   // If the ulp error is exactly 0.5 (i.e a tie), we would check that the result
461   // is rounded to the nearest even.
462   MPFRNumber mpfrResult = unaryOperation(op, input);
463   double ulp = mpfrResult.ulp(libcResult);
464   bool bitsAreEven = ((FPBits<T>(libcResult).bitsAsUInt() & 1) == 0);
465   return (ulp < ulpError) ||
466          ((ulp == ulpError) && ((ulp != 0.5) || bitsAreEven));
467 }
468 
469 template bool compareUnaryOperationSingleOutput<float>(Operation, float, float,
470                                                        double);
471 template bool compareUnaryOperationSingleOutput<double>(Operation, double,
472                                                         double, double);
473 template bool compareUnaryOperationSingleOutput<long double>(Operation,
474                                                              long double,
475                                                              long double,
476                                                              double);
477 
478 template <typename T>
479 bool compareUnaryOperationTwoOutputs(Operation op, T input,
480                                      const BinaryOutput<T> &libcResult,
481                                      double ulpError) {
482   int mpfrIntResult;
483   MPFRNumber mpfrResult = unaryOperationTwoOutputs(op, input, mpfrIntResult);
484   double ulp = mpfrResult.ulp(libcResult.f);
485 
486   if (mpfrIntResult != libcResult.i)
487     return false;
488 
489   bool bitsAreEven = ((FPBits<T>(libcResult.f).bitsAsUInt() & 1) == 0);
490   return (ulp < ulpError) ||
491          ((ulp == ulpError) && ((ulp != 0.5) || bitsAreEven));
492 }
493 
494 template bool
495 compareUnaryOperationTwoOutputs<float>(Operation, float,
496                                        const BinaryOutput<float> &, double);
497 template bool
498 compareUnaryOperationTwoOutputs<double>(Operation, double,
499                                         const BinaryOutput<double> &, double);
500 template bool compareUnaryOperationTwoOutputs<long double>(
501     Operation, long double, const BinaryOutput<long double> &, double);
502 
503 template <typename T>
504 bool compareBinaryOperationTwoOutputs(Operation op, const BinaryInput<T> &input,
505                                       const BinaryOutput<T> &libcResult,
506                                       double ulpError) {
507   int mpfrIntResult;
508   MPFRNumber mpfrResult =
509       binaryOperationTwoOutputs(op, input.x, input.y, mpfrIntResult);
510   double ulp = mpfrResult.ulp(libcResult.f);
511 
512   if (mpfrIntResult != libcResult.i) {
513     if (op == Operation::RemQuo) {
514       if ((0x7 & mpfrIntResult) != (0x7 & libcResult.i))
515         return false;
516     } else {
517       return false;
518     }
519   }
520 
521   bool bitsAreEven = ((FPBits<T>(libcResult.f).bitsAsUInt() & 1) == 0);
522   return (ulp < ulpError) ||
523          ((ulp == ulpError) && ((ulp != 0.5) || bitsAreEven));
524 }
525 
526 template bool
527 compareBinaryOperationTwoOutputs<float>(Operation, const BinaryInput<float> &,
528                                         const BinaryOutput<float> &, double);
529 template bool
530 compareBinaryOperationTwoOutputs<double>(Operation, const BinaryInput<double> &,
531                                          const BinaryOutput<double> &, double);
532 template bool compareBinaryOperationTwoOutputs<long double>(
533     Operation, const BinaryInput<long double> &,
534     const BinaryOutput<long double> &, double);
535 
536 template <typename T>
537 bool compareBinaryOperationOneOutput(Operation op, const BinaryInput<T> &input,
538                                      T libcResult, double ulpError) {
539   MPFRNumber mpfrResult = binaryOperationOneOutput(op, input.x, input.y);
540   double ulp = mpfrResult.ulp(libcResult);
541 
542   bool bitsAreEven = ((FPBits<T>(libcResult).bitsAsUInt() & 1) == 0);
543   return (ulp < ulpError) ||
544          ((ulp == ulpError) && ((ulp != 0.5) || bitsAreEven));
545 }
546 
547 template bool compareBinaryOperationOneOutput<float>(Operation,
548                                                      const BinaryInput<float> &,
549                                                      float, double);
550 template bool
551 compareBinaryOperationOneOutput<double>(Operation, const BinaryInput<double> &,
552                                         double, double);
553 template bool compareBinaryOperationOneOutput<long double>(
554     Operation, const BinaryInput<long double> &, long double, double);
555 
556 } // namespace internal
557 
558 } // namespace mpfr
559 } // namespace testing
560 } // namespace __llvm_libc
561