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 "src/__support/CPP/StringView.h"
12 #include "src/__support/FPUtil/FPBits.h"
13 #include "src/__support/FPUtil/PlatformDefs.h"
14 #include "utils/UnitTest/FPMatcher.h"
15 
16 #include <cmath>
17 #include <fenv.h>
18 #include <memory>
19 #include <sstream>
20 #include <stdint.h>
21 #include <string>
22 
23 #ifdef CUSTOM_MPFR_INCLUDER
24 // Some downstream repos are monoliths carrying MPFR sources in their third
25 // party directory. In such repos, including the MPFR header as
26 // `#include <mpfr.h>` is either disallowed or not possible. If that is the
27 // case, a file named `CustomMPFRIncluder.h` should be added through which the
28 // MPFR header can be included in manner allowed in that repo.
29 #include "CustomMPFRIncluder.h"
30 #else
31 #include <mpfr.h>
32 #endif
33 
34 template <typename T> using FPBits = __llvm_libc::fputil::FPBits<T>;
35 
36 namespace __llvm_libc {
37 namespace testing {
38 namespace mpfr {
39 
40 template <typename T> struct Precision;
41 
42 template <> struct Precision<float> {
43   static constexpr unsigned int VALUE = 24;
44 };
45 
46 template <> struct Precision<double> {
47   static constexpr unsigned int VALUE = 53;
48 };
49 
50 #if defined(LONG_DOUBLE_IS_DOUBLE)
51 template <> struct Precision<long double> {
52   static constexpr unsigned int VALUE = 53;
53 };
54 #elif defined(SPECIAL_X86_LONG_DOUBLE)
55 template <> struct Precision<long double> {
56   static constexpr unsigned int VALUE = 64;
57 };
58 #else
59 template <> struct Precision<long double> {
60   static constexpr unsigned int VALUE = 113;
61 };
62 #endif
63 
64 // A precision value which allows sufficiently large additional
65 // precision compared to the floating point precision.
66 template <typename T> struct ExtraPrecision;
67 
68 template <> struct ExtraPrecision<float> {
69   static constexpr unsigned int VALUE = 128;
70 };
71 
72 template <> struct ExtraPrecision<double> {
73   static constexpr unsigned int VALUE = 256;
74 };
75 
76 template <> struct ExtraPrecision<long double> {
77   static constexpr unsigned int VALUE = 256;
78 };
79 
80 // If the ulp tolerance is less than or equal to 0.5, we would check that the
81 // result is rounded correctly with respect to the rounding mode by using the
82 // same precision as the inputs.
83 template <typename T>
84 static inline unsigned int get_precision(double ulp_tolerance) {
85   if (ulp_tolerance <= 0.5) {
86     return Precision<T>::VALUE;
87   } else {
88     return ExtraPrecision<T>::VALUE;
89   }
90 }
91 
92 static inline mpfr_rnd_t get_mpfr_rounding_mode(RoundingMode mode) {
93   switch (mode) {
94   case RoundingMode::Upward:
95     return MPFR_RNDU;
96     break;
97   case RoundingMode::Downward:
98     return MPFR_RNDD;
99     break;
100   case RoundingMode::TowardZero:
101     return MPFR_RNDZ;
102     break;
103   case RoundingMode::Nearest:
104     return MPFR_RNDN;
105     break;
106   }
107 }
108 
109 int get_fe_rounding(RoundingMode mode) {
110   switch (mode) {
111   case RoundingMode::Upward:
112     return FE_UPWARD;
113     break;
114   case RoundingMode::Downward:
115     return FE_DOWNWARD;
116     break;
117   case RoundingMode::TowardZero:
118     return FE_TOWARDZERO;
119     break;
120   case RoundingMode::Nearest:
121     return FE_TONEAREST;
122     break;
123   }
124 }
125 
126 ForceRoundingMode::ForceRoundingMode(RoundingMode mode) {
127   old_rounding_mode = fegetround();
128   rounding_mode = get_fe_rounding(mode);
129   if (old_rounding_mode != rounding_mode)
130     fesetround(rounding_mode);
131 }
132 
133 ForceRoundingMode::~ForceRoundingMode() {
134   if (old_rounding_mode != rounding_mode)
135     fesetround(old_rounding_mode);
136 }
137 
138 class MPFRNumber {
139   unsigned int mpfr_precision;
140   mpfr_rnd_t mpfr_rounding;
141 
142   mpfr_t value;
143 
144 public:
145   MPFRNumber() : mpfr_precision(256), mpfr_rounding(MPFR_RNDN) {
146     mpfr_init2(value, mpfr_precision);
147   }
148 
149   // We use explicit EnableIf specializations to disallow implicit
150   // conversions. Implicit conversions can potentially lead to loss of
151   // precision.
152   template <typename XType,
153             cpp::EnableIfType<cpp::IsSame<float, XType>::Value, int> = 0>
154   explicit MPFRNumber(XType x, int precision = ExtraPrecision<XType>::VALUE,
155                       RoundingMode rounding = RoundingMode::Nearest)
156       : mpfr_precision(precision),
157         mpfr_rounding(get_mpfr_rounding_mode(rounding)) {
158     mpfr_init2(value, mpfr_precision);
159     mpfr_set_flt(value, x, mpfr_rounding);
160   }
161 
162   template <typename XType,
163             cpp::EnableIfType<cpp::IsSame<double, XType>::Value, int> = 0>
164   explicit MPFRNumber(XType x, int precision = ExtraPrecision<XType>::VALUE,
165                       RoundingMode rounding = RoundingMode::Nearest)
166       : mpfr_precision(precision),
167         mpfr_rounding(get_mpfr_rounding_mode(rounding)) {
168     mpfr_init2(value, mpfr_precision);
169     mpfr_set_d(value, x, mpfr_rounding);
170   }
171 
172   template <typename XType,
173             cpp::EnableIfType<cpp::IsSame<long double, XType>::Value, int> = 0>
174   explicit MPFRNumber(XType x, int precision = ExtraPrecision<XType>::VALUE,
175                       RoundingMode rounding = RoundingMode::Nearest)
176       : mpfr_precision(precision),
177         mpfr_rounding(get_mpfr_rounding_mode(rounding)) {
178     mpfr_init2(value, mpfr_precision);
179     mpfr_set_ld(value, x, mpfr_rounding);
180   }
181 
182   template <typename XType,
183             cpp::EnableIfType<cpp::IsIntegral<XType>::Value, int> = 0>
184   explicit MPFRNumber(XType x, int precision = ExtraPrecision<float>::VALUE,
185                       RoundingMode rounding = RoundingMode::Nearest)
186       : mpfr_precision(precision),
187         mpfr_rounding(get_mpfr_rounding_mode(rounding)) {
188     mpfr_init2(value, mpfr_precision);
189     mpfr_set_sj(value, x, mpfr_rounding);
190   }
191 
192   MPFRNumber(const MPFRNumber &other)
193       : mpfr_precision(other.mpfr_precision),
194         mpfr_rounding(other.mpfr_rounding) {
195     mpfr_init2(value, mpfr_precision);
196     mpfr_set(value, other.value, mpfr_rounding);
197   }
198 
199   ~MPFRNumber() { mpfr_clear(value); }
200 
201   MPFRNumber &operator=(const MPFRNumber &rhs) {
202     mpfr_precision = rhs.mpfr_precision;
203     mpfr_rounding = rhs.mpfr_rounding;
204     mpfr_set(value, rhs.value, mpfr_rounding);
205     return *this;
206   }
207 
208   MPFRNumber abs() const {
209     MPFRNumber result(*this);
210     mpfr_abs(result.value, value, mpfr_rounding);
211     return result;
212   }
213 
214   MPFRNumber ceil() const {
215     MPFRNumber result(*this);
216     mpfr_ceil(result.value, value);
217     return result;
218   }
219 
220   MPFRNumber cos() const {
221     MPFRNumber result(*this);
222     mpfr_cos(result.value, value, mpfr_rounding);
223     return result;
224   }
225 
226   MPFRNumber exp() const {
227     MPFRNumber result(*this);
228     mpfr_exp(result.value, value, mpfr_rounding);
229     return result;
230   }
231 
232   MPFRNumber exp2() const {
233     MPFRNumber result(*this);
234     mpfr_exp2(result.value, value, mpfr_rounding);
235     return result;
236   }
237 
238   MPFRNumber expm1() const {
239     MPFRNumber result(*this);
240     mpfr_expm1(result.value, value, mpfr_rounding);
241     return result;
242   }
243 
244   MPFRNumber floor() const {
245     MPFRNumber result(*this);
246     mpfr_floor(result.value, value);
247     return result;
248   }
249 
250   MPFRNumber frexp(int &exp) {
251     MPFRNumber result(*this);
252     mpfr_exp_t resultExp;
253     mpfr_frexp(&resultExp, result.value, value, mpfr_rounding);
254     exp = resultExp;
255     return result;
256   }
257 
258   MPFRNumber hypot(const MPFRNumber &b) {
259     MPFRNumber result(*this);
260     mpfr_hypot(result.value, value, b.value, mpfr_rounding);
261     return result;
262   }
263 
264   MPFRNumber log() const {
265     MPFRNumber result(*this);
266     mpfr_log(result.value, value, mpfr_rounding);
267     return result;
268   }
269 
270   MPFRNumber log2() const {
271     MPFRNumber result(*this);
272     mpfr_log2(result.value, value, mpfr_rounding);
273     return result;
274   }
275 
276   MPFRNumber log10() const {
277     MPFRNumber result(*this);
278     mpfr_log10(result.value, value, mpfr_rounding);
279     return result;
280   }
281 
282   MPFRNumber remquo(const MPFRNumber &divisor, int &quotient) {
283     MPFRNumber remainder(*this);
284     long q;
285     mpfr_remquo(remainder.value, &q, value, divisor.value, mpfr_rounding);
286     quotient = q;
287     return remainder;
288   }
289 
290   MPFRNumber round() const {
291     MPFRNumber result(*this);
292     mpfr_round(result.value, value);
293     return result;
294   }
295 
296   bool round_to_long(long &result) const {
297     // We first calculate the rounded value. This way, when converting
298     // to long using mpfr_get_si, the rounding direction of MPFR_RNDN
299     // (or any other rounding mode), does not have an influence.
300     MPFRNumber roundedValue = round();
301     mpfr_clear_erangeflag();
302     result = mpfr_get_si(roundedValue.value, MPFR_RNDN);
303     return mpfr_erangeflag_p();
304   }
305 
306   bool round_to_long(mpfr_rnd_t rnd, long &result) const {
307     MPFRNumber rint_result(*this);
308     mpfr_rint(rint_result.value, value, rnd);
309     return rint_result.round_to_long(result);
310   }
311 
312   MPFRNumber rint(mpfr_rnd_t rnd) const {
313     MPFRNumber result(*this);
314     mpfr_rint(result.value, value, rnd);
315     return result;
316   }
317 
318   MPFRNumber mod_2pi() const {
319     MPFRNumber result(0.0, 1280);
320     MPFRNumber _2pi(0.0, 1280);
321     mpfr_const_pi(_2pi.value, MPFR_RNDN);
322     mpfr_mul_si(_2pi.value, _2pi.value, 2, MPFR_RNDN);
323     mpfr_fmod(result.value, value, _2pi.value, MPFR_RNDN);
324     return result;
325   }
326 
327   MPFRNumber mod_pi_over_2() const {
328     MPFRNumber result(0.0, 1280);
329     MPFRNumber pi_over_2(0.0, 1280);
330     mpfr_const_pi(pi_over_2.value, MPFR_RNDN);
331     mpfr_mul_d(pi_over_2.value, pi_over_2.value, 0.5, MPFR_RNDN);
332     mpfr_fmod(result.value, value, pi_over_2.value, MPFR_RNDN);
333     return result;
334   }
335 
336   MPFRNumber mod_pi_over_4() const {
337     MPFRNumber result(0.0, 1280);
338     MPFRNumber pi_over_4(0.0, 1280);
339     mpfr_const_pi(pi_over_4.value, MPFR_RNDN);
340     mpfr_mul_d(pi_over_4.value, pi_over_4.value, 0.25, MPFR_RNDN);
341     mpfr_fmod(result.value, value, pi_over_4.value, MPFR_RNDN);
342     return result;
343   }
344 
345   MPFRNumber sin() const {
346     MPFRNumber result(*this);
347     mpfr_sin(result.value, value, mpfr_rounding);
348     return result;
349   }
350 
351   MPFRNumber sqrt() const {
352     MPFRNumber result(*this);
353     mpfr_sqrt(result.value, value, mpfr_rounding);
354     return result;
355   }
356 
357   MPFRNumber tan() const {
358     MPFRNumber result(*this);
359     mpfr_tan(result.value, value, mpfr_rounding);
360     return result;
361   }
362 
363   MPFRNumber trunc() const {
364     MPFRNumber result(*this);
365     mpfr_trunc(result.value, value);
366     return result;
367   }
368 
369   MPFRNumber fma(const MPFRNumber &b, const MPFRNumber &c) {
370     MPFRNumber result(*this);
371     mpfr_fma(result.value, value, b.value, c.value, mpfr_rounding);
372     return result;
373   }
374 
375   std::string str() const {
376     // 200 bytes should be more than sufficient to hold a 100-digit number
377     // plus additional bytes for the decimal point, '-' sign etc.
378     constexpr size_t printBufSize = 200;
379     char buffer[printBufSize];
380     mpfr_snprintf(buffer, printBufSize, "%100.50Rf", value);
381     cpp::StringView view(buffer);
382     view = view.trim(' ');
383     return std::string(view.data());
384   }
385 
386   // These functions are useful for debugging.
387   template <typename T> T as() const;
388 
389   template <> float as<float>() const {
390     return mpfr_get_flt(value, mpfr_rounding);
391   }
392   template <> double as<double>() const {
393     return mpfr_get_d(value, mpfr_rounding);
394   }
395   template <> long double as<long double>() const {
396     return mpfr_get_ld(value, mpfr_rounding);
397   }
398 
399   void dump(const char *msg) const { mpfr_printf("%s%.128Rf\n", msg, value); }
400 
401   // Return the ULP (units-in-the-last-place) difference between the
402   // stored MPFR and a floating point number.
403   //
404   // We define ULP difference as follows:
405   //   If exponents of this value and the |input| are same, then:
406   //     ULP(this_value, input) = abs(this_value - input) / eps(input)
407   //   else:
408   //     max = max(abs(this_value), abs(input))
409   //     min = min(abs(this_value), abs(input))
410   //     maxExponent = exponent(max)
411   //     ULP(this_value, input) = (max - 2^maxExponent) / eps(max) +
412   //                              (2^maxExponent - min) / eps(min)
413   //
414   // Remarks:
415   // 1. A ULP of 0.0 will imply that the value is correctly rounded.
416   // 2. We expect that this value and the value to be compared (the [input]
417   //    argument) are reasonable close, and we will provide an upper bound
418   //    of ULP value for testing.  Morever, most of the fractional parts of
419   //    ULP value do not matter much, so using double as the return type
420   //    should be good enough.
421   // 3. For close enough values (values which don't diff in their exponent by
422   //    not more than 1), a ULP difference of N indicates a bit distance
423   //    of N between this number and [input].
424   // 4. A values of +0.0 and -0.0 are treated as equal.
425   template <typename T>
426   cpp::EnableIfType<cpp::IsFloatingPointType<T>::Value, double> ulp(T input) {
427     T thisAsT = as<T>();
428     if (thisAsT == input)
429       return T(0.0);
430 
431     int thisExponent = fputil::FPBits<T>(thisAsT).get_exponent();
432     int inputExponent = fputil::FPBits<T>(input).get_exponent();
433     // Adjust the exponents for denormal numbers.
434     if (fputil::FPBits<T>(thisAsT).get_unbiased_exponent() == 0)
435       ++thisExponent;
436     if (fputil::FPBits<T>(input).get_unbiased_exponent() == 0)
437       ++inputExponent;
438 
439     if (thisAsT * input < 0 || thisExponent == inputExponent) {
440       MPFRNumber inputMPFR(input);
441       mpfr_sub(inputMPFR.value, value, inputMPFR.value, MPFR_RNDN);
442       mpfr_abs(inputMPFR.value, inputMPFR.value, MPFR_RNDN);
443       mpfr_mul_2si(inputMPFR.value, inputMPFR.value,
444                    -thisExponent + int(fputil::MantissaWidth<T>::VALUE),
445                    MPFR_RNDN);
446       return inputMPFR.as<double>();
447     }
448 
449     // If the control reaches here, it means that this number and input are
450     // of the same sign but different exponent. In such a case, ULP error is
451     // calculated as sum of two parts.
452     thisAsT = std::abs(thisAsT);
453     input = std::abs(input);
454     T min = thisAsT > input ? input : thisAsT;
455     T max = thisAsT > input ? thisAsT : input;
456     int minExponent = fputil::FPBits<T>(min).get_exponent();
457     int maxExponent = fputil::FPBits<T>(max).get_exponent();
458     // Adjust the exponents for denormal numbers.
459     if (fputil::FPBits<T>(min).get_unbiased_exponent() == 0)
460       ++minExponent;
461     if (fputil::FPBits<T>(max).get_unbiased_exponent() == 0)
462       ++maxExponent;
463 
464     MPFRNumber minMPFR(min);
465     MPFRNumber maxMPFR(max);
466 
467     MPFRNumber pivot(uint32_t(1));
468     mpfr_mul_2si(pivot.value, pivot.value, maxExponent, MPFR_RNDN);
469 
470     mpfr_sub(minMPFR.value, pivot.value, minMPFR.value, MPFR_RNDN);
471     mpfr_mul_2si(minMPFR.value, minMPFR.value,
472                  -minExponent + int(fputil::MantissaWidth<T>::VALUE),
473                  MPFR_RNDN);
474 
475     mpfr_sub(maxMPFR.value, maxMPFR.value, pivot.value, MPFR_RNDN);
476     mpfr_mul_2si(maxMPFR.value, maxMPFR.value,
477                  -maxExponent + int(fputil::MantissaWidth<T>::VALUE),
478                  MPFR_RNDN);
479 
480     mpfr_add(minMPFR.value, minMPFR.value, maxMPFR.value, MPFR_RNDN);
481     return minMPFR.as<double>();
482   }
483 };
484 
485 namespace internal {
486 
487 template <typename InputType>
488 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
489 unary_operation(Operation op, InputType input, unsigned int precision,
490                 RoundingMode rounding) {
491   MPFRNumber mpfrInput(input, precision, rounding);
492   switch (op) {
493   case Operation::Abs:
494     return mpfrInput.abs();
495   case Operation::Ceil:
496     return mpfrInput.ceil();
497   case Operation::Cos:
498     return mpfrInput.cos();
499   case Operation::Exp:
500     return mpfrInput.exp();
501   case Operation::Exp2:
502     return mpfrInput.exp2();
503   case Operation::Expm1:
504     return mpfrInput.expm1();
505   case Operation::Floor:
506     return mpfrInput.floor();
507   case Operation::Log:
508     return mpfrInput.log();
509   case Operation::Log2:
510     return mpfrInput.log2();
511   case Operation::Log10:
512     return mpfrInput.log10();
513   case Operation::Mod2PI:
514     return mpfrInput.mod_2pi();
515   case Operation::ModPIOver2:
516     return mpfrInput.mod_pi_over_2();
517   case Operation::ModPIOver4:
518     return mpfrInput.mod_pi_over_4();
519   case Operation::Round:
520     return mpfrInput.round();
521   case Operation::Sin:
522     return mpfrInput.sin();
523   case Operation::Sqrt:
524     return mpfrInput.sqrt();
525   case Operation::Tan:
526     return mpfrInput.tan();
527   case Operation::Trunc:
528     return mpfrInput.trunc();
529   default:
530     __builtin_unreachable();
531   }
532 }
533 
534 template <typename InputType>
535 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
536 unary_operation_two_outputs(Operation op, InputType input, int &output,
537                             unsigned int precision, RoundingMode rounding) {
538   MPFRNumber mpfrInput(input, precision, rounding);
539   switch (op) {
540   case Operation::Frexp:
541     return mpfrInput.frexp(output);
542   default:
543     __builtin_unreachable();
544   }
545 }
546 
547 template <typename InputType>
548 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
549 binary_operation_one_output(Operation op, InputType x, InputType y,
550                             unsigned int precision, RoundingMode rounding) {
551   MPFRNumber inputX(x, precision, rounding);
552   MPFRNumber inputY(y, precision, rounding);
553   switch (op) {
554   case Operation::Hypot:
555     return inputX.hypot(inputY);
556   default:
557     __builtin_unreachable();
558   }
559 }
560 
561 template <typename InputType>
562 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
563 binary_operation_two_outputs(Operation op, InputType x, InputType y,
564                              int &output, unsigned int precision,
565                              RoundingMode rounding) {
566   MPFRNumber inputX(x, precision, rounding);
567   MPFRNumber inputY(y, precision, rounding);
568   switch (op) {
569   case Operation::RemQuo:
570     return inputX.remquo(inputY, output);
571   default:
572     __builtin_unreachable();
573   }
574 }
575 
576 template <typename InputType>
577 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
578 ternary_operation_one_output(Operation op, InputType x, InputType y,
579                              InputType z, unsigned int precision,
580                              RoundingMode rounding) {
581   // For FMA function, we just need to compare with the mpfr_fma with the same
582   // precision as InputType.  Using higher precision as the intermediate results
583   // to compare might incorrectly fail due to double-rounding errors.
584   MPFRNumber inputX(x, precision, rounding);
585   MPFRNumber inputY(y, precision, rounding);
586   MPFRNumber inputZ(z, precision, rounding);
587   switch (op) {
588   case Operation::Fma:
589     return inputX.fma(inputY, inputZ);
590   default:
591     __builtin_unreachable();
592   }
593 }
594 
595 // Remark: For all the explain_*_error functions, we will use std::stringstream
596 // to build the complete error messages before sending it to the outstream `OS`
597 // once at the end.  This will stop the error messages from interleaving when
598 // the tests are running concurrently.
599 template <typename T>
600 void explain_unary_operation_single_output_error(Operation op, T input,
601                                                  T matchValue,
602                                                  double ulp_tolerance,
603                                                  RoundingMode rounding,
604                                                  testutils::StreamWrapper &OS) {
605   unsigned int precision = get_precision<T>(ulp_tolerance);
606   MPFRNumber mpfrInput(input, precision);
607   MPFRNumber mpfr_result;
608   mpfr_result = unary_operation(op, input, precision, rounding);
609   MPFRNumber mpfrMatchValue(matchValue);
610   std::stringstream ss;
611   ss << "Match value not within tolerance value of MPFR result:\n"
612      << "  Input decimal: " << mpfrInput.str() << '\n';
613   __llvm_libc::fputil::testing::describeValue("     Input bits: ", input, ss);
614   ss << '\n' << "  Match decimal: " << mpfrMatchValue.str() << '\n';
615   __llvm_libc::fputil::testing::describeValue("     Match bits: ", matchValue,
616                                               ss);
617   ss << '\n' << "    MPFR result: " << mpfr_result.str() << '\n';
618   __llvm_libc::fputil::testing::describeValue(
619       "   MPFR rounded: ", mpfr_result.as<T>(), ss);
620   ss << '\n';
621   ss << "      ULP error: " << std::to_string(mpfr_result.ulp(matchValue))
622      << '\n';
623   OS << ss.str();
624 }
625 
626 template void
627 explain_unary_operation_single_output_error<float>(Operation op, float, float,
628                                                    double, RoundingMode,
629                                                    testutils::StreamWrapper &);
630 template void explain_unary_operation_single_output_error<double>(
631     Operation op, double, double, double, RoundingMode,
632     testutils::StreamWrapper &);
633 template void explain_unary_operation_single_output_error<long double>(
634     Operation op, long double, long double, double, RoundingMode,
635     testutils::StreamWrapper &);
636 
637 template <typename T>
638 void explain_unary_operation_two_outputs_error(
639     Operation op, T input, const BinaryOutput<T> &libc_result,
640     double ulp_tolerance, RoundingMode rounding, testutils::StreamWrapper &OS) {
641   unsigned int precision = get_precision<T>(ulp_tolerance);
642   MPFRNumber mpfrInput(input, precision);
643   int mpfrIntResult;
644   MPFRNumber mpfr_result = unary_operation_two_outputs(op, input, mpfrIntResult,
645                                                        precision, rounding);
646   std::stringstream ss;
647 
648   if (mpfrIntResult != libc_result.i) {
649     ss << "MPFR integral result: " << mpfrIntResult << '\n'
650        << "Libc integral result: " << libc_result.i << '\n';
651   } else {
652     ss << "Integral result from libc matches integral result from MPFR.\n";
653   }
654 
655   MPFRNumber mpfrMatchValue(libc_result.f);
656   ss << "Libc floating point result is not within tolerance value of the MPFR "
657      << "result.\n\n";
658 
659   ss << "            Input decimal: " << mpfrInput.str() << "\n\n";
660 
661   ss << "Libc floating point value: " << mpfrMatchValue.str() << '\n';
662   __llvm_libc::fputil::testing::describeValue(
663       " Libc floating point bits: ", libc_result.f, ss);
664   ss << "\n\n";
665 
666   ss << "              MPFR result: " << mpfr_result.str() << '\n';
667   __llvm_libc::fputil::testing::describeValue(
668       "             MPFR rounded: ", mpfr_result.as<T>(), ss);
669   ss << '\n'
670      << "                ULP error: "
671      << std::to_string(mpfr_result.ulp(libc_result.f)) << '\n';
672   OS << ss.str();
673 }
674 
675 template void explain_unary_operation_two_outputs_error<float>(
676     Operation, float, const BinaryOutput<float> &, double, RoundingMode,
677     testutils::StreamWrapper &);
678 template void explain_unary_operation_two_outputs_error<double>(
679     Operation, double, const BinaryOutput<double> &, double, RoundingMode,
680     testutils::StreamWrapper &);
681 template void explain_unary_operation_two_outputs_error<long double>(
682     Operation, long double, const BinaryOutput<long double> &, double,
683     RoundingMode, testutils::StreamWrapper &);
684 
685 template <typename T>
686 void explain_binary_operation_two_outputs_error(
687     Operation op, const BinaryInput<T> &input,
688     const BinaryOutput<T> &libc_result, double ulp_tolerance,
689     RoundingMode rounding, testutils::StreamWrapper &OS) {
690   unsigned int precision = get_precision<T>(ulp_tolerance);
691   MPFRNumber mpfrX(input.x, precision);
692   MPFRNumber mpfrY(input.y, precision);
693   int mpfrIntResult;
694   MPFRNumber mpfr_result = binary_operation_two_outputs(
695       op, input.x, input.y, mpfrIntResult, precision, rounding);
696   MPFRNumber mpfrMatchValue(libc_result.f);
697   std::stringstream ss;
698 
699   ss << "Input decimal: x: " << mpfrX.str() << " y: " << mpfrY.str() << '\n'
700      << "MPFR integral result: " << mpfrIntResult << '\n'
701      << "Libc integral result: " << libc_result.i << '\n'
702      << "Libc floating point result: " << mpfrMatchValue.str() << '\n'
703      << "               MPFR result: " << mpfr_result.str() << '\n';
704   __llvm_libc::fputil::testing::describeValue(
705       "Libc floating point result bits: ", libc_result.f, ss);
706   __llvm_libc::fputil::testing::describeValue(
707       "              MPFR rounded bits: ", mpfr_result.as<T>(), ss);
708   ss << "ULP error: " << std::to_string(mpfr_result.ulp(libc_result.f)) << '\n';
709   OS << ss.str();
710 }
711 
712 template void explain_binary_operation_two_outputs_error<float>(
713     Operation, const BinaryInput<float> &, const BinaryOutput<float> &, double,
714     RoundingMode, testutils::StreamWrapper &);
715 template void explain_binary_operation_two_outputs_error<double>(
716     Operation, const BinaryInput<double> &, const BinaryOutput<double> &,
717     double, RoundingMode, testutils::StreamWrapper &);
718 template void explain_binary_operation_two_outputs_error<long double>(
719     Operation, const BinaryInput<long double> &,
720     const BinaryOutput<long double> &, double, RoundingMode,
721     testutils::StreamWrapper &);
722 
723 template <typename T>
724 void explain_binary_operation_one_output_error(
725     Operation op, const BinaryInput<T> &input, T libc_result,
726     double ulp_tolerance, RoundingMode rounding, testutils::StreamWrapper &OS) {
727   unsigned int precision = get_precision<T>(ulp_tolerance);
728   MPFRNumber mpfrX(input.x, precision);
729   MPFRNumber mpfrY(input.y, precision);
730   FPBits<T> xbits(input.x);
731   FPBits<T> ybits(input.y);
732   MPFRNumber mpfr_result =
733       binary_operation_one_output(op, input.x, input.y, precision, rounding);
734   MPFRNumber mpfrMatchValue(libc_result);
735   std::stringstream ss;
736 
737   ss << "Input decimal: x: " << mpfrX.str() << " y: " << mpfrY.str() << '\n';
738   __llvm_libc::fputil::testing::describeValue("First input bits: ", input.x,
739                                               ss);
740   __llvm_libc::fputil::testing::describeValue("Second input bits: ", input.y,
741                                               ss);
742 
743   ss << "Libc result: " << mpfrMatchValue.str() << '\n'
744      << "MPFR result: " << mpfr_result.str() << '\n';
745   __llvm_libc::fputil::testing::describeValue(
746       "Libc floating point result bits: ", libc_result, ss);
747   __llvm_libc::fputil::testing::describeValue(
748       "              MPFR rounded bits: ", mpfr_result.as<T>(), ss);
749   ss << "ULP error: " << std::to_string(mpfr_result.ulp(libc_result)) << '\n';
750   OS << ss.str();
751 }
752 
753 template void explain_binary_operation_one_output_error<float>(
754     Operation, const BinaryInput<float> &, float, double, RoundingMode,
755     testutils::StreamWrapper &);
756 template void explain_binary_operation_one_output_error<double>(
757     Operation, const BinaryInput<double> &, double, double, RoundingMode,
758     testutils::StreamWrapper &);
759 template void explain_binary_operation_one_output_error<long double>(
760     Operation, const BinaryInput<long double> &, long double, double,
761     RoundingMode, testutils::StreamWrapper &);
762 
763 template <typename T>
764 void explain_ternary_operation_one_output_error(
765     Operation op, const TernaryInput<T> &input, T libc_result,
766     double ulp_tolerance, RoundingMode rounding, testutils::StreamWrapper &OS) {
767   unsigned int precision = get_precision<T>(ulp_tolerance);
768   MPFRNumber mpfrX(input.x, precision);
769   MPFRNumber mpfrY(input.y, precision);
770   MPFRNumber mpfrZ(input.z, precision);
771   FPBits<T> xbits(input.x);
772   FPBits<T> ybits(input.y);
773   FPBits<T> zbits(input.z);
774   MPFRNumber mpfr_result = ternary_operation_one_output(
775       op, input.x, input.y, input.z, precision, rounding);
776   MPFRNumber mpfrMatchValue(libc_result);
777   std::stringstream ss;
778 
779   ss << "Input decimal: x: " << mpfrX.str() << " y: " << mpfrY.str()
780      << " z: " << mpfrZ.str() << '\n';
781   __llvm_libc::fputil::testing::describeValue("First input bits: ", input.x,
782                                               ss);
783   __llvm_libc::fputil::testing::describeValue("Second input bits: ", input.y,
784                                               ss);
785   __llvm_libc::fputil::testing::describeValue("Third input bits: ", input.z,
786                                               ss);
787 
788   ss << "Libc result: " << mpfrMatchValue.str() << '\n'
789      << "MPFR result: " << mpfr_result.str() << '\n';
790   __llvm_libc::fputil::testing::describeValue(
791       "Libc floating point result bits: ", libc_result, ss);
792   __llvm_libc::fputil::testing::describeValue(
793       "              MPFR rounded bits: ", mpfr_result.as<T>(), ss);
794   ss << "ULP error: " << std::to_string(mpfr_result.ulp(libc_result)) << '\n';
795   OS << ss.str();
796 }
797 
798 template void explain_ternary_operation_one_output_error<float>(
799     Operation, const TernaryInput<float> &, float, double, RoundingMode,
800     testutils::StreamWrapper &);
801 template void explain_ternary_operation_one_output_error<double>(
802     Operation, const TernaryInput<double> &, double, double, RoundingMode,
803     testutils::StreamWrapper &);
804 template void explain_ternary_operation_one_output_error<long double>(
805     Operation, const TernaryInput<long double> &, long double, double,
806     RoundingMode, testutils::StreamWrapper &);
807 
808 template <typename T>
809 bool compare_unary_operation_single_output(Operation op, T input, T libc_result,
810                                            double ulp_tolerance,
811                                            RoundingMode rounding) {
812   unsigned int precision = get_precision<T>(ulp_tolerance);
813   MPFRNumber mpfr_result;
814   mpfr_result = unary_operation(op, input, precision, rounding);
815   double ulp = mpfr_result.ulp(libc_result);
816   return (ulp <= ulp_tolerance);
817 }
818 
819 template bool compare_unary_operation_single_output<float>(Operation, float,
820                                                            float, double,
821                                                            RoundingMode);
822 template bool compare_unary_operation_single_output<double>(Operation, double,
823                                                             double, double,
824                                                             RoundingMode);
825 template bool compare_unary_operation_single_output<long double>(
826     Operation, long double, long double, double, RoundingMode);
827 
828 template <typename T>
829 bool compare_unary_operation_two_outputs(Operation op, T input,
830                                          const BinaryOutput<T> &libc_result,
831                                          double ulp_tolerance,
832                                          RoundingMode rounding) {
833   int mpfrIntResult;
834   unsigned int precision = get_precision<T>(ulp_tolerance);
835   MPFRNumber mpfr_result = unary_operation_two_outputs(op, input, mpfrIntResult,
836                                                        precision, rounding);
837   double ulp = mpfr_result.ulp(libc_result.f);
838 
839   if (mpfrIntResult != libc_result.i)
840     return false;
841 
842   return (ulp <= ulp_tolerance);
843 }
844 
845 template bool compare_unary_operation_two_outputs<float>(
846     Operation, float, const BinaryOutput<float> &, double, RoundingMode);
847 template bool compare_unary_operation_two_outputs<double>(
848     Operation, double, const BinaryOutput<double> &, double, RoundingMode);
849 template bool compare_unary_operation_two_outputs<long double>(
850     Operation, long double, const BinaryOutput<long double> &, double,
851     RoundingMode);
852 
853 template <typename T>
854 bool compare_binary_operation_two_outputs(Operation op,
855                                           const BinaryInput<T> &input,
856                                           const BinaryOutput<T> &libc_result,
857                                           double ulp_tolerance,
858                                           RoundingMode rounding) {
859   int mpfrIntResult;
860   unsigned int precision = get_precision<T>(ulp_tolerance);
861   MPFRNumber mpfr_result = binary_operation_two_outputs(
862       op, input.x, input.y, mpfrIntResult, precision, rounding);
863   double ulp = mpfr_result.ulp(libc_result.f);
864 
865   if (mpfrIntResult != libc_result.i) {
866     if (op == Operation::RemQuo) {
867       if ((0x7 & mpfrIntResult) != (0x7 & libc_result.i))
868         return false;
869     } else {
870       return false;
871     }
872   }
873 
874   return (ulp <= ulp_tolerance);
875 }
876 
877 template bool compare_binary_operation_two_outputs<float>(
878     Operation, const BinaryInput<float> &, const BinaryOutput<float> &, double,
879     RoundingMode);
880 template bool compare_binary_operation_two_outputs<double>(
881     Operation, const BinaryInput<double> &, const BinaryOutput<double> &,
882     double, RoundingMode);
883 template bool compare_binary_operation_two_outputs<long double>(
884     Operation, const BinaryInput<long double> &,
885     const BinaryOutput<long double> &, double, RoundingMode);
886 
887 template <typename T>
888 bool compare_binary_operation_one_output(Operation op,
889                                          const BinaryInput<T> &input,
890                                          T libc_result, double ulp_tolerance,
891                                          RoundingMode rounding) {
892   unsigned int precision = get_precision<T>(ulp_tolerance);
893   MPFRNumber mpfr_result =
894       binary_operation_one_output(op, input.x, input.y, precision, rounding);
895   double ulp = mpfr_result.ulp(libc_result);
896 
897   return (ulp <= ulp_tolerance);
898 }
899 
900 template bool compare_binary_operation_one_output<float>(
901     Operation, const BinaryInput<float> &, float, double, RoundingMode);
902 template bool compare_binary_operation_one_output<double>(
903     Operation, const BinaryInput<double> &, double, double, RoundingMode);
904 template bool compare_binary_operation_one_output<long double>(
905     Operation, const BinaryInput<long double> &, long double, double,
906     RoundingMode);
907 
908 template <typename T>
909 bool compare_ternary_operation_one_output(Operation op,
910                                           const TernaryInput<T> &input,
911                                           T libc_result, double ulp_tolerance,
912                                           RoundingMode rounding) {
913   unsigned int precision = get_precision<T>(ulp_tolerance);
914   MPFRNumber mpfr_result = ternary_operation_one_output(
915       op, input.x, input.y, input.z, precision, rounding);
916   double ulp = mpfr_result.ulp(libc_result);
917 
918   return (ulp <= ulp_tolerance);
919 }
920 
921 template bool compare_ternary_operation_one_output<float>(
922     Operation, const TernaryInput<float> &, float, double, RoundingMode);
923 template bool compare_ternary_operation_one_output<double>(
924     Operation, const TernaryInput<double> &, double, double, RoundingMode);
925 template bool compare_ternary_operation_one_output<long double>(
926     Operation, const TernaryInput<long double> &, long double, double,
927     RoundingMode);
928 
929 } // namespace internal
930 
931 template <typename T> bool round_to_long(T x, long &result) {
932   MPFRNumber mpfr(x);
933   return mpfr.round_to_long(result);
934 }
935 
936 template bool round_to_long<float>(float, long &);
937 template bool round_to_long<double>(double, long &);
938 template bool round_to_long<long double>(long double, long &);
939 
940 template <typename T> bool round_to_long(T x, RoundingMode mode, long &result) {
941   MPFRNumber mpfr(x);
942   return mpfr.round_to_long(get_mpfr_rounding_mode(mode), result);
943 }
944 
945 template bool round_to_long<float>(float, RoundingMode, long &);
946 template bool round_to_long<double>(double, RoundingMode, long &);
947 template bool round_to_long<long double>(long double, RoundingMode, long &);
948 
949 template <typename T> T round(T x, RoundingMode mode) {
950   MPFRNumber mpfr(x);
951   MPFRNumber result = mpfr.rint(get_mpfr_rounding_mode(mode));
952   return result.as<T>();
953 }
954 
955 template float round<float>(float, RoundingMode);
956 template double round<double>(double, RoundingMode);
957 template long double round<long double>(long double, RoundingMode);
958 
959 } // namespace mpfr
960 } // namespace testing
961 } // namespace __llvm_libc
962