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     // When calculating error, we first round this number to the floating
274     // point type T and calculate the ulp error against this rounded number.
275     // The input is always either exact or rounded. So, we if compare
276     // with this number directly, we can end up with a large ULP error.
277     MPFRNumber thisRoundedToFloat(as<T>());
278 
279     // abs(thisRoundedToFloat - input)
280     mpfr_sub(mpfrInput.value, thisRoundedToFloat.value, mpfrInput.value,
281              MPFR_RNDN);
282     mpfr_abs(mpfrInput.value, mpfrInput.value, MPFR_RNDN);
283 
284     // get eps(input)
285     int epsExponent = bits.encoding.exponent - fputil::FPBits<T>::exponentBias -
286                       fputil::MantissaWidth<T>::value;
287     if (bits.encoding.exponent == 0) {
288       // correcting denormal exponent
289       ++epsExponent;
290     } else if ((bits.encoding.mantissa == 0) && (bits.encoding.exponent > 1) &&
291                mpfr_less_p(thisRoundedToFloat.value, mpfrInput.value)) {
292       // when the input is exactly 2^n, distance (epsilon) between the input
293       // and the next floating point number is different from the distance to
294       // the previous floating point number.  So in that case, if the correct
295       // value from MPFR is smaller than the input, we use the smaller epsilon
296       --epsExponent;
297     }
298 
299     // Since eps(value) is of the form 2^e, instead of dividing such number,
300     // we multiply by its inverse 2^{-e}.
301     mpfr_mul_2si(mpfrInput.value, mpfrInput.value, -epsExponent, MPFR_RNDN);
302 
303     return mpfrInput.as<double>();
304   }
305 };
306 
307 namespace internal {
308 
309 template <typename InputType>
310 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
311 unaryOperation(Operation op, InputType input) {
312   MPFRNumber mpfrInput(input);
313   switch (op) {
314   case Operation::Abs:
315     return mpfrInput.abs();
316   case Operation::Ceil:
317     return mpfrInput.ceil();
318   case Operation::Cos:
319     return mpfrInput.cos();
320   case Operation::Exp:
321     return mpfrInput.exp();
322   case Operation::Exp2:
323     return mpfrInput.exp2();
324   case Operation::Expm1:
325     return mpfrInput.expm1();
326   case Operation::Floor:
327     return mpfrInput.floor();
328   case Operation::Round:
329     return mpfrInput.round();
330   case Operation::Sin:
331     return mpfrInput.sin();
332   case Operation::Sqrt:
333     return mpfrInput.sqrt();
334   case Operation::Tan:
335     return mpfrInput.tan();
336   case Operation::Trunc:
337     return mpfrInput.trunc();
338   default:
339     __builtin_unreachable();
340   }
341 }
342 
343 template <typename InputType>
344 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
345 unaryOperationTwoOutputs(Operation op, InputType input, int &output) {
346   MPFRNumber mpfrInput(input);
347   switch (op) {
348   case Operation::Frexp:
349     return mpfrInput.frexp(output);
350   default:
351     __builtin_unreachable();
352   }
353 }
354 
355 template <typename InputType>
356 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
357 binaryOperationOneOutput(Operation op, InputType x, InputType y) {
358   MPFRNumber inputX(x), inputY(y);
359   switch (op) {
360   case Operation::Hypot:
361     return inputX.hypot(inputY);
362   default:
363     __builtin_unreachable();
364   }
365 }
366 
367 template <typename InputType>
368 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
369 binaryOperationTwoOutputs(Operation op, InputType x, InputType y, int &output) {
370   MPFRNumber inputX(x), inputY(y);
371   switch (op) {
372   case Operation::RemQuo:
373     return inputX.remquo(inputY, output);
374   default:
375     __builtin_unreachable();
376   }
377 }
378 
379 template <typename InputType>
380 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
381 ternaryOperationOneOutput(Operation op, InputType x, InputType y, InputType z) {
382   // For FMA function, we just need to compare with the mpfr_fma with the same
383   // precision as InputType.  Using higher precision as the intermediate results
384   // to compare might incorrectly fail due to double-rounding errors.
385   constexpr unsigned int prec = Precision<InputType>::value;
386   MPFRNumber inputX(x, prec), inputY(y, prec), inputZ(z, prec);
387   switch (op) {
388   case Operation::Fma:
389     return inputX.fma(inputY, inputZ);
390   default:
391     __builtin_unreachable();
392   }
393 }
394 
395 template <typename T>
396 void explainUnaryOperationSingleOutputError(Operation op, T input, T matchValue,
397                                             testutils::StreamWrapper &OS) {
398   MPFRNumber mpfrInput(input);
399   MPFRNumber mpfrResult = unaryOperation(op, input);
400   MPFRNumber mpfrMatchValue(matchValue);
401   FPBits<T> inputBits(input);
402   FPBits<T> matchBits(matchValue);
403   FPBits<T> mpfrResultBits(mpfrResult.as<T>());
404   OS << "Match value not within tolerance value of MPFR result:\n"
405      << "  Input decimal: " << mpfrInput.str() << '\n';
406   __llvm_libc::fputil::testing::describeValue("     Input bits: ", input, OS);
407   OS << '\n' << "  Match decimal: " << mpfrMatchValue.str() << '\n';
408   __llvm_libc::fputil::testing::describeValue("     Match bits: ", matchValue,
409                                               OS);
410   OS << '\n' << "    MPFR result: " << mpfrResult.str() << '\n';
411   __llvm_libc::fputil::testing::describeValue(
412       "   MPFR rounded: ", mpfrResult.as<T>(), OS);
413   OS << '\n';
414   OS << "      ULP error: " << std::to_string(mpfrResult.ulp(matchValue))
415      << '\n';
416 }
417 
418 template void
419 explainUnaryOperationSingleOutputError<float>(Operation op, float, float,
420                                               testutils::StreamWrapper &);
421 template void
422 explainUnaryOperationSingleOutputError<double>(Operation op, double, double,
423                                                testutils::StreamWrapper &);
424 template void explainUnaryOperationSingleOutputError<long double>(
425     Operation op, long double, long double, testutils::StreamWrapper &);
426 
427 template <typename T>
428 void explainUnaryOperationTwoOutputsError(Operation op, T input,
429                                           const BinaryOutput<T> &libcResult,
430                                           testutils::StreamWrapper &OS) {
431   MPFRNumber mpfrInput(input);
432   FPBits<T> inputBits(input);
433   int mpfrIntResult;
434   MPFRNumber mpfrResult = unaryOperationTwoOutputs(op, input, mpfrIntResult);
435 
436   if (mpfrIntResult != libcResult.i) {
437     OS << "MPFR integral result: " << mpfrIntResult << '\n'
438        << "Libc integral result: " << libcResult.i << '\n';
439   } else {
440     OS << "Integral result from libc matches integral result from MPFR.\n";
441   }
442 
443   MPFRNumber mpfrMatchValue(libcResult.f);
444   OS << "Libc floating point result is not within tolerance value of the MPFR "
445      << "result.\n\n";
446 
447   OS << "            Input decimal: " << mpfrInput.str() << "\n\n";
448 
449   OS << "Libc floating point value: " << mpfrMatchValue.str() << '\n';
450   __llvm_libc::fputil::testing::describeValue(
451       " Libc floating point bits: ", libcResult.f, OS);
452   OS << "\n\n";
453 
454   OS << "              MPFR result: " << mpfrResult.str() << '\n';
455   __llvm_libc::fputil::testing::describeValue(
456       "             MPFR rounded: ", mpfrResult.as<T>(), OS);
457   OS << '\n'
458      << "                ULP error: "
459      << std::to_string(mpfrResult.ulp(libcResult.f)) << '\n';
460 }
461 
462 template void explainUnaryOperationTwoOutputsError<float>(
463     Operation, float, const BinaryOutput<float> &, testutils::StreamWrapper &);
464 template void
465 explainUnaryOperationTwoOutputsError<double>(Operation, double,
466                                              const BinaryOutput<double> &,
467                                              testutils::StreamWrapper &);
468 template void explainUnaryOperationTwoOutputsError<long double>(
469     Operation, long double, const BinaryOutput<long double> &,
470     testutils::StreamWrapper &);
471 
472 template <typename T>
473 void explainBinaryOperationTwoOutputsError(Operation op,
474                                            const BinaryInput<T> &input,
475                                            const BinaryOutput<T> &libcResult,
476                                            testutils::StreamWrapper &OS) {
477   MPFRNumber mpfrX(input.x);
478   MPFRNumber mpfrY(input.y);
479   FPBits<T> xbits(input.x);
480   FPBits<T> ybits(input.y);
481   int mpfrIntResult;
482   MPFRNumber mpfrResult =
483       binaryOperationTwoOutputs(op, input.x, input.y, mpfrIntResult);
484   MPFRNumber mpfrMatchValue(libcResult.f);
485 
486   OS << "Input decimal: x: " << mpfrX.str() << " y: " << mpfrY.str() << '\n'
487      << "MPFR integral result: " << mpfrIntResult << '\n'
488      << "Libc integral result: " << libcResult.i << '\n'
489      << "Libc floating point result: " << mpfrMatchValue.str() << '\n'
490      << "               MPFR result: " << mpfrResult.str() << '\n';
491   __llvm_libc::fputil::testing::describeValue(
492       "Libc floating point result bits: ", libcResult.f, OS);
493   __llvm_libc::fputil::testing::describeValue(
494       "              MPFR rounded bits: ", mpfrResult.as<T>(), OS);
495   OS << "ULP error: " << std::to_string(mpfrResult.ulp(libcResult.f)) << '\n';
496 }
497 
498 template void explainBinaryOperationTwoOutputsError<float>(
499     Operation, const BinaryInput<float> &, const BinaryOutput<float> &,
500     testutils::StreamWrapper &);
501 template void explainBinaryOperationTwoOutputsError<double>(
502     Operation, const BinaryInput<double> &, const BinaryOutput<double> &,
503     testutils::StreamWrapper &);
504 template void explainBinaryOperationTwoOutputsError<long double>(
505     Operation, const BinaryInput<long double> &,
506     const BinaryOutput<long double> &, testutils::StreamWrapper &);
507 
508 template <typename T>
509 void explainBinaryOperationOneOutputError(Operation op,
510                                           const BinaryInput<T> &input,
511                                           T libcResult,
512                                           testutils::StreamWrapper &OS) {
513   MPFRNumber mpfrX(input.x);
514   MPFRNumber mpfrY(input.y);
515   FPBits<T> xbits(input.x);
516   FPBits<T> ybits(input.y);
517   MPFRNumber mpfrResult = binaryOperationOneOutput(op, input.x, input.y);
518   MPFRNumber mpfrMatchValue(libcResult);
519 
520   OS << "Input decimal: x: " << mpfrX.str() << " y: " << mpfrY.str() << '\n';
521   __llvm_libc::fputil::testing::describeValue("First input bits: ", input.x,
522                                               OS);
523   __llvm_libc::fputil::testing::describeValue("Second input bits: ", input.y,
524                                               OS);
525 
526   OS << "Libc result: " << mpfrMatchValue.str() << '\n'
527      << "MPFR result: " << mpfrResult.str() << '\n';
528   __llvm_libc::fputil::testing::describeValue(
529       "Libc floating point result bits: ", libcResult, OS);
530   __llvm_libc::fputil::testing::describeValue(
531       "              MPFR rounded bits: ", mpfrResult.as<T>(), OS);
532   OS << "ULP error: " << std::to_string(mpfrResult.ulp(libcResult)) << '\n';
533 }
534 
535 template void explainBinaryOperationOneOutputError<float>(
536     Operation, const BinaryInput<float> &, float, testutils::StreamWrapper &);
537 template void explainBinaryOperationOneOutputError<double>(
538     Operation, const BinaryInput<double> &, double, testutils::StreamWrapper &);
539 template void explainBinaryOperationOneOutputError<long double>(
540     Operation, const BinaryInput<long double> &, long double,
541     testutils::StreamWrapper &);
542 
543 template <typename T>
544 void explainTernaryOperationOneOutputError(Operation op,
545                                            const TernaryInput<T> &input,
546                                            T libcResult,
547                                            testutils::StreamWrapper &OS) {
548   MPFRNumber mpfrX(input.x, Precision<T>::value);
549   MPFRNumber mpfrY(input.y, Precision<T>::value);
550   MPFRNumber mpfrZ(input.z, Precision<T>::value);
551   FPBits<T> xbits(input.x);
552   FPBits<T> ybits(input.y);
553   FPBits<T> zbits(input.z);
554   MPFRNumber mpfrResult =
555       ternaryOperationOneOutput(op, input.x, input.y, input.z);
556   MPFRNumber mpfrMatchValue(libcResult);
557 
558   OS << "Input decimal: x: " << mpfrX.str() << " y: " << mpfrY.str()
559      << " z: " << mpfrZ.str() << '\n';
560   __llvm_libc::fputil::testing::describeValue("First input bits: ", input.x,
561                                               OS);
562   __llvm_libc::fputil::testing::describeValue("Second input bits: ", input.y,
563                                               OS);
564   __llvm_libc::fputil::testing::describeValue("Third input bits: ", input.z,
565                                               OS);
566 
567   OS << "Libc result: " << mpfrMatchValue.str() << '\n'
568      << "MPFR result: " << mpfrResult.str() << '\n';
569   __llvm_libc::fputil::testing::describeValue(
570       "Libc floating point result bits: ", libcResult, OS);
571   __llvm_libc::fputil::testing::describeValue(
572       "              MPFR rounded bits: ", mpfrResult.as<T>(), OS);
573   OS << "ULP error: " << std::to_string(mpfrResult.ulp(libcResult)) << '\n';
574 }
575 
576 template void explainTernaryOperationOneOutputError<float>(
577     Operation, const TernaryInput<float> &, float, testutils::StreamWrapper &);
578 template void explainTernaryOperationOneOutputError<double>(
579     Operation, const TernaryInput<double> &, double,
580     testutils::StreamWrapper &);
581 template void explainTernaryOperationOneOutputError<long double>(
582     Operation, const TernaryInput<long double> &, long double,
583     testutils::StreamWrapper &);
584 
585 template <typename T>
586 bool compareUnaryOperationSingleOutput(Operation op, T input, T libcResult,
587                                        double ulpError) {
588   // If the ulp error is exactly 0.5 (i.e a tie), we would check that the result
589   // is rounded to the nearest even.
590   MPFRNumber mpfrResult = unaryOperation(op, input);
591   double ulp = mpfrResult.ulp(libcResult);
592   bool bitsAreEven = ((FPBits<T>(libcResult).uintval() & 1) == 0);
593   return (ulp < ulpError) ||
594          ((ulp == ulpError) && ((ulp != 0.5) || bitsAreEven));
595 }
596 
597 template bool compareUnaryOperationSingleOutput<float>(Operation, float, float,
598                                                        double);
599 template bool compareUnaryOperationSingleOutput<double>(Operation, double,
600                                                         double, double);
601 template bool compareUnaryOperationSingleOutput<long double>(Operation,
602                                                              long double,
603                                                              long double,
604                                                              double);
605 
606 template <typename T>
607 bool compareUnaryOperationTwoOutputs(Operation op, T input,
608                                      const BinaryOutput<T> &libcResult,
609                                      double ulpError) {
610   int mpfrIntResult;
611   MPFRNumber mpfrResult = unaryOperationTwoOutputs(op, input, mpfrIntResult);
612   double ulp = mpfrResult.ulp(libcResult.f);
613 
614   if (mpfrIntResult != libcResult.i)
615     return false;
616 
617   bool bitsAreEven = ((FPBits<T>(libcResult.f).uintval() & 1) == 0);
618   return (ulp < ulpError) ||
619          ((ulp == ulpError) && ((ulp != 0.5) || bitsAreEven));
620 }
621 
622 template bool
623 compareUnaryOperationTwoOutputs<float>(Operation, float,
624                                        const BinaryOutput<float> &, double);
625 template bool
626 compareUnaryOperationTwoOutputs<double>(Operation, double,
627                                         const BinaryOutput<double> &, double);
628 template bool compareUnaryOperationTwoOutputs<long double>(
629     Operation, long double, const BinaryOutput<long double> &, double);
630 
631 template <typename T>
632 bool compareBinaryOperationTwoOutputs(Operation op, const BinaryInput<T> &input,
633                                       const BinaryOutput<T> &libcResult,
634                                       double ulpError) {
635   int mpfrIntResult;
636   MPFRNumber mpfrResult =
637       binaryOperationTwoOutputs(op, input.x, input.y, mpfrIntResult);
638   double ulp = mpfrResult.ulp(libcResult.f);
639 
640   if (mpfrIntResult != libcResult.i) {
641     if (op == Operation::RemQuo) {
642       if ((0x7 & mpfrIntResult) != (0x7 & libcResult.i))
643         return false;
644     } else {
645       return false;
646     }
647   }
648 
649   bool bitsAreEven = ((FPBits<T>(libcResult.f).uintval() & 1) == 0);
650   return (ulp < ulpError) ||
651          ((ulp == ulpError) && ((ulp != 0.5) || bitsAreEven));
652 }
653 
654 template bool
655 compareBinaryOperationTwoOutputs<float>(Operation, const BinaryInput<float> &,
656                                         const BinaryOutput<float> &, double);
657 template bool
658 compareBinaryOperationTwoOutputs<double>(Operation, const BinaryInput<double> &,
659                                          const BinaryOutput<double> &, double);
660 template bool compareBinaryOperationTwoOutputs<long double>(
661     Operation, const BinaryInput<long double> &,
662     const BinaryOutput<long double> &, double);
663 
664 template <typename T>
665 bool compareBinaryOperationOneOutput(Operation op, const BinaryInput<T> &input,
666                                      T libcResult, double ulpError) {
667   MPFRNumber mpfrResult = binaryOperationOneOutput(op, input.x, input.y);
668   double ulp = mpfrResult.ulp(libcResult);
669 
670   bool bitsAreEven = ((FPBits<T>(libcResult).uintval() & 1) == 0);
671   return (ulp < ulpError) ||
672          ((ulp == ulpError) && ((ulp != 0.5) || bitsAreEven));
673 }
674 
675 template bool compareBinaryOperationOneOutput<float>(Operation,
676                                                      const BinaryInput<float> &,
677                                                      float, double);
678 template bool
679 compareBinaryOperationOneOutput<double>(Operation, const BinaryInput<double> &,
680                                         double, double);
681 template bool compareBinaryOperationOneOutput<long double>(
682     Operation, const BinaryInput<long double> &, long double, double);
683 
684 template <typename T>
685 bool compareTernaryOperationOneOutput(Operation op,
686                                       const TernaryInput<T> &input,
687                                       T libcResult, double ulpError) {
688   MPFRNumber mpfrResult =
689       ternaryOperationOneOutput(op, input.x, input.y, input.z);
690   double ulp = mpfrResult.ulp(libcResult);
691 
692   bool bitsAreEven = ((FPBits<T>(libcResult).uintval() & 1) == 0);
693   return (ulp < ulpError) ||
694          ((ulp == ulpError) && ((ulp != 0.5) || bitsAreEven));
695 }
696 
697 template bool
698 compareTernaryOperationOneOutput<float>(Operation, const TernaryInput<float> &,
699                                         float, double);
700 template bool compareTernaryOperationOneOutput<double>(
701     Operation, const TernaryInput<double> &, double, double);
702 template bool compareTernaryOperationOneOutput<long double>(
703     Operation, const TernaryInput<long double> &, long double, double);
704 
705 static mpfr_rnd_t getMPFRRoundingMode(RoundingMode mode) {
706   switch (mode) {
707   case RoundingMode::Upward:
708     return MPFR_RNDU;
709     break;
710   case RoundingMode::Downward:
711     return MPFR_RNDD;
712     break;
713   case RoundingMode::TowardZero:
714     return MPFR_RNDZ;
715     break;
716   case RoundingMode::Nearest:
717     return MPFR_RNDN;
718     break;
719   }
720 }
721 
722 } // namespace internal
723 
724 template <typename T> bool RoundToLong(T x, long &result) {
725   MPFRNumber mpfr(x);
726   return mpfr.roundToLong(result);
727 }
728 
729 template bool RoundToLong<float>(float, long &);
730 template bool RoundToLong<double>(double, long &);
731 template bool RoundToLong<long double>(long double, long &);
732 
733 template <typename T> bool RoundToLong(T x, RoundingMode mode, long &result) {
734   MPFRNumber mpfr(x);
735   return mpfr.roundToLong(internal::getMPFRRoundingMode(mode), result);
736 }
737 
738 template bool RoundToLong<float>(float, RoundingMode, long &);
739 template bool RoundToLong<double>(double, RoundingMode, long &);
740 template bool RoundToLong<long double>(long double, RoundingMode, long &);
741 
742 template <typename T> T Round(T x, RoundingMode mode) {
743   MPFRNumber mpfr(x);
744   MPFRNumber result = mpfr.rint(internal::getMPFRRoundingMode(mode));
745   return result.as<T>();
746 }
747 
748 template float Round<float>(float, RoundingMode);
749 template double Round<double>(double, RoundingMode);
750 template long double Round<long double>(long double, RoundingMode);
751 
752 } // namespace mpfr
753 } // namespace testing
754 } // namespace __llvm_libc
755