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/CPP/StringView.h"
12 #include "utils/FPUtil/FPBits.h"
13 #include "utils/FPUtil/TestHelpers.h"
14 
15 #include <memory>
16 #include <stdint.h>
17 #include <string>
18 
19 #ifdef CUSTOM_MPFR_INCLUDER
20 // Some downstream repos are monoliths carrying MPFR sources in their third
21 // party directory. In such repos, including the MPFR header as
22 // `#include <mpfr.h>` is either disallowed or not possible. If that is the
23 // case, a file named `CustomMPFRIncluder.h` should be added through which the
24 // MPFR header can be included in manner allowed in that repo.
25 #include "CustomMPFRIncluder.h"
26 #else
27 #include <mpfr.h>
28 #endif
29 
30 template <typename T> using FPBits = __llvm_libc::fputil::FPBits<T>;
31 
32 namespace __llvm_libc {
33 namespace testing {
34 namespace mpfr {
35 
36 template <typename T> struct Precision;
37 
38 template <> struct Precision<float> {
39   static constexpr unsigned int value = 24;
40 };
41 
42 template <> struct Precision<double> {
43   static constexpr unsigned int value = 53;
44 };
45 
46 #if !(defined(__x86_64__) || defined(__i386__))
47 template <> struct Precision<long double> {
48   static constexpr unsigned int value = 64;
49 };
50 #else
51 template <> struct Precision<long double> {
52   static constexpr unsigned int value = 113;
53 };
54 #endif
55 
56 class MPFRNumber {
57   // A precision value which allows sufficiently large additional
58   // precision even compared to quad-precision floating point values.
59   unsigned int mpfrPrecision;
60 
61   mpfr_t value;
62 
63 public:
64   MPFRNumber() : mpfrPrecision(128) { mpfr_init2(value, mpfrPrecision); }
65 
66   // We use explicit EnableIf specializations to disallow implicit
67   // conversions. Implicit conversions can potentially lead to loss of
68   // precision.
69   template <typename XType,
70             cpp::EnableIfType<cpp::IsSame<float, XType>::Value, int> = 0>
71   explicit MPFRNumber(XType x, int precision = 128) : mpfrPrecision(precision) {
72     mpfr_init2(value, mpfrPrecision);
73     mpfr_set_flt(value, x, MPFR_RNDN);
74   }
75 
76   template <typename XType,
77             cpp::EnableIfType<cpp::IsSame<double, XType>::Value, int> = 0>
78   explicit MPFRNumber(XType x, int precision = 128) : mpfrPrecision(precision) {
79     mpfr_init2(value, mpfrPrecision);
80     mpfr_set_d(value, x, MPFR_RNDN);
81   }
82 
83   template <typename XType,
84             cpp::EnableIfType<cpp::IsSame<long double, XType>::Value, int> = 0>
85   explicit MPFRNumber(XType x, int precision = 128) : mpfrPrecision(precision) {
86     mpfr_init2(value, mpfrPrecision);
87     mpfr_set_ld(value, x, MPFR_RNDN);
88   }
89 
90   template <typename XType,
91             cpp::EnableIfType<cpp::IsIntegral<XType>::Value, int> = 0>
92   explicit MPFRNumber(XType x, int precision = 128) : mpfrPrecision(precision) {
93     mpfr_init2(value, mpfrPrecision);
94     mpfr_set_sj(value, x, MPFR_RNDN);
95   }
96 
97   MPFRNumber(const MPFRNumber &other) : mpfrPrecision(other.mpfrPrecision) {
98     mpfr_init2(value, mpfrPrecision);
99     mpfr_set(value, other.value, MPFR_RNDN);
100   }
101 
102   ~MPFRNumber() {
103     mpfr_clear(value);
104   }
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 sin() const {
205     MPFRNumber result;
206     mpfr_sin(result.value, value, MPFR_RNDN);
207     return result;
208   }
209 
210   MPFRNumber sqrt() const {
211     MPFRNumber result;
212     mpfr_sqrt(result.value, value, MPFR_RNDN);
213     return result;
214   }
215 
216   MPFRNumber tan() const {
217     MPFRNumber result;
218     mpfr_tan(result.value, value, MPFR_RNDN);
219     return result;
220   }
221 
222   MPFRNumber trunc() const {
223     MPFRNumber result;
224     mpfr_trunc(result.value, value);
225     return result;
226   }
227 
228   MPFRNumber fma(const MPFRNumber &b, const MPFRNumber &c) {
229     MPFRNumber result(*this);
230     mpfr_fma(result.value, value, b.value, c.value, MPFR_RNDN);
231     return result;
232   }
233 
234   std::string str() const {
235     // 200 bytes should be more than sufficient to hold a 100-digit number
236     // plus additional bytes for the decimal point, '-' sign etc.
237     constexpr size_t printBufSize = 200;
238     char buffer[printBufSize];
239     mpfr_snprintf(buffer, printBufSize, "%100.50Rf", value);
240     cpp::StringView view(buffer);
241     view = view.trim(' ');
242     return std::string(view.data());
243   }
244 
245   // These functions are useful for debugging.
246   template <typename T> T as() const;
247 
248   template <> float as<float>() const { return mpfr_get_flt(value, MPFR_RNDN); }
249   template <> double as<double>() const { return mpfr_get_d(value, MPFR_RNDN); }
250   template <> long double as<long double>() const {
251     return mpfr_get_ld(value, MPFR_RNDN);
252   }
253 
254   void dump(const char *msg) const { mpfr_printf("%s%.128Rf\n", msg, value); }
255 
256   // Return the ULP (units-in-the-last-place) difference between the
257   // stored MPFR and a floating point number.
258   //
259   // We define:
260   //   ULP(mpfr_value, value) = abs(mpfr_value - value) / eps(value)
261   //
262   // Remarks:
263   // 1. ULP < 0.5 will imply that the value is correctly rounded.
264   // 2. We expect that this value and the value to be compared (the [input]
265   //    argument) are reasonable close, and we will provide an upper bound
266   //    of ULP value for testing.  Morever, most of the fractional parts of
267   //    ULP value do not matter much, so using double as the return type
268   //    should be good enough.
269   template <typename T>
270   cpp::EnableIfType<cpp::IsFloatingPointType<T>::Value, double> ulp(T input) {
271     fputil::FPBits<T> bits(input);
272     MPFRNumber mpfrInput(input);
273 
274     // abs(value - input)
275     mpfr_sub(mpfrInput.value, value, mpfrInput.value, MPFR_RNDN);
276     mpfr_abs(mpfrInput.value, mpfrInput.value, MPFR_RNDN);
277 
278     // get eps(input)
279     int epsExponent = bits.encoding.exponent - fputil::FPBits<T>::exponentBias -
280                       fputil::MantissaWidth<T>::value;
281     if (bits.encoding.exponent == 0) {
282       // correcting denormal exponent
283       ++epsExponent;
284     } else if ((bits.encoding.mantissa == 0) && (bits.encoding.exponent > 1) &&
285                mpfr_less_p(value, mpfrInput.value)) {
286       // when the input is exactly 2^n, distance (epsilon) between the input
287       // and the next floating point number is different from the distance to
288       // the previous floating point number.  So in that case, if the correct
289       // value from MPFR is smaller than the input, we use the smaller epsilon
290       --epsExponent;
291     }
292 
293     // Since eps(value) is of the form 2^e, instead of dividing such number,
294     // we multiply by its inverse 2^{-e}.
295     mpfr_mul_2si(mpfrInput.value, mpfrInput.value, -epsExponent, MPFR_RNDN);
296 
297     return mpfrInput.as<double>();
298   }
299 };
300 
301 namespace internal {
302 
303 template <typename InputType>
304 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
305 unaryOperation(Operation op, InputType input) {
306   MPFRNumber mpfrInput(input);
307   switch (op) {
308   case Operation::Abs:
309     return mpfrInput.abs();
310   case Operation::Ceil:
311     return mpfrInput.ceil();
312   case Operation::Cos:
313     return mpfrInput.cos();
314   case Operation::Exp:
315     return mpfrInput.exp();
316   case Operation::Exp2:
317     return mpfrInput.exp2();
318   case Operation::Expm1:
319     return mpfrInput.expm1();
320   case Operation::Floor:
321     return mpfrInput.floor();
322   case Operation::Round:
323     return mpfrInput.round();
324   case Operation::Sin:
325     return mpfrInput.sin();
326   case Operation::Sqrt:
327     return mpfrInput.sqrt();
328   case Operation::Tan:
329     return mpfrInput.tan();
330   case Operation::Trunc:
331     return mpfrInput.trunc();
332   default:
333     __builtin_unreachable();
334   }
335 }
336 
337 template <typename InputType>
338 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
339 unaryOperationTwoOutputs(Operation op, InputType input, int &output) {
340   MPFRNumber mpfrInput(input);
341   switch (op) {
342   case Operation::Frexp:
343     return mpfrInput.frexp(output);
344   default:
345     __builtin_unreachable();
346   }
347 }
348 
349 template <typename InputType>
350 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
351 binaryOperationOneOutput(Operation op, InputType x, InputType y) {
352   MPFRNumber inputX(x), inputY(y);
353   switch (op) {
354   case Operation::Hypot:
355     return inputX.hypot(inputY);
356   default:
357     __builtin_unreachable();
358   }
359 }
360 
361 template <typename InputType>
362 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
363 binaryOperationTwoOutputs(Operation op, InputType x, InputType y, int &output) {
364   MPFRNumber inputX(x), inputY(y);
365   switch (op) {
366   case Operation::RemQuo:
367     return inputX.remquo(inputY, output);
368   default:
369     __builtin_unreachable();
370   }
371 }
372 
373 template <typename InputType>
374 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
375 ternaryOperationOneOutput(Operation op, InputType x, InputType y, InputType z) {
376   // For FMA function, we just need to compare with the mpfr_fma with the same
377   // precision as InputType.  Using higher precision as the intermediate results
378   // to compare might incorrectly fail due to double-rounding errors.
379   constexpr unsigned int prec = Precision<InputType>::value;
380   MPFRNumber inputX(x, prec), inputY(y, prec), inputZ(z, prec);
381   switch (op) {
382   case Operation::Fma:
383     return inputX.fma(inputY, inputZ);
384   default:
385     __builtin_unreachable();
386   }
387 }
388 
389 template <typename T>
390 void explainUnaryOperationSingleOutputError(Operation op, T input, T matchValue,
391                                             testutils::StreamWrapper &OS) {
392   MPFRNumber mpfrInput(input);
393   MPFRNumber mpfrResult = unaryOperation(op, input);
394   MPFRNumber mpfrMatchValue(matchValue);
395   FPBits<T> inputBits(input);
396   FPBits<T> matchBits(matchValue);
397   FPBits<T> mpfrResultBits(mpfrResult.as<T>());
398   OS << "Match value not within tolerance value of MPFR result:\n"
399      << "  Input decimal: " << mpfrInput.str() << '\n';
400   __llvm_libc::fputil::testing::describeValue("     Input bits: ", input, OS);
401   OS << '\n' << "  Match decimal: " << mpfrMatchValue.str() << '\n';
402   __llvm_libc::fputil::testing::describeValue("     Match bits: ", matchValue,
403                                               OS);
404   OS << '\n' << "    MPFR result: " << mpfrResult.str() << '\n';
405   __llvm_libc::fputil::testing::describeValue(
406       "   MPFR rounded: ", mpfrResult.as<T>(), OS);
407   OS << '\n';
408   OS << "      ULP error: " << std::to_string(mpfrResult.ulp(matchValue))
409      << '\n';
410 }
411 
412 template void
413 explainUnaryOperationSingleOutputError<float>(Operation op, float, float,
414                                               testutils::StreamWrapper &);
415 template void
416 explainUnaryOperationSingleOutputError<double>(Operation op, double, double,
417                                                testutils::StreamWrapper &);
418 template void explainUnaryOperationSingleOutputError<long double>(
419     Operation op, long double, long double, testutils::StreamWrapper &);
420 
421 template <typename T>
422 void explainUnaryOperationTwoOutputsError(Operation op, T input,
423                                           const BinaryOutput<T> &libcResult,
424                                           testutils::StreamWrapper &OS) {
425   MPFRNumber mpfrInput(input);
426   FPBits<T> inputBits(input);
427   int mpfrIntResult;
428   MPFRNumber mpfrResult = unaryOperationTwoOutputs(op, input, mpfrIntResult);
429 
430   if (mpfrIntResult != libcResult.i) {
431     OS << "MPFR integral result: " << mpfrIntResult << '\n'
432        << "Libc integral result: " << libcResult.i << '\n';
433   } else {
434     OS << "Integral result from libc matches integral result from MPFR.\n";
435   }
436 
437   MPFRNumber mpfrMatchValue(libcResult.f);
438   OS << "Libc floating point result is not within tolerance value of the MPFR "
439      << "result.\n\n";
440 
441   OS << "            Input decimal: " << mpfrInput.str() << "\n\n";
442 
443   OS << "Libc floating point value: " << mpfrMatchValue.str() << '\n';
444   __llvm_libc::fputil::testing::describeValue(
445       " Libc floating point bits: ", libcResult.f, OS);
446   OS << "\n\n";
447 
448   OS << "              MPFR result: " << mpfrResult.str() << '\n';
449   __llvm_libc::fputil::testing::describeValue(
450       "             MPFR rounded: ", mpfrResult.as<T>(), OS);
451   OS << '\n'
452      << "                ULP error: "
453      << std::to_string(mpfrResult.ulp(libcResult.f)) << '\n';
454 }
455 
456 template void explainUnaryOperationTwoOutputsError<float>(
457     Operation, float, const BinaryOutput<float> &, testutils::StreamWrapper &);
458 template void
459 explainUnaryOperationTwoOutputsError<double>(Operation, double,
460                                              const BinaryOutput<double> &,
461                                              testutils::StreamWrapper &);
462 template void explainUnaryOperationTwoOutputsError<long double>(
463     Operation, long double, const BinaryOutput<long double> &,
464     testutils::StreamWrapper &);
465 
466 template <typename T>
467 void explainBinaryOperationTwoOutputsError(Operation op,
468                                            const BinaryInput<T> &input,
469                                            const BinaryOutput<T> &libcResult,
470                                            testutils::StreamWrapper &OS) {
471   MPFRNumber mpfrX(input.x);
472   MPFRNumber mpfrY(input.y);
473   FPBits<T> xbits(input.x);
474   FPBits<T> ybits(input.y);
475   int mpfrIntResult;
476   MPFRNumber mpfrResult =
477       binaryOperationTwoOutputs(op, input.x, input.y, mpfrIntResult);
478   MPFRNumber mpfrMatchValue(libcResult.f);
479 
480   OS << "Input decimal: x: " << mpfrX.str() << " y: " << mpfrY.str() << '\n'
481      << "MPFR integral result: " << mpfrIntResult << '\n'
482      << "Libc integral result: " << libcResult.i << '\n'
483      << "Libc floating point result: " << mpfrMatchValue.str() << '\n'
484      << "               MPFR result: " << mpfrResult.str() << '\n';
485   __llvm_libc::fputil::testing::describeValue(
486       "Libc floating point result bits: ", libcResult.f, OS);
487   __llvm_libc::fputil::testing::describeValue(
488       "              MPFR rounded bits: ", mpfrResult.as<T>(), OS);
489   OS << "ULP error: " << std::to_string(mpfrResult.ulp(libcResult.f)) << '\n';
490 }
491 
492 template void explainBinaryOperationTwoOutputsError<float>(
493     Operation, const BinaryInput<float> &, const BinaryOutput<float> &,
494     testutils::StreamWrapper &);
495 template void explainBinaryOperationTwoOutputsError<double>(
496     Operation, const BinaryInput<double> &, const BinaryOutput<double> &,
497     testutils::StreamWrapper &);
498 template void explainBinaryOperationTwoOutputsError<long double>(
499     Operation, const BinaryInput<long double> &,
500     const BinaryOutput<long double> &, testutils::StreamWrapper &);
501 
502 template <typename T>
503 void explainBinaryOperationOneOutputError(Operation op,
504                                           const BinaryInput<T> &input,
505                                           T libcResult,
506                                           testutils::StreamWrapper &OS) {
507   MPFRNumber mpfrX(input.x);
508   MPFRNumber mpfrY(input.y);
509   FPBits<T> xbits(input.x);
510   FPBits<T> ybits(input.y);
511   MPFRNumber mpfrResult = binaryOperationOneOutput(op, input.x, input.y);
512   MPFRNumber mpfrMatchValue(libcResult);
513 
514   OS << "Input decimal: x: " << mpfrX.str() << " y: " << mpfrY.str() << '\n';
515   __llvm_libc::fputil::testing::describeValue("First input bits: ", input.x,
516                                               OS);
517   __llvm_libc::fputil::testing::describeValue("Second input bits: ", input.y,
518                                               OS);
519 
520   OS << "Libc result: " << mpfrMatchValue.str() << '\n'
521      << "MPFR result: " << mpfrResult.str() << '\n';
522   __llvm_libc::fputil::testing::describeValue(
523       "Libc floating point result bits: ", libcResult, OS);
524   __llvm_libc::fputil::testing::describeValue(
525       "              MPFR rounded bits: ", mpfrResult.as<T>(), OS);
526   OS << "ULP error: " << std::to_string(mpfrResult.ulp(libcResult)) << '\n';
527 }
528 
529 template void explainBinaryOperationOneOutputError<float>(
530     Operation, const BinaryInput<float> &, float, testutils::StreamWrapper &);
531 template void explainBinaryOperationOneOutputError<double>(
532     Operation, const BinaryInput<double> &, double, testutils::StreamWrapper &);
533 template void explainBinaryOperationOneOutputError<long double>(
534     Operation, const BinaryInput<long double> &, long double,
535     testutils::StreamWrapper &);
536 
537 template <typename T>
538 void explainTernaryOperationOneOutputError(Operation op,
539                                            const TernaryInput<T> &input,
540                                            T libcResult,
541                                            testutils::StreamWrapper &OS) {
542   MPFRNumber mpfrX(input.x, Precision<T>::value);
543   MPFRNumber mpfrY(input.y, Precision<T>::value);
544   MPFRNumber mpfrZ(input.z, Precision<T>::value);
545   FPBits<T> xbits(input.x);
546   FPBits<T> ybits(input.y);
547   FPBits<T> zbits(input.z);
548   MPFRNumber mpfrResult =
549       ternaryOperationOneOutput(op, input.x, input.y, input.z);
550   MPFRNumber mpfrMatchValue(libcResult);
551 
552   OS << "Input decimal: x: " << mpfrX.str() << " y: " << mpfrY.str()
553      << " z: " << mpfrZ.str() << '\n';
554   __llvm_libc::fputil::testing::describeValue("First input bits: ", input.x,
555                                               OS);
556   __llvm_libc::fputil::testing::describeValue("Second input bits: ", input.y,
557                                               OS);
558   __llvm_libc::fputil::testing::describeValue("Third input bits: ", input.z,
559                                               OS);
560 
561   OS << "Libc result: " << mpfrMatchValue.str() << '\n'
562      << "MPFR result: " << mpfrResult.str() << '\n';
563   __llvm_libc::fputil::testing::describeValue(
564       "Libc floating point result bits: ", libcResult, OS);
565   __llvm_libc::fputil::testing::describeValue(
566       "              MPFR rounded bits: ", mpfrResult.as<T>(), OS);
567   OS << "ULP error: " << std::to_string(mpfrResult.ulp(libcResult)) << '\n';
568 }
569 
570 template void explainTernaryOperationOneOutputError<float>(
571     Operation, const TernaryInput<float> &, float, testutils::StreamWrapper &);
572 template void explainTernaryOperationOneOutputError<double>(
573     Operation, const TernaryInput<double> &, double,
574     testutils::StreamWrapper &);
575 template void explainTernaryOperationOneOutputError<long double>(
576     Operation, const TernaryInput<long double> &, long double,
577     testutils::StreamWrapper &);
578 
579 template <typename T>
580 bool compareUnaryOperationSingleOutput(Operation op, T input, T libcResult,
581                                        double ulpError) {
582   // If the ulp error is exactly 0.5 (i.e a tie), we would check that the result
583   // is rounded to the nearest even.
584   MPFRNumber mpfrResult = unaryOperation(op, input);
585   double ulp = mpfrResult.ulp(libcResult);
586   bool bitsAreEven = ((FPBits<T>(libcResult).uintval() & 1) == 0);
587   return (ulp < ulpError) ||
588          ((ulp == ulpError) && ((ulp != 0.5) || bitsAreEven));
589 }
590 
591 template bool compareUnaryOperationSingleOutput<float>(Operation, float, float,
592                                                        double);
593 template bool compareUnaryOperationSingleOutput<double>(Operation, double,
594                                                         double, double);
595 template bool compareUnaryOperationSingleOutput<long double>(Operation,
596                                                              long double,
597                                                              long double,
598                                                              double);
599 
600 template <typename T>
601 bool compareUnaryOperationTwoOutputs(Operation op, T input,
602                                      const BinaryOutput<T> &libcResult,
603                                      double ulpError) {
604   int mpfrIntResult;
605   MPFRNumber mpfrResult = unaryOperationTwoOutputs(op, input, mpfrIntResult);
606   double ulp = mpfrResult.ulp(libcResult.f);
607 
608   if (mpfrIntResult != libcResult.i)
609     return false;
610 
611   bool bitsAreEven = ((FPBits<T>(libcResult.f).uintval() & 1) == 0);
612   return (ulp < ulpError) ||
613          ((ulp == ulpError) && ((ulp != 0.5) || bitsAreEven));
614 }
615 
616 template bool
617 compareUnaryOperationTwoOutputs<float>(Operation, float,
618                                        const BinaryOutput<float> &, double);
619 template bool
620 compareUnaryOperationTwoOutputs<double>(Operation, double,
621                                         const BinaryOutput<double> &, double);
622 template bool compareUnaryOperationTwoOutputs<long double>(
623     Operation, long double, const BinaryOutput<long double> &, double);
624 
625 template <typename T>
626 bool compareBinaryOperationTwoOutputs(Operation op, const BinaryInput<T> &input,
627                                       const BinaryOutput<T> &libcResult,
628                                       double ulpError) {
629   int mpfrIntResult;
630   MPFRNumber mpfrResult =
631       binaryOperationTwoOutputs(op, input.x, input.y, mpfrIntResult);
632   double ulp = mpfrResult.ulp(libcResult.f);
633 
634   if (mpfrIntResult != libcResult.i) {
635     if (op == Operation::RemQuo) {
636       if ((0x7 & mpfrIntResult) != (0x7 & libcResult.i))
637         return false;
638     } else {
639       return false;
640     }
641   }
642 
643   bool bitsAreEven = ((FPBits<T>(libcResult.f).uintval() & 1) == 0);
644   return (ulp < ulpError) ||
645          ((ulp == ulpError) && ((ulp != 0.5) || bitsAreEven));
646 }
647 
648 template bool
649 compareBinaryOperationTwoOutputs<float>(Operation, const BinaryInput<float> &,
650                                         const BinaryOutput<float> &, double);
651 template bool
652 compareBinaryOperationTwoOutputs<double>(Operation, const BinaryInput<double> &,
653                                          const BinaryOutput<double> &, double);
654 template bool compareBinaryOperationTwoOutputs<long double>(
655     Operation, const BinaryInput<long double> &,
656     const BinaryOutput<long double> &, double);
657 
658 template <typename T>
659 bool compareBinaryOperationOneOutput(Operation op, const BinaryInput<T> &input,
660                                      T libcResult, double ulpError) {
661   MPFRNumber mpfrResult = binaryOperationOneOutput(op, input.x, input.y);
662   double ulp = mpfrResult.ulp(libcResult);
663 
664   bool bitsAreEven = ((FPBits<T>(libcResult).uintval() & 1) == 0);
665   return (ulp < ulpError) ||
666          ((ulp == ulpError) && ((ulp != 0.5) || bitsAreEven));
667 }
668 
669 template bool compareBinaryOperationOneOutput<float>(Operation,
670                                                      const BinaryInput<float> &,
671                                                      float, double);
672 template bool
673 compareBinaryOperationOneOutput<double>(Operation, const BinaryInput<double> &,
674                                         double, double);
675 template bool compareBinaryOperationOneOutput<long double>(
676     Operation, const BinaryInput<long double> &, long double, double);
677 
678 template <typename T>
679 bool compareTernaryOperationOneOutput(Operation op,
680                                       const TernaryInput<T> &input,
681                                       T libcResult, double ulpError) {
682   MPFRNumber mpfrResult =
683       ternaryOperationOneOutput(op, input.x, input.y, input.z);
684   double ulp = mpfrResult.ulp(libcResult);
685 
686   bool bitsAreEven = ((FPBits<T>(libcResult).uintval() & 1) == 0);
687   return (ulp < ulpError) ||
688          ((ulp == ulpError) && ((ulp != 0.5) || bitsAreEven));
689 }
690 
691 template bool
692 compareTernaryOperationOneOutput<float>(Operation, const TernaryInput<float> &,
693                                         float, double);
694 template bool compareTernaryOperationOneOutput<double>(
695     Operation, const TernaryInput<double> &, double, double);
696 template bool compareTernaryOperationOneOutput<long double>(
697     Operation, const TernaryInput<long double> &, long double, double);
698 
699 static mpfr_rnd_t getMPFRRoundingMode(RoundingMode mode) {
700   switch (mode) {
701   case RoundingMode::Upward:
702     return MPFR_RNDU;
703     break;
704   case RoundingMode::Downward:
705     return MPFR_RNDD;
706     break;
707   case RoundingMode::TowardZero:
708     return MPFR_RNDZ;
709     break;
710   case RoundingMode::Nearest:
711     return MPFR_RNDN;
712     break;
713   }
714 }
715 
716 } // namespace internal
717 
718 template <typename T> bool RoundToLong(T x, long &result) {
719   MPFRNumber mpfr(x);
720   return mpfr.roundToLong(result);
721 }
722 
723 template bool RoundToLong<float>(float, long &);
724 template bool RoundToLong<double>(double, long &);
725 template bool RoundToLong<long double>(long double, long &);
726 
727 template <typename T> bool RoundToLong(T x, RoundingMode mode, long &result) {
728   MPFRNumber mpfr(x);
729   return mpfr.roundToLong(internal::getMPFRRoundingMode(mode), result);
730 }
731 
732 template bool RoundToLong<float>(float, RoundingMode, long &);
733 template bool RoundToLong<double>(double, RoundingMode, long &);
734 template bool RoundToLong<long double>(long double, RoundingMode, long &);
735 
736 template <typename T> T Round(T x, RoundingMode mode) {
737   MPFRNumber mpfr(x);
738   MPFRNumber result = mpfr.rint(internal::getMPFRRoundingMode(mode));
739   return result.as<T>();
740 }
741 
742 template float Round<float>(float, RoundingMode);
743 template double Round<double>(double, RoundingMode);
744 template long double Round<long double>(long double, RoundingMode);
745 
746 } // namespace mpfr
747 } // namespace testing
748 } // namespace __llvm_libc
749