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