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 floor() const {
143     MPFRNumber result;
144     mpfr_floor(result.value, value);
145     return result;
146   }
147 
148   MPFRNumber frexp(int &exp) {
149     MPFRNumber result;
150     mpfr_exp_t resultExp;
151     mpfr_frexp(&resultExp, result.value, value, MPFR_RNDN);
152     exp = resultExp;
153     return result;
154   }
155 
156   MPFRNumber hypot(const MPFRNumber &b) {
157     MPFRNumber result;
158     mpfr_hypot(result.value, value, b.value, MPFR_RNDN);
159     return result;
160   }
161 
162   MPFRNumber remquo(const MPFRNumber &divisor, int &quotient) {
163     MPFRNumber remainder;
164     long q;
165     mpfr_remquo(remainder.value, &q, value, divisor.value, MPFR_RNDN);
166     quotient = q;
167     return remainder;
168   }
169 
170   MPFRNumber round() const {
171     MPFRNumber result;
172     mpfr_round(result.value, value);
173     return result;
174   }
175 
176   bool roundToLong(long &result) const {
177     // We first calculate the rounded value. This way, when converting
178     // to long using mpfr_get_si, the rounding direction of MPFR_RNDN
179     // (or any other rounding mode), does not have an influence.
180     MPFRNumber roundedValue = round();
181     mpfr_clear_erangeflag();
182     result = mpfr_get_si(roundedValue.value, MPFR_RNDN);
183     return mpfr_erangeflag_p();
184   }
185 
186   bool roundToLong(mpfr_rnd_t rnd, long &result) const {
187     MPFRNumber rint_result;
188     mpfr_rint(rint_result.value, value, rnd);
189     return rint_result.roundToLong(result);
190   }
191 
192   MPFRNumber rint(mpfr_rnd_t rnd) const {
193     MPFRNumber result;
194     mpfr_rint(result.value, value, rnd);
195     return result;
196   }
197 
198   MPFRNumber sin() const {
199     MPFRNumber result;
200     mpfr_sin(result.value, value, MPFR_RNDN);
201     return result;
202   }
203 
204   MPFRNumber sqrt() const {
205     MPFRNumber result;
206     mpfr_sqrt(result.value, value, MPFR_RNDN);
207     return result;
208   }
209 
210   MPFRNumber tan() const {
211     MPFRNumber result;
212     mpfr_tan(result.value, value, MPFR_RNDN);
213     return result;
214   }
215 
216   MPFRNumber trunc() const {
217     MPFRNumber result;
218     mpfr_trunc(result.value, value);
219     return result;
220   }
221 
222   MPFRNumber fma(const MPFRNumber &b, const MPFRNumber &c) {
223     MPFRNumber result(*this);
224     mpfr_fma(result.value, value, b.value, c.value, MPFR_RNDN);
225     return result;
226   }
227 
228   std::string str() const {
229     // 200 bytes should be more than sufficient to hold a 100-digit number
230     // plus additional bytes for the decimal point, '-' sign etc.
231     constexpr size_t printBufSize = 200;
232     char buffer[printBufSize];
233     mpfr_snprintf(buffer, printBufSize, "%100.50Rf", value);
234     cpp::StringView view(buffer);
235     view = view.trim(' ');
236     return std::string(view.data());
237   }
238 
239   // These functions are useful for debugging.
240   template <typename T> T as() const;
241 
242   template <> float as<float>() const { return mpfr_get_flt(value, MPFR_RNDN); }
243   template <> double as<double>() const { return mpfr_get_d(value, MPFR_RNDN); }
244   template <> long double as<long double>() const {
245     return mpfr_get_ld(value, MPFR_RNDN);
246   }
247 
248   void dump(const char *msg) const { mpfr_printf("%s%.128Rf\n", msg, value); }
249 
250   // Return the ULP (units-in-the-last-place) difference between the
251   // stored MPFR and a floating point number.
252   //
253   // We define:
254   //   ULP(mpfr_value, value) = abs(mpfr_value - value) / eps(value)
255   //
256   // Remarks:
257   // 1. ULP < 0.5 will imply that the value is correctly rounded.
258   // 2. We expect that this value and the value to be compared (the [input]
259   //    argument) are reasonable close, and we will provide an upper bound
260   //    of ULP value for testing.  Morever, most of the fractional parts of
261   //    ULP value do not matter much, so using double as the return type
262   //    should be good enough.
263   template <typename T>
264   cpp::EnableIfType<cpp::IsFloatingPointType<T>::Value, double> ulp(T input) {
265     fputil::FPBits<T> bits(input);
266     MPFRNumber mpfrInput(input);
267 
268     // abs(value - input)
269     mpfr_sub(mpfrInput.value, value, mpfrInput.value, MPFR_RNDN);
270     mpfr_abs(mpfrInput.value, mpfrInput.value, MPFR_RNDN);
271 
272     // get eps(input)
273     int epsExponent = bits.encoding.exponent - fputil::FPBits<T>::exponentBias -
274                       fputil::MantissaWidth<T>::value;
275     if (bits.encoding.exponent == 0) {
276       // correcting denormal exponent
277       ++epsExponent;
278     } else if ((bits.encoding.mantissa == 0) && (bits.encoding.exponent > 1) &&
279                mpfr_less_p(value, mpfrInput.value)) {
280       // when the input is exactly 2^n, distance (epsilon) between the input
281       // and the next floating point number is different from the distance to
282       // the previous floating point number.  So in that case, if the correct
283       // value from MPFR is smaller than the input, we use the smaller epsilon
284       --epsExponent;
285     }
286 
287     // Since eps(value) is of the form 2^e, instead of dividing such number,
288     // we multiply by its inverse 2^{-e}.
289     mpfr_mul_2si(mpfrInput.value, mpfrInput.value, -epsExponent, MPFR_RNDN);
290 
291     return mpfrInput.as<double>();
292   }
293 };
294 
295 namespace internal {
296 
297 template <typename InputType>
298 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
299 unaryOperation(Operation op, InputType input) {
300   MPFRNumber mpfrInput(input);
301   switch (op) {
302   case Operation::Abs:
303     return mpfrInput.abs();
304   case Operation::Ceil:
305     return mpfrInput.ceil();
306   case Operation::Cos:
307     return mpfrInput.cos();
308   case Operation::Exp:
309     return mpfrInput.exp();
310   case Operation::Exp2:
311     return mpfrInput.exp2();
312   case Operation::Floor:
313     return mpfrInput.floor();
314   case Operation::Round:
315     return mpfrInput.round();
316   case Operation::Sin:
317     return mpfrInput.sin();
318   case Operation::Sqrt:
319     return mpfrInput.sqrt();
320   case Operation::Tan:
321     return mpfrInput.tan();
322   case Operation::Trunc:
323     return mpfrInput.trunc();
324   default:
325     __builtin_unreachable();
326   }
327 }
328 
329 template <typename InputType>
330 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
331 unaryOperationTwoOutputs(Operation op, InputType input, int &output) {
332   MPFRNumber mpfrInput(input);
333   switch (op) {
334   case Operation::Frexp:
335     return mpfrInput.frexp(output);
336   default:
337     __builtin_unreachable();
338   }
339 }
340 
341 template <typename InputType>
342 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
343 binaryOperationOneOutput(Operation op, InputType x, InputType y) {
344   MPFRNumber inputX(x), inputY(y);
345   switch (op) {
346   case Operation::Hypot:
347     return inputX.hypot(inputY);
348   default:
349     __builtin_unreachable();
350   }
351 }
352 
353 template <typename InputType>
354 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
355 binaryOperationTwoOutputs(Operation op, InputType x, InputType y, int &output) {
356   MPFRNumber inputX(x), inputY(y);
357   switch (op) {
358   case Operation::RemQuo:
359     return inputX.remquo(inputY, output);
360   default:
361     __builtin_unreachable();
362   }
363 }
364 
365 template <typename InputType>
366 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
367 ternaryOperationOneOutput(Operation op, InputType x, InputType y, InputType z) {
368   // For FMA function, we just need to compare with the mpfr_fma with the same
369   // precision as InputType.  Using higher precision as the intermediate results
370   // to compare might incorrectly fail due to double-rounding errors.
371   constexpr unsigned int prec = Precision<InputType>::value;
372   MPFRNumber inputX(x, prec), inputY(y, prec), inputZ(z, prec);
373   switch (op) {
374   case Operation::Fma:
375     return inputX.fma(inputY, inputZ);
376   default:
377     __builtin_unreachable();
378   }
379 }
380 
381 template <typename T>
382 void explainUnaryOperationSingleOutputError(Operation op, T input, T matchValue,
383                                             testutils::StreamWrapper &OS) {
384   MPFRNumber mpfrInput(input);
385   MPFRNumber mpfrResult = unaryOperation(op, input);
386   MPFRNumber mpfrMatchValue(matchValue);
387   FPBits<T> inputBits(input);
388   FPBits<T> matchBits(matchValue);
389   FPBits<T> mpfrResultBits(mpfrResult.as<T>());
390   OS << "Match value not within tolerance value of MPFR result:\n"
391      << "  Input decimal: " << mpfrInput.str() << '\n';
392   __llvm_libc::fputil::testing::describeValue("     Input bits: ", input, OS);
393   OS << '\n' << "  Match decimal: " << mpfrMatchValue.str() << '\n';
394   __llvm_libc::fputil::testing::describeValue("     Match bits: ", matchValue,
395                                               OS);
396   OS << '\n' << "    MPFR result: " << mpfrResult.str() << '\n';
397   __llvm_libc::fputil::testing::describeValue(
398       "   MPFR rounded: ", mpfrResult.as<T>(), OS);
399   OS << '\n';
400   OS << "      ULP error: " << std::to_string(mpfrResult.ulp(matchValue))
401      << '\n';
402 }
403 
404 template void
405 explainUnaryOperationSingleOutputError<float>(Operation op, float, float,
406                                               testutils::StreamWrapper &);
407 template void
408 explainUnaryOperationSingleOutputError<double>(Operation op, double, double,
409                                                testutils::StreamWrapper &);
410 template void explainUnaryOperationSingleOutputError<long double>(
411     Operation op, long double, long double, testutils::StreamWrapper &);
412 
413 template <typename T>
414 void explainUnaryOperationTwoOutputsError(Operation op, T input,
415                                           const BinaryOutput<T> &libcResult,
416                                           testutils::StreamWrapper &OS) {
417   MPFRNumber mpfrInput(input);
418   FPBits<T> inputBits(input);
419   int mpfrIntResult;
420   MPFRNumber mpfrResult = unaryOperationTwoOutputs(op, input, mpfrIntResult);
421 
422   if (mpfrIntResult != libcResult.i) {
423     OS << "MPFR integral result: " << mpfrIntResult << '\n'
424        << "Libc integral result: " << libcResult.i << '\n';
425   } else {
426     OS << "Integral result from libc matches integral result from MPFR.\n";
427   }
428 
429   MPFRNumber mpfrMatchValue(libcResult.f);
430   OS << "Libc floating point result is not within tolerance value of the MPFR "
431      << "result.\n\n";
432 
433   OS << "            Input decimal: " << mpfrInput.str() << "\n\n";
434 
435   OS << "Libc floating point value: " << mpfrMatchValue.str() << '\n';
436   __llvm_libc::fputil::testing::describeValue(
437       " Libc floating point bits: ", libcResult.f, OS);
438   OS << "\n\n";
439 
440   OS << "              MPFR result: " << mpfrResult.str() << '\n';
441   __llvm_libc::fputil::testing::describeValue(
442       "             MPFR rounded: ", mpfrResult.as<T>(), OS);
443   OS << '\n'
444      << "                ULP error: "
445      << std::to_string(mpfrResult.ulp(libcResult.f)) << '\n';
446 }
447 
448 template void explainUnaryOperationTwoOutputsError<float>(
449     Operation, float, const BinaryOutput<float> &, testutils::StreamWrapper &);
450 template void
451 explainUnaryOperationTwoOutputsError<double>(Operation, double,
452                                              const BinaryOutput<double> &,
453                                              testutils::StreamWrapper &);
454 template void explainUnaryOperationTwoOutputsError<long double>(
455     Operation, long double, const BinaryOutput<long double> &,
456     testutils::StreamWrapper &);
457 
458 template <typename T>
459 void explainBinaryOperationTwoOutputsError(Operation op,
460                                            const BinaryInput<T> &input,
461                                            const BinaryOutput<T> &libcResult,
462                                            testutils::StreamWrapper &OS) {
463   MPFRNumber mpfrX(input.x);
464   MPFRNumber mpfrY(input.y);
465   FPBits<T> xbits(input.x);
466   FPBits<T> ybits(input.y);
467   int mpfrIntResult;
468   MPFRNumber mpfrResult =
469       binaryOperationTwoOutputs(op, input.x, input.y, mpfrIntResult);
470   MPFRNumber mpfrMatchValue(libcResult.f);
471 
472   OS << "Input decimal: x: " << mpfrX.str() << " y: " << mpfrY.str() << '\n'
473      << "MPFR integral result: " << mpfrIntResult << '\n'
474      << "Libc integral result: " << libcResult.i << '\n'
475      << "Libc floating point result: " << mpfrMatchValue.str() << '\n'
476      << "               MPFR result: " << mpfrResult.str() << '\n';
477   __llvm_libc::fputil::testing::describeValue(
478       "Libc floating point result bits: ", libcResult.f, OS);
479   __llvm_libc::fputil::testing::describeValue(
480       "              MPFR rounded bits: ", mpfrResult.as<T>(), OS);
481   OS << "ULP error: " << std::to_string(mpfrResult.ulp(libcResult.f)) << '\n';
482 }
483 
484 template void explainBinaryOperationTwoOutputsError<float>(
485     Operation, const BinaryInput<float> &, const BinaryOutput<float> &,
486     testutils::StreamWrapper &);
487 template void explainBinaryOperationTwoOutputsError<double>(
488     Operation, const BinaryInput<double> &, const BinaryOutput<double> &,
489     testutils::StreamWrapper &);
490 template void explainBinaryOperationTwoOutputsError<long double>(
491     Operation, const BinaryInput<long double> &,
492     const BinaryOutput<long double> &, testutils::StreamWrapper &);
493 
494 template <typename T>
495 void explainBinaryOperationOneOutputError(Operation op,
496                                           const BinaryInput<T> &input,
497                                           T libcResult,
498                                           testutils::StreamWrapper &OS) {
499   MPFRNumber mpfrX(input.x);
500   MPFRNumber mpfrY(input.y);
501   FPBits<T> xbits(input.x);
502   FPBits<T> ybits(input.y);
503   MPFRNumber mpfrResult = binaryOperationOneOutput(op, input.x, input.y);
504   MPFRNumber mpfrMatchValue(libcResult);
505 
506   OS << "Input decimal: x: " << mpfrX.str() << " y: " << mpfrY.str() << '\n';
507   __llvm_libc::fputil::testing::describeValue("First input bits: ", input.x,
508                                               OS);
509   __llvm_libc::fputil::testing::describeValue("Second input bits: ", input.y,
510                                               OS);
511 
512   OS << "Libc result: " << mpfrMatchValue.str() << '\n'
513      << "MPFR result: " << mpfrResult.str() << '\n';
514   __llvm_libc::fputil::testing::describeValue(
515       "Libc floating point result bits: ", libcResult, OS);
516   __llvm_libc::fputil::testing::describeValue(
517       "              MPFR rounded bits: ", mpfrResult.as<T>(), OS);
518   OS << "ULP error: " << std::to_string(mpfrResult.ulp(libcResult)) << '\n';
519 }
520 
521 template void explainBinaryOperationOneOutputError<float>(
522     Operation, const BinaryInput<float> &, float, testutils::StreamWrapper &);
523 template void explainBinaryOperationOneOutputError<double>(
524     Operation, const BinaryInput<double> &, double, testutils::StreamWrapper &);
525 template void explainBinaryOperationOneOutputError<long double>(
526     Operation, const BinaryInput<long double> &, long double,
527     testutils::StreamWrapper &);
528 
529 template <typename T>
530 void explainTernaryOperationOneOutputError(Operation op,
531                                            const TernaryInput<T> &input,
532                                            T libcResult,
533                                            testutils::StreamWrapper &OS) {
534   MPFRNumber mpfrX(input.x, Precision<T>::value);
535   MPFRNumber mpfrY(input.y, Precision<T>::value);
536   MPFRNumber mpfrZ(input.z, Precision<T>::value);
537   FPBits<T> xbits(input.x);
538   FPBits<T> ybits(input.y);
539   FPBits<T> zbits(input.z);
540   MPFRNumber mpfrResult =
541       ternaryOperationOneOutput(op, input.x, input.y, input.z);
542   MPFRNumber mpfrMatchValue(libcResult);
543 
544   OS << "Input decimal: x: " << mpfrX.str() << " y: " << mpfrY.str()
545      << " z: " << mpfrZ.str() << '\n';
546   __llvm_libc::fputil::testing::describeValue("First input bits: ", input.x,
547                                               OS);
548   __llvm_libc::fputil::testing::describeValue("Second input bits: ", input.y,
549                                               OS);
550   __llvm_libc::fputil::testing::describeValue("Third input bits: ", input.z,
551                                               OS);
552 
553   OS << "Libc result: " << mpfrMatchValue.str() << '\n'
554      << "MPFR result: " << mpfrResult.str() << '\n';
555   __llvm_libc::fputil::testing::describeValue(
556       "Libc floating point result bits: ", libcResult, OS);
557   __llvm_libc::fputil::testing::describeValue(
558       "              MPFR rounded bits: ", mpfrResult.as<T>(), OS);
559   OS << "ULP error: " << std::to_string(mpfrResult.ulp(libcResult)) << '\n';
560 }
561 
562 template void explainTernaryOperationOneOutputError<float>(
563     Operation, const TernaryInput<float> &, float, testutils::StreamWrapper &);
564 template void explainTernaryOperationOneOutputError<double>(
565     Operation, const TernaryInput<double> &, double,
566     testutils::StreamWrapper &);
567 template void explainTernaryOperationOneOutputError<long double>(
568     Operation, const TernaryInput<long double> &, long double,
569     testutils::StreamWrapper &);
570 
571 template <typename T>
572 bool compareUnaryOperationSingleOutput(Operation op, T input, T libcResult,
573                                        double ulpError) {
574   // If the ulp error is exactly 0.5 (i.e a tie), we would check that the result
575   // is rounded to the nearest even.
576   MPFRNumber mpfrResult = unaryOperation(op, input);
577   double ulp = mpfrResult.ulp(libcResult);
578   bool bitsAreEven = ((FPBits<T>(libcResult).uintval() & 1) == 0);
579   return (ulp < ulpError) ||
580          ((ulp == ulpError) && ((ulp != 0.5) || bitsAreEven));
581 }
582 
583 template bool compareUnaryOperationSingleOutput<float>(Operation, float, float,
584                                                        double);
585 template bool compareUnaryOperationSingleOutput<double>(Operation, double,
586                                                         double, double);
587 template bool compareUnaryOperationSingleOutput<long double>(Operation,
588                                                              long double,
589                                                              long double,
590                                                              double);
591 
592 template <typename T>
593 bool compareUnaryOperationTwoOutputs(Operation op, T input,
594                                      const BinaryOutput<T> &libcResult,
595                                      double ulpError) {
596   int mpfrIntResult;
597   MPFRNumber mpfrResult = unaryOperationTwoOutputs(op, input, mpfrIntResult);
598   double ulp = mpfrResult.ulp(libcResult.f);
599 
600   if (mpfrIntResult != libcResult.i)
601     return false;
602 
603   bool bitsAreEven = ((FPBits<T>(libcResult.f).uintval() & 1) == 0);
604   return (ulp < ulpError) ||
605          ((ulp == ulpError) && ((ulp != 0.5) || bitsAreEven));
606 }
607 
608 template bool
609 compareUnaryOperationTwoOutputs<float>(Operation, float,
610                                        const BinaryOutput<float> &, double);
611 template bool
612 compareUnaryOperationTwoOutputs<double>(Operation, double,
613                                         const BinaryOutput<double> &, double);
614 template bool compareUnaryOperationTwoOutputs<long double>(
615     Operation, long double, const BinaryOutput<long double> &, double);
616 
617 template <typename T>
618 bool compareBinaryOperationTwoOutputs(Operation op, const BinaryInput<T> &input,
619                                       const BinaryOutput<T> &libcResult,
620                                       double ulpError) {
621   int mpfrIntResult;
622   MPFRNumber mpfrResult =
623       binaryOperationTwoOutputs(op, input.x, input.y, mpfrIntResult);
624   double ulp = mpfrResult.ulp(libcResult.f);
625 
626   if (mpfrIntResult != libcResult.i) {
627     if (op == Operation::RemQuo) {
628       if ((0x7 & mpfrIntResult) != (0x7 & libcResult.i))
629         return false;
630     } else {
631       return false;
632     }
633   }
634 
635   bool bitsAreEven = ((FPBits<T>(libcResult.f).uintval() & 1) == 0);
636   return (ulp < ulpError) ||
637          ((ulp == ulpError) && ((ulp != 0.5) || bitsAreEven));
638 }
639 
640 template bool
641 compareBinaryOperationTwoOutputs<float>(Operation, const BinaryInput<float> &,
642                                         const BinaryOutput<float> &, double);
643 template bool
644 compareBinaryOperationTwoOutputs<double>(Operation, const BinaryInput<double> &,
645                                          const BinaryOutput<double> &, double);
646 template bool compareBinaryOperationTwoOutputs<long double>(
647     Operation, const BinaryInput<long double> &,
648     const BinaryOutput<long double> &, double);
649 
650 template <typename T>
651 bool compareBinaryOperationOneOutput(Operation op, const BinaryInput<T> &input,
652                                      T libcResult, double ulpError) {
653   MPFRNumber mpfrResult = binaryOperationOneOutput(op, input.x, input.y);
654   double ulp = mpfrResult.ulp(libcResult);
655 
656   bool bitsAreEven = ((FPBits<T>(libcResult).uintval() & 1) == 0);
657   return (ulp < ulpError) ||
658          ((ulp == ulpError) && ((ulp != 0.5) || bitsAreEven));
659 }
660 
661 template bool compareBinaryOperationOneOutput<float>(Operation,
662                                                      const BinaryInput<float> &,
663                                                      float, double);
664 template bool
665 compareBinaryOperationOneOutput<double>(Operation, const BinaryInput<double> &,
666                                         double, double);
667 template bool compareBinaryOperationOneOutput<long double>(
668     Operation, const BinaryInput<long double> &, long double, double);
669 
670 template <typename T>
671 bool compareTernaryOperationOneOutput(Operation op,
672                                       const TernaryInput<T> &input,
673                                       T libcResult, double ulpError) {
674   MPFRNumber mpfrResult =
675       ternaryOperationOneOutput(op, input.x, input.y, input.z);
676   double ulp = mpfrResult.ulp(libcResult);
677 
678   bool bitsAreEven = ((FPBits<T>(libcResult).uintval() & 1) == 0);
679   return (ulp < ulpError) ||
680          ((ulp == ulpError) && ((ulp != 0.5) || bitsAreEven));
681 }
682 
683 template bool
684 compareTernaryOperationOneOutput<float>(Operation, const TernaryInput<float> &,
685                                         float, double);
686 template bool compareTernaryOperationOneOutput<double>(
687     Operation, const TernaryInput<double> &, double, double);
688 template bool compareTernaryOperationOneOutput<long double>(
689     Operation, const TernaryInput<long double> &, long double, double);
690 
691 static mpfr_rnd_t getMPFRRoundingMode(RoundingMode mode) {
692   switch (mode) {
693   case RoundingMode::Upward:
694     return MPFR_RNDU;
695     break;
696   case RoundingMode::Downward:
697     return MPFR_RNDD;
698     break;
699   case RoundingMode::TowardZero:
700     return MPFR_RNDZ;
701     break;
702   case RoundingMode::Nearest:
703     return MPFR_RNDN;
704     break;
705   }
706 }
707 
708 } // namespace internal
709 
710 template <typename T> bool RoundToLong(T x, long &result) {
711   MPFRNumber mpfr(x);
712   return mpfr.roundToLong(result);
713 }
714 
715 template bool RoundToLong<float>(float, long &);
716 template bool RoundToLong<double>(double, long &);
717 template bool RoundToLong<long double>(long double, long &);
718 
719 template <typename T> bool RoundToLong(T x, RoundingMode mode, long &result) {
720   MPFRNumber mpfr(x);
721   return mpfr.roundToLong(internal::getMPFRRoundingMode(mode), result);
722 }
723 
724 template bool RoundToLong<float>(float, RoundingMode, long &);
725 template bool RoundToLong<double>(double, RoundingMode, long &);
726 template bool RoundToLong<long double>(long double, RoundingMode, long &);
727 
728 template <typename T> T Round(T x, RoundingMode mode) {
729   MPFRNumber mpfr(x);
730   MPFRNumber result = mpfr.rint(internal::getMPFRRoundingMode(mode));
731   return result.as<T>();
732 }
733 
734 template float Round<float>(float, RoundingMode);
735 template double Round<double>(double, RoundingMode);
736 template long double Round<long double>(long double, RoundingMode);
737 
738 } // namespace mpfr
739 } // namespace testing
740 } // namespace __llvm_libc
741