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