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