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