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