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 "src/__support/FPUtil/FPBits.h"
12 #include "src/__support/FPUtil/TestHelpers.h"
13 #include "utils/CPP/StringView.h"
14 
15 #include <cmath>
16 #include <memory>
17 #include <stdint.h>
18 #include <string>
19 
20 #ifdef CUSTOM_MPFR_INCLUDER
21 // Some downstream repos are monoliths carrying MPFR sources in their third
22 // party directory. In such repos, including the MPFR header as
23 // `#include <mpfr.h>` is either disallowed or not possible. If that is the
24 // case, a file named `CustomMPFRIncluder.h` should be added through which the
25 // MPFR header can be included in manner allowed in that repo.
26 #include "CustomMPFRIncluder.h"
27 #else
28 #include <mpfr.h>
29 #endif
30 
31 template <typename T> using FPBits = __llvm_libc::fputil::FPBits<T>;
32 
33 namespace __llvm_libc {
34 namespace testing {
35 namespace mpfr {
36 
37 template <typename T> struct Precision;
38 
39 template <> struct Precision<float> {
40   static constexpr unsigned int value = 24;
41 };
42 
43 template <> struct Precision<double> {
44   static constexpr unsigned int value = 53;
45 };
46 
47 #if !(defined(__x86_64__) || defined(__i386__))
48 template <> struct Precision<long double> {
49   static constexpr unsigned int value = 64;
50 };
51 #else
52 template <> struct Precision<long double> {
53   static constexpr unsigned int value = 113;
54 };
55 #endif
56 
57 class MPFRNumber {
58   // A precision value which allows sufficiently large additional
59   // precision even compared to quad-precision floating point values.
60   unsigned int mpfrPrecision;
61 
62   mpfr_t value;
63 
64 public:
65   MPFRNumber() : mpfrPrecision(256) { mpfr_init2(value, mpfrPrecision); }
66 
67   // We use explicit EnableIf specializations to disallow implicit
68   // conversions. Implicit conversions can potentially lead to loss of
69   // precision.
70   template <typename XType,
71             cpp::EnableIfType<cpp::IsSame<float, XType>::Value, int> = 0>
72   explicit MPFRNumber(XType x, int precision = 128) : mpfrPrecision(precision) {
73     mpfr_init2(value, mpfrPrecision);
74     mpfr_set_flt(value, x, MPFR_RNDN);
75   }
76 
77   template <typename XType,
78             cpp::EnableIfType<cpp::IsSame<double, XType>::Value, int> = 0>
79   explicit MPFRNumber(XType x, int precision = 128) : mpfrPrecision(precision) {
80     mpfr_init2(value, mpfrPrecision);
81     mpfr_set_d(value, x, MPFR_RNDN);
82   }
83 
84   template <typename XType,
85             cpp::EnableIfType<cpp::IsSame<long double, XType>::Value, int> = 0>
86   explicit MPFRNumber(XType x, int precision = 128) : mpfrPrecision(precision) {
87     mpfr_init2(value, mpfrPrecision);
88     mpfr_set_ld(value, x, MPFR_RNDN);
89   }
90 
91   template <typename XType,
92             cpp::EnableIfType<cpp::IsIntegral<XType>::Value, int> = 0>
93   explicit MPFRNumber(XType x, int precision = 128) : mpfrPrecision(precision) {
94     mpfr_init2(value, mpfrPrecision);
95     mpfr_set_sj(value, x, MPFR_RNDN);
96   }
97 
98   MPFRNumber(const MPFRNumber &other) : mpfrPrecision(other.mpfrPrecision) {
99     mpfr_init2(value, mpfrPrecision);
100     mpfr_set(value, other.value, MPFR_RNDN);
101   }
102 
103   ~MPFRNumber() {
104     mpfr_clear(value);
105   }
106 
107   MPFRNumber &operator=(const MPFRNumber &rhs) {
108     mpfrPrecision = rhs.mpfrPrecision;
109     mpfr_set(value, rhs.value, MPFR_RNDN);
110     return *this;
111   }
112 
113   MPFRNumber abs() const {
114     MPFRNumber result;
115     mpfr_abs(result.value, value, MPFR_RNDN);
116     return result;
117   }
118 
119   MPFRNumber ceil() const {
120     MPFRNumber result;
121     mpfr_ceil(result.value, value);
122     return result;
123   }
124 
125   MPFRNumber cos() const {
126     MPFRNumber result;
127     mpfr_cos(result.value, value, MPFR_RNDN);
128     return result;
129   }
130 
131   MPFRNumber exp() const {
132     MPFRNumber result;
133     mpfr_exp(result.value, value, MPFR_RNDN);
134     return result;
135   }
136 
137   MPFRNumber exp2() const {
138     MPFRNumber result;
139     mpfr_exp2(result.value, value, MPFR_RNDN);
140     return result;
141   }
142 
143   MPFRNumber expm1() const {
144     MPFRNumber result;
145     mpfr_expm1(result.value, value, MPFR_RNDN);
146     return result;
147   }
148 
149   MPFRNumber floor() const {
150     MPFRNumber result;
151     mpfr_floor(result.value, value);
152     return result;
153   }
154 
155   MPFRNumber frexp(int &exp) {
156     MPFRNumber result;
157     mpfr_exp_t resultExp;
158     mpfr_frexp(&resultExp, result.value, value, MPFR_RNDN);
159     exp = resultExp;
160     return result;
161   }
162 
163   MPFRNumber hypot(const MPFRNumber &b) {
164     MPFRNumber result;
165     mpfr_hypot(result.value, value, b.value, MPFR_RNDN);
166     return result;
167   }
168 
169   MPFRNumber remquo(const MPFRNumber &divisor, int &quotient) {
170     MPFRNumber remainder;
171     long q;
172     mpfr_remquo(remainder.value, &q, value, divisor.value, MPFR_RNDN);
173     quotient = q;
174     return remainder;
175   }
176 
177   MPFRNumber round() const {
178     MPFRNumber result;
179     mpfr_round(result.value, value);
180     return result;
181   }
182 
183   bool roundToLong(long &result) const {
184     // We first calculate the rounded value. This way, when converting
185     // to long using mpfr_get_si, the rounding direction of MPFR_RNDN
186     // (or any other rounding mode), does not have an influence.
187     MPFRNumber roundedValue = round();
188     mpfr_clear_erangeflag();
189     result = mpfr_get_si(roundedValue.value, MPFR_RNDN);
190     return mpfr_erangeflag_p();
191   }
192 
193   bool roundToLong(mpfr_rnd_t rnd, long &result) const {
194     MPFRNumber rint_result;
195     mpfr_rint(rint_result.value, value, rnd);
196     return rint_result.roundToLong(result);
197   }
198 
199   MPFRNumber rint(mpfr_rnd_t rnd) const {
200     MPFRNumber result;
201     mpfr_rint(result.value, value, rnd);
202     return result;
203   }
204 
205   MPFRNumber mod_2pi() const {
206     MPFRNumber result(0.0, 1280);
207     MPFRNumber _2pi(0.0, 1280);
208     mpfr_const_pi(_2pi.value, MPFR_RNDN);
209     mpfr_mul_si(_2pi.value, _2pi.value, 2, MPFR_RNDN);
210     mpfr_fmod(result.value, value, _2pi.value, MPFR_RNDN);
211     return result;
212   }
213 
214   MPFRNumber mod_pi_over_2() const {
215     MPFRNumber result(0.0, 1280);
216     MPFRNumber pi_over_2(0.0, 1280);
217     mpfr_const_pi(pi_over_2.value, MPFR_RNDN);
218     mpfr_mul_d(pi_over_2.value, pi_over_2.value, 0.5, MPFR_RNDN);
219     mpfr_fmod(result.value, value, pi_over_2.value, MPFR_RNDN);
220     return result;
221   }
222 
223   MPFRNumber mod_pi_over_4() const {
224     MPFRNumber result(0.0, 1280);
225     MPFRNumber pi_over_4(0.0, 1280);
226     mpfr_const_pi(pi_over_4.value, MPFR_RNDN);
227     mpfr_mul_d(pi_over_4.value, pi_over_4.value, 0.25, MPFR_RNDN);
228     mpfr_fmod(result.value, value, pi_over_4.value, MPFR_RNDN);
229     return result;
230   }
231 
232   MPFRNumber sin() const {
233     MPFRNumber result;
234     mpfr_sin(result.value, value, MPFR_RNDN);
235     return result;
236   }
237 
238   MPFRNumber sqrt() const {
239     MPFRNumber result;
240     mpfr_sqrt(result.value, value, MPFR_RNDN);
241     return result;
242   }
243 
244   MPFRNumber tan() const {
245     MPFRNumber result;
246     mpfr_tan(result.value, value, MPFR_RNDN);
247     return result;
248   }
249 
250   MPFRNumber trunc() const {
251     MPFRNumber result;
252     mpfr_trunc(result.value, value);
253     return result;
254   }
255 
256   MPFRNumber fma(const MPFRNumber &b, const MPFRNumber &c) {
257     MPFRNumber result(*this);
258     mpfr_fma(result.value, value, b.value, c.value, MPFR_RNDN);
259     return result;
260   }
261 
262   std::string str() const {
263     // 200 bytes should be more than sufficient to hold a 100-digit number
264     // plus additional bytes for the decimal point, '-' sign etc.
265     constexpr size_t printBufSize = 200;
266     char buffer[printBufSize];
267     mpfr_snprintf(buffer, printBufSize, "%100.50Rf", value);
268     cpp::StringView view(buffer);
269     view = view.trim(' ');
270     return std::string(view.data());
271   }
272 
273   // These functions are useful for debugging.
274   template <typename T> T as() const;
275 
276   template <> float as<float>() const { return mpfr_get_flt(value, MPFR_RNDN); }
277   template <> double as<double>() const { return mpfr_get_d(value, MPFR_RNDN); }
278   template <> long double as<long double>() const {
279     return mpfr_get_ld(value, MPFR_RNDN);
280   }
281 
282   void dump(const char *msg) const { mpfr_printf("%s%.128Rf\n", msg, value); }
283 
284   // Return the ULP (units-in-the-last-place) difference between the
285   // stored MPFR and a floating point number.
286   //
287   // We define ULP difference as follows:
288   //   If exponents of this value and the |input| are same, then:
289   //     ULP(this_value, input) = abs(this_value - input) / eps(input)
290   //   else:
291   //     max = max(abs(this_value), abs(input))
292   //     min = min(abs(this_value), abs(input))
293   //     maxExponent = exponent(max)
294   //     ULP(this_value, input) = (max - 2^maxExponent) / eps(max) +
295   //                              (2^maxExponent - min) / eps(min)
296   //
297   // Remarks:
298   // 1. A ULP of 0.0 will imply that the value is correctly rounded.
299   // 2. We expect that this value and the value to be compared (the [input]
300   //    argument) are reasonable close, and we will provide an upper bound
301   //    of ULP value for testing.  Morever, most of the fractional parts of
302   //    ULP value do not matter much, so using double as the return type
303   //    should be good enough.
304   // 3. For close enough values (values which don't diff in their exponent by
305   //    not more than 1), a ULP difference of N indicates a bit distance
306   //    of N between this number and [input].
307   // 4. A values of +0.0 and -0.0 are treated as equal.
308   template <typename T>
309   cpp::EnableIfType<cpp::IsFloatingPointType<T>::Value, double> ulp(T input) {
310     T thisAsT = as<T>();
311     if (thisAsT == input)
312       return T(0.0);
313 
314     int thisExponent = fputil::FPBits<T>(thisAsT).getExponent();
315     int inputExponent = fputil::FPBits<T>(input).getExponent();
316     if (thisAsT * input < 0 || thisExponent == inputExponent) {
317       MPFRNumber inputMPFR(input);
318       mpfr_sub(inputMPFR.value, value, inputMPFR.value, MPFR_RNDN);
319       mpfr_abs(inputMPFR.value, inputMPFR.value, MPFR_RNDN);
320       mpfr_mul_2si(inputMPFR.value, inputMPFR.value, -thisExponent, MPFR_RNDN);
321       return inputMPFR.as<double>();
322     }
323 
324     // If the control reaches here, it means that this number and input are
325     // of the same sign but different exponent. In such a case, ULP error is
326     // calculated as sum of two parts.
327     thisAsT = std::abs(thisAsT);
328     input = std::abs(input);
329     T min = thisAsT > input ? input : thisAsT;
330     T max = thisAsT > input ? thisAsT : input;
331     int minExponent = fputil::FPBits<T>(min).getExponent();
332     int maxExponent = fputil::FPBits<T>(max).getExponent();
333 
334     MPFRNumber minMPFR(min);
335     MPFRNumber maxMPFR(max);
336 
337     MPFRNumber pivot(uint32_t(1));
338     mpfr_mul_2si(pivot.value, pivot.value, maxExponent, MPFR_RNDN);
339 
340     mpfr_sub(minMPFR.value, pivot.value, minMPFR.value, MPFR_RNDN);
341     mpfr_mul_2si(minMPFR.value, minMPFR.value, -minExponent, MPFR_RNDN);
342 
343     mpfr_sub(maxMPFR.value, maxMPFR.value, pivot.value, MPFR_RNDN);
344     mpfr_mul_2si(maxMPFR.value, maxMPFR.value, -maxExponent, MPFR_RNDN);
345 
346     mpfr_add(minMPFR.value, minMPFR.value, maxMPFR.value, MPFR_RNDN);
347     return minMPFR.as<double>();
348   }
349 };
350 
351 namespace internal {
352 
353 template <typename InputType>
354 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
355 unaryOperation(Operation op, InputType input) {
356   MPFRNumber mpfrInput(input);
357   switch (op) {
358   case Operation::Abs:
359     return mpfrInput.abs();
360   case Operation::Ceil:
361     return mpfrInput.ceil();
362   case Operation::Cos:
363     return mpfrInput.cos();
364   case Operation::Exp:
365     return mpfrInput.exp();
366   case Operation::Exp2:
367     return mpfrInput.exp2();
368   case Operation::Expm1:
369     return mpfrInput.expm1();
370   case Operation::Floor:
371     return mpfrInput.floor();
372   case Operation::Mod2PI:
373     return mpfrInput.mod_2pi();
374   case Operation::ModPIOver2:
375     return mpfrInput.mod_pi_over_2();
376   case Operation::ModPIOver4:
377     return mpfrInput.mod_pi_over_4();
378   case Operation::Round:
379     return mpfrInput.round();
380   case Operation::Sin:
381     return mpfrInput.sin();
382   case Operation::Sqrt:
383     return mpfrInput.sqrt();
384   case Operation::Tan:
385     return mpfrInput.tan();
386   case Operation::Trunc:
387     return mpfrInput.trunc();
388   default:
389     __builtin_unreachable();
390   }
391 }
392 
393 template <typename InputType>
394 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
395 unaryOperationTwoOutputs(Operation op, InputType input, int &output) {
396   MPFRNumber mpfrInput(input);
397   switch (op) {
398   case Operation::Frexp:
399     return mpfrInput.frexp(output);
400   default:
401     __builtin_unreachable();
402   }
403 }
404 
405 template <typename InputType>
406 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
407 binaryOperationOneOutput(Operation op, InputType x, InputType y) {
408   MPFRNumber inputX(x), inputY(y);
409   switch (op) {
410   case Operation::Hypot:
411     return inputX.hypot(inputY);
412   default:
413     __builtin_unreachable();
414   }
415 }
416 
417 template <typename InputType>
418 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
419 binaryOperationTwoOutputs(Operation op, InputType x, InputType y, int &output) {
420   MPFRNumber inputX(x), inputY(y);
421   switch (op) {
422   case Operation::RemQuo:
423     return inputX.remquo(inputY, output);
424   default:
425     __builtin_unreachable();
426   }
427 }
428 
429 template <typename InputType>
430 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
431 ternaryOperationOneOutput(Operation op, InputType x, InputType y, InputType z) {
432   // For FMA function, we just need to compare with the mpfr_fma with the same
433   // precision as InputType.  Using higher precision as the intermediate results
434   // to compare might incorrectly fail due to double-rounding errors.
435   constexpr unsigned int prec = Precision<InputType>::value;
436   MPFRNumber inputX(x, prec), inputY(y, prec), inputZ(z, prec);
437   switch (op) {
438   case Operation::Fma:
439     return inputX.fma(inputY, inputZ);
440   default:
441     __builtin_unreachable();
442   }
443 }
444 
445 template <typename T>
446 void explainUnaryOperationSingleOutputError(Operation op, T input, T matchValue,
447                                             testutils::StreamWrapper &OS) {
448   MPFRNumber mpfrInput(input);
449   MPFRNumber mpfrResult = unaryOperation(op, input);
450   MPFRNumber mpfrMatchValue(matchValue);
451   FPBits<T> inputBits(input);
452   FPBits<T> matchBits(matchValue);
453   FPBits<T> mpfrResultBits(mpfrResult.as<T>());
454   OS << "Match value not within tolerance value of MPFR result:\n"
455      << "  Input decimal: " << mpfrInput.str() << '\n';
456   __llvm_libc::fputil::testing::describeValue("     Input bits: ", input, OS);
457   OS << '\n' << "  Match decimal: " << mpfrMatchValue.str() << '\n';
458   __llvm_libc::fputil::testing::describeValue("     Match bits: ", matchValue,
459                                               OS);
460   OS << '\n' << "    MPFR result: " << mpfrResult.str() << '\n';
461   __llvm_libc::fputil::testing::describeValue(
462       "   MPFR rounded: ", mpfrResult.as<T>(), OS);
463   OS << '\n';
464   OS << "      ULP error: " << std::to_string(mpfrResult.ulp(matchValue))
465      << '\n';
466 }
467 
468 template void
469 explainUnaryOperationSingleOutputError<float>(Operation op, float, float,
470                                               testutils::StreamWrapper &);
471 template void
472 explainUnaryOperationSingleOutputError<double>(Operation op, double, double,
473                                                testutils::StreamWrapper &);
474 template void explainUnaryOperationSingleOutputError<long double>(
475     Operation op, long double, long double, testutils::StreamWrapper &);
476 
477 template <typename T>
478 void explainUnaryOperationTwoOutputsError(Operation op, T input,
479                                           const BinaryOutput<T> &libcResult,
480                                           testutils::StreamWrapper &OS) {
481   MPFRNumber mpfrInput(input);
482   FPBits<T> inputBits(input);
483   int mpfrIntResult;
484   MPFRNumber mpfrResult = unaryOperationTwoOutputs(op, input, mpfrIntResult);
485 
486   if (mpfrIntResult != libcResult.i) {
487     OS << "MPFR integral result: " << mpfrIntResult << '\n'
488        << "Libc integral result: " << libcResult.i << '\n';
489   } else {
490     OS << "Integral result from libc matches integral result from MPFR.\n";
491   }
492 
493   MPFRNumber mpfrMatchValue(libcResult.f);
494   OS << "Libc floating point result is not within tolerance value of the MPFR "
495      << "result.\n\n";
496 
497   OS << "            Input decimal: " << mpfrInput.str() << "\n\n";
498 
499   OS << "Libc floating point value: " << mpfrMatchValue.str() << '\n';
500   __llvm_libc::fputil::testing::describeValue(
501       " Libc floating point bits: ", libcResult.f, OS);
502   OS << "\n\n";
503 
504   OS << "              MPFR result: " << mpfrResult.str() << '\n';
505   __llvm_libc::fputil::testing::describeValue(
506       "             MPFR rounded: ", mpfrResult.as<T>(), OS);
507   OS << '\n'
508      << "                ULP error: "
509      << std::to_string(mpfrResult.ulp(libcResult.f)) << '\n';
510 }
511 
512 template void explainUnaryOperationTwoOutputsError<float>(
513     Operation, float, const BinaryOutput<float> &, testutils::StreamWrapper &);
514 template void
515 explainUnaryOperationTwoOutputsError<double>(Operation, double,
516                                              const BinaryOutput<double> &,
517                                              testutils::StreamWrapper &);
518 template void explainUnaryOperationTwoOutputsError<long double>(
519     Operation, long double, const BinaryOutput<long double> &,
520     testutils::StreamWrapper &);
521 
522 template <typename T>
523 void explainBinaryOperationTwoOutputsError(Operation op,
524                                            const BinaryInput<T> &input,
525                                            const BinaryOutput<T> &libcResult,
526                                            testutils::StreamWrapper &OS) {
527   MPFRNumber mpfrX(input.x);
528   MPFRNumber mpfrY(input.y);
529   FPBits<T> xbits(input.x);
530   FPBits<T> ybits(input.y);
531   int mpfrIntResult;
532   MPFRNumber mpfrResult =
533       binaryOperationTwoOutputs(op, input.x, input.y, mpfrIntResult);
534   MPFRNumber mpfrMatchValue(libcResult.f);
535 
536   OS << "Input decimal: x: " << mpfrX.str() << " y: " << mpfrY.str() << '\n'
537      << "MPFR integral result: " << mpfrIntResult << '\n'
538      << "Libc integral result: " << libcResult.i << '\n'
539      << "Libc floating point result: " << mpfrMatchValue.str() << '\n'
540      << "               MPFR result: " << mpfrResult.str() << '\n';
541   __llvm_libc::fputil::testing::describeValue(
542       "Libc floating point result bits: ", libcResult.f, OS);
543   __llvm_libc::fputil::testing::describeValue(
544       "              MPFR rounded bits: ", mpfrResult.as<T>(), OS);
545   OS << "ULP error: " << std::to_string(mpfrResult.ulp(libcResult.f)) << '\n';
546 }
547 
548 template void explainBinaryOperationTwoOutputsError<float>(
549     Operation, const BinaryInput<float> &, const BinaryOutput<float> &,
550     testutils::StreamWrapper &);
551 template void explainBinaryOperationTwoOutputsError<double>(
552     Operation, const BinaryInput<double> &, const BinaryOutput<double> &,
553     testutils::StreamWrapper &);
554 template void explainBinaryOperationTwoOutputsError<long double>(
555     Operation, const BinaryInput<long double> &,
556     const BinaryOutput<long double> &, testutils::StreamWrapper &);
557 
558 template <typename T>
559 void explainBinaryOperationOneOutputError(Operation op,
560                                           const BinaryInput<T> &input,
561                                           T libcResult,
562                                           testutils::StreamWrapper &OS) {
563   MPFRNumber mpfrX(input.x);
564   MPFRNumber mpfrY(input.y);
565   FPBits<T> xbits(input.x);
566   FPBits<T> ybits(input.y);
567   MPFRNumber mpfrResult = binaryOperationOneOutput(op, input.x, input.y);
568   MPFRNumber mpfrMatchValue(libcResult);
569 
570   OS << "Input decimal: x: " << mpfrX.str() << " y: " << mpfrY.str() << '\n';
571   __llvm_libc::fputil::testing::describeValue("First input bits: ", input.x,
572                                               OS);
573   __llvm_libc::fputil::testing::describeValue("Second input bits: ", input.y,
574                                               OS);
575 
576   OS << "Libc result: " << mpfrMatchValue.str() << '\n'
577      << "MPFR result: " << mpfrResult.str() << '\n';
578   __llvm_libc::fputil::testing::describeValue(
579       "Libc floating point result bits: ", libcResult, OS);
580   __llvm_libc::fputil::testing::describeValue(
581       "              MPFR rounded bits: ", mpfrResult.as<T>(), OS);
582   OS << "ULP error: " << std::to_string(mpfrResult.ulp(libcResult)) << '\n';
583 }
584 
585 template void explainBinaryOperationOneOutputError<float>(
586     Operation, const BinaryInput<float> &, float, testutils::StreamWrapper &);
587 template void explainBinaryOperationOneOutputError<double>(
588     Operation, const BinaryInput<double> &, double, testutils::StreamWrapper &);
589 template void explainBinaryOperationOneOutputError<long double>(
590     Operation, const BinaryInput<long double> &, long double,
591     testutils::StreamWrapper &);
592 
593 template <typename T>
594 void explainTernaryOperationOneOutputError(Operation op,
595                                            const TernaryInput<T> &input,
596                                            T libcResult,
597                                            testutils::StreamWrapper &OS) {
598   MPFRNumber mpfrX(input.x, Precision<T>::value);
599   MPFRNumber mpfrY(input.y, Precision<T>::value);
600   MPFRNumber mpfrZ(input.z, Precision<T>::value);
601   FPBits<T> xbits(input.x);
602   FPBits<T> ybits(input.y);
603   FPBits<T> zbits(input.z);
604   MPFRNumber mpfrResult =
605       ternaryOperationOneOutput(op, input.x, input.y, input.z);
606   MPFRNumber mpfrMatchValue(libcResult);
607 
608   OS << "Input decimal: x: " << mpfrX.str() << " y: " << mpfrY.str()
609      << " z: " << mpfrZ.str() << '\n';
610   __llvm_libc::fputil::testing::describeValue("First input bits: ", input.x,
611                                               OS);
612   __llvm_libc::fputil::testing::describeValue("Second input bits: ", input.y,
613                                               OS);
614   __llvm_libc::fputil::testing::describeValue("Third input bits: ", input.z,
615                                               OS);
616 
617   OS << "Libc result: " << mpfrMatchValue.str() << '\n'
618      << "MPFR result: " << mpfrResult.str() << '\n';
619   __llvm_libc::fputil::testing::describeValue(
620       "Libc floating point result bits: ", libcResult, OS);
621   __llvm_libc::fputil::testing::describeValue(
622       "              MPFR rounded bits: ", mpfrResult.as<T>(), OS);
623   OS << "ULP error: " << std::to_string(mpfrResult.ulp(libcResult)) << '\n';
624 }
625 
626 template void explainTernaryOperationOneOutputError<float>(
627     Operation, const TernaryInput<float> &, float, testutils::StreamWrapper &);
628 template void explainTernaryOperationOneOutputError<double>(
629     Operation, const TernaryInput<double> &, double,
630     testutils::StreamWrapper &);
631 template void explainTernaryOperationOneOutputError<long double>(
632     Operation, const TernaryInput<long double> &, long double,
633     testutils::StreamWrapper &);
634 
635 template <typename T>
636 bool compareUnaryOperationSingleOutput(Operation op, T input, T libcResult,
637                                        double ulpError) {
638   // If the ulp error is exactly 0.5 (i.e a tie), we would check that the result
639   // is rounded to the nearest even.
640   MPFRNumber mpfrResult = unaryOperation(op, input);
641   double ulp = mpfrResult.ulp(libcResult);
642   bool bitsAreEven = ((FPBits<T>(libcResult).uintval() & 1) == 0);
643   return (ulp < ulpError) ||
644          ((ulp == ulpError) && ((ulp != 0.5) || bitsAreEven));
645 }
646 
647 template bool compareUnaryOperationSingleOutput<float>(Operation, float, float,
648                                                        double);
649 template bool compareUnaryOperationSingleOutput<double>(Operation, double,
650                                                         double, double);
651 template bool compareUnaryOperationSingleOutput<long double>(Operation,
652                                                              long double,
653                                                              long double,
654                                                              double);
655 
656 template <typename T>
657 bool compareUnaryOperationTwoOutputs(Operation op, T input,
658                                      const BinaryOutput<T> &libcResult,
659                                      double ulpError) {
660   int mpfrIntResult;
661   MPFRNumber mpfrResult = unaryOperationTwoOutputs(op, input, mpfrIntResult);
662   double ulp = mpfrResult.ulp(libcResult.f);
663 
664   if (mpfrIntResult != libcResult.i)
665     return false;
666 
667   bool bitsAreEven = ((FPBits<T>(libcResult.f).uintval() & 1) == 0);
668   return (ulp < ulpError) ||
669          ((ulp == ulpError) && ((ulp != 0.5) || bitsAreEven));
670 }
671 
672 template bool
673 compareUnaryOperationTwoOutputs<float>(Operation, float,
674                                        const BinaryOutput<float> &, double);
675 template bool
676 compareUnaryOperationTwoOutputs<double>(Operation, double,
677                                         const BinaryOutput<double> &, double);
678 template bool compareUnaryOperationTwoOutputs<long double>(
679     Operation, long double, const BinaryOutput<long double> &, double);
680 
681 template <typename T>
682 bool compareBinaryOperationTwoOutputs(Operation op, const BinaryInput<T> &input,
683                                       const BinaryOutput<T> &libcResult,
684                                       double ulpError) {
685   int mpfrIntResult;
686   MPFRNumber mpfrResult =
687       binaryOperationTwoOutputs(op, input.x, input.y, mpfrIntResult);
688   double ulp = mpfrResult.ulp(libcResult.f);
689 
690   if (mpfrIntResult != libcResult.i) {
691     if (op == Operation::RemQuo) {
692       if ((0x7 & mpfrIntResult) != (0x7 & libcResult.i))
693         return false;
694     } else {
695       return false;
696     }
697   }
698 
699   bool bitsAreEven = ((FPBits<T>(libcResult.f).uintval() & 1) == 0);
700   return (ulp < ulpError) ||
701          ((ulp == ulpError) && ((ulp != 0.5) || bitsAreEven));
702 }
703 
704 template bool
705 compareBinaryOperationTwoOutputs<float>(Operation, const BinaryInput<float> &,
706                                         const BinaryOutput<float> &, double);
707 template bool
708 compareBinaryOperationTwoOutputs<double>(Operation, const BinaryInput<double> &,
709                                          const BinaryOutput<double> &, double);
710 template bool compareBinaryOperationTwoOutputs<long double>(
711     Operation, const BinaryInput<long double> &,
712     const BinaryOutput<long double> &, double);
713 
714 template <typename T>
715 bool compareBinaryOperationOneOutput(Operation op, const BinaryInput<T> &input,
716                                      T libcResult, double ulpError) {
717   MPFRNumber mpfrResult = binaryOperationOneOutput(op, input.x, input.y);
718   double ulp = mpfrResult.ulp(libcResult);
719 
720   bool bitsAreEven = ((FPBits<T>(libcResult).uintval() & 1) == 0);
721   return (ulp < ulpError) ||
722          ((ulp == ulpError) && ((ulp != 0.5) || bitsAreEven));
723 }
724 
725 template bool compareBinaryOperationOneOutput<float>(Operation,
726                                                      const BinaryInput<float> &,
727                                                      float, double);
728 template bool
729 compareBinaryOperationOneOutput<double>(Operation, const BinaryInput<double> &,
730                                         double, double);
731 template bool compareBinaryOperationOneOutput<long double>(
732     Operation, const BinaryInput<long double> &, long double, double);
733 
734 template <typename T>
735 bool compareTernaryOperationOneOutput(Operation op,
736                                       const TernaryInput<T> &input,
737                                       T libcResult, double ulpError) {
738   MPFRNumber mpfrResult =
739       ternaryOperationOneOutput(op, input.x, input.y, input.z);
740   double ulp = mpfrResult.ulp(libcResult);
741 
742   bool bitsAreEven = ((FPBits<T>(libcResult).uintval() & 1) == 0);
743   return (ulp < ulpError) ||
744          ((ulp == ulpError) && ((ulp != 0.5) || bitsAreEven));
745 }
746 
747 template bool
748 compareTernaryOperationOneOutput<float>(Operation, const TernaryInput<float> &,
749                                         float, double);
750 template bool compareTernaryOperationOneOutput<double>(
751     Operation, const TernaryInput<double> &, double, double);
752 template bool compareTernaryOperationOneOutput<long double>(
753     Operation, const TernaryInput<long double> &, long double, double);
754 
755 static mpfr_rnd_t getMPFRRoundingMode(RoundingMode mode) {
756   switch (mode) {
757   case RoundingMode::Upward:
758     return MPFR_RNDU;
759     break;
760   case RoundingMode::Downward:
761     return MPFR_RNDD;
762     break;
763   case RoundingMode::TowardZero:
764     return MPFR_RNDZ;
765     break;
766   case RoundingMode::Nearest:
767     return MPFR_RNDN;
768     break;
769   }
770 }
771 
772 } // namespace internal
773 
774 template <typename T> bool RoundToLong(T x, long &result) {
775   MPFRNumber mpfr(x);
776   return mpfr.roundToLong(result);
777 }
778 
779 template bool RoundToLong<float>(float, long &);
780 template bool RoundToLong<double>(double, long &);
781 template bool RoundToLong<long double>(long double, long &);
782 
783 template <typename T> bool RoundToLong(T x, RoundingMode mode, long &result) {
784   MPFRNumber mpfr(x);
785   return mpfr.roundToLong(internal::getMPFRRoundingMode(mode), result);
786 }
787 
788 template bool RoundToLong<float>(float, RoundingMode, long &);
789 template bool RoundToLong<double>(double, RoundingMode, long &);
790 template bool RoundToLong<long double>(long double, RoundingMode, long &);
791 
792 template <typename T> T Round(T x, RoundingMode mode) {
793   MPFRNumber mpfr(x);
794   MPFRNumber result = mpfr.rint(internal::getMPFRRoundingMode(mode));
795   return result.as<T>();
796 }
797 
798 template float Round<float>(float, RoundingMode);
799 template double Round<double>(double, RoundingMode);
800 template long double Round<long double>(long double, RoundingMode);
801 
802 } // namespace mpfr
803 } // namespace testing
804 } // namespace __llvm_libc
805