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