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 remquo(const MPFRNumber &divisor, int &quotient) {
271     MPFRNumber remainder(*this);
272     long q;
273     mpfr_remquo(remainder.value, &q, value, divisor.value, mpfr_rounding);
274     quotient = q;
275     return remainder;
276   }
277 
278   MPFRNumber round() const {
279     MPFRNumber result(*this);
280     mpfr_round(result.value, value);
281     return result;
282   }
283 
284   bool round_to_long(long &result) const {
285     // We first calculate the rounded value. This way, when converting
286     // to long using mpfr_get_si, the rounding direction of MPFR_RNDN
287     // (or any other rounding mode), does not have an influence.
288     MPFRNumber roundedValue = round();
289     mpfr_clear_erangeflag();
290     result = mpfr_get_si(roundedValue.value, MPFR_RNDN);
291     return mpfr_erangeflag_p();
292   }
293 
294   bool round_to_long(mpfr_rnd_t rnd, long &result) const {
295     MPFRNumber rint_result(*this);
296     mpfr_rint(rint_result.value, value, rnd);
297     return rint_result.round_to_long(result);
298   }
299 
300   MPFRNumber rint(mpfr_rnd_t rnd) const {
301     MPFRNumber result(*this);
302     mpfr_rint(result.value, value, rnd);
303     return result;
304   }
305 
306   MPFRNumber mod_2pi() const {
307     MPFRNumber result(0.0, 1280);
308     MPFRNumber _2pi(0.0, 1280);
309     mpfr_const_pi(_2pi.value, MPFR_RNDN);
310     mpfr_mul_si(_2pi.value, _2pi.value, 2, MPFR_RNDN);
311     mpfr_fmod(result.value, value, _2pi.value, MPFR_RNDN);
312     return result;
313   }
314 
315   MPFRNumber mod_pi_over_2() const {
316     MPFRNumber result(0.0, 1280);
317     MPFRNumber pi_over_2(0.0, 1280);
318     mpfr_const_pi(pi_over_2.value, MPFR_RNDN);
319     mpfr_mul_d(pi_over_2.value, pi_over_2.value, 0.5, MPFR_RNDN);
320     mpfr_fmod(result.value, value, pi_over_2.value, MPFR_RNDN);
321     return result;
322   }
323 
324   MPFRNumber mod_pi_over_4() const {
325     MPFRNumber result(0.0, 1280);
326     MPFRNumber pi_over_4(0.0, 1280);
327     mpfr_const_pi(pi_over_4.value, MPFR_RNDN);
328     mpfr_mul_d(pi_over_4.value, pi_over_4.value, 0.25, MPFR_RNDN);
329     mpfr_fmod(result.value, value, pi_over_4.value, MPFR_RNDN);
330     return result;
331   }
332 
333   MPFRNumber sin() const {
334     MPFRNumber result(*this);
335     mpfr_sin(result.value, value, mpfr_rounding);
336     return result;
337   }
338 
339   MPFRNumber sqrt() const {
340     MPFRNumber result(*this);
341     mpfr_sqrt(result.value, value, mpfr_rounding);
342     return result;
343   }
344 
345   MPFRNumber tan() const {
346     MPFRNumber result(*this);
347     mpfr_tan(result.value, value, mpfr_rounding);
348     return result;
349   }
350 
351   MPFRNumber trunc() const {
352     MPFRNumber result(*this);
353     mpfr_trunc(result.value, value);
354     return result;
355   }
356 
357   MPFRNumber fma(const MPFRNumber &b, const MPFRNumber &c) {
358     MPFRNumber result(*this);
359     mpfr_fma(result.value, value, b.value, c.value, mpfr_rounding);
360     return result;
361   }
362 
363   std::string str() const {
364     // 200 bytes should be more than sufficient to hold a 100-digit number
365     // plus additional bytes for the decimal point, '-' sign etc.
366     constexpr size_t printBufSize = 200;
367     char buffer[printBufSize];
368     mpfr_snprintf(buffer, printBufSize, "%100.50Rf", value);
369     cpp::StringView view(buffer);
370     view = view.trim(' ');
371     return std::string(view.data());
372   }
373 
374   // These functions are useful for debugging.
375   template <typename T> T as() const;
376 
377   template <> float as<float>() const {
378     return mpfr_get_flt(value, mpfr_rounding);
379   }
380   template <> double as<double>() const {
381     return mpfr_get_d(value, mpfr_rounding);
382   }
383   template <> long double as<long double>() const {
384     return mpfr_get_ld(value, mpfr_rounding);
385   }
386 
387   void dump(const char *msg) const { mpfr_printf("%s%.128Rf\n", msg, value); }
388 
389   // Return the ULP (units-in-the-last-place) difference between the
390   // stored MPFR and a floating point number.
391   //
392   // We define ULP difference as follows:
393   //   If exponents of this value and the |input| are same, then:
394   //     ULP(this_value, input) = abs(this_value - input) / eps(input)
395   //   else:
396   //     max = max(abs(this_value), abs(input))
397   //     min = min(abs(this_value), abs(input))
398   //     maxExponent = exponent(max)
399   //     ULP(this_value, input) = (max - 2^maxExponent) / eps(max) +
400   //                              (2^maxExponent - min) / eps(min)
401   //
402   // Remarks:
403   // 1. A ULP of 0.0 will imply that the value is correctly rounded.
404   // 2. We expect that this value and the value to be compared (the [input]
405   //    argument) are reasonable close, and we will provide an upper bound
406   //    of ULP value for testing.  Morever, most of the fractional parts of
407   //    ULP value do not matter much, so using double as the return type
408   //    should be good enough.
409   // 3. For close enough values (values which don't diff in their exponent by
410   //    not more than 1), a ULP difference of N indicates a bit distance
411   //    of N between this number and [input].
412   // 4. A values of +0.0 and -0.0 are treated as equal.
413   template <typename T>
414   cpp::EnableIfType<cpp::IsFloatingPointType<T>::Value, double> ulp(T input) {
415     T thisAsT = as<T>();
416     if (thisAsT == input)
417       return T(0.0);
418 
419     int thisExponent = fputil::FPBits<T>(thisAsT).get_exponent();
420     int inputExponent = fputil::FPBits<T>(input).get_exponent();
421     // Adjust the exponents for denormal numbers.
422     if (fputil::FPBits<T>(thisAsT).get_unbiased_exponent() == 0)
423       ++thisExponent;
424     if (fputil::FPBits<T>(input).get_unbiased_exponent() == 0)
425       ++inputExponent;
426 
427     if (thisAsT * input < 0 || thisExponent == inputExponent) {
428       MPFRNumber inputMPFR(input);
429       mpfr_sub(inputMPFR.value, value, inputMPFR.value, MPFR_RNDN);
430       mpfr_abs(inputMPFR.value, inputMPFR.value, MPFR_RNDN);
431       mpfr_mul_2si(inputMPFR.value, inputMPFR.value,
432                    -thisExponent + int(fputil::MantissaWidth<T>::VALUE),
433                    MPFR_RNDN);
434       return inputMPFR.as<double>();
435     }
436 
437     // If the control reaches here, it means that this number and input are
438     // of the same sign but different exponent. In such a case, ULP error is
439     // calculated as sum of two parts.
440     thisAsT = std::abs(thisAsT);
441     input = std::abs(input);
442     T min = thisAsT > input ? input : thisAsT;
443     T max = thisAsT > input ? thisAsT : input;
444     int minExponent = fputil::FPBits<T>(min).get_exponent();
445     int maxExponent = fputil::FPBits<T>(max).get_exponent();
446     // Adjust the exponents for denormal numbers.
447     if (fputil::FPBits<T>(min).get_unbiased_exponent() == 0)
448       ++minExponent;
449     if (fputil::FPBits<T>(max).get_unbiased_exponent() == 0)
450       ++maxExponent;
451 
452     MPFRNumber minMPFR(min);
453     MPFRNumber maxMPFR(max);
454 
455     MPFRNumber pivot(uint32_t(1));
456     mpfr_mul_2si(pivot.value, pivot.value, maxExponent, MPFR_RNDN);
457 
458     mpfr_sub(minMPFR.value, pivot.value, minMPFR.value, MPFR_RNDN);
459     mpfr_mul_2si(minMPFR.value, minMPFR.value,
460                  -minExponent + int(fputil::MantissaWidth<T>::VALUE),
461                  MPFR_RNDN);
462 
463     mpfr_sub(maxMPFR.value, maxMPFR.value, pivot.value, MPFR_RNDN);
464     mpfr_mul_2si(maxMPFR.value, maxMPFR.value,
465                  -maxExponent + int(fputil::MantissaWidth<T>::VALUE),
466                  MPFR_RNDN);
467 
468     mpfr_add(minMPFR.value, minMPFR.value, maxMPFR.value, MPFR_RNDN);
469     return minMPFR.as<double>();
470   }
471 };
472 
473 namespace internal {
474 
475 template <typename InputType>
476 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
477 unary_operation(Operation op, InputType input, unsigned int precision,
478                 RoundingMode rounding) {
479   MPFRNumber mpfrInput(input, precision, rounding);
480   switch (op) {
481   case Operation::Abs:
482     return mpfrInput.abs();
483   case Operation::Ceil:
484     return mpfrInput.ceil();
485   case Operation::Cos:
486     return mpfrInput.cos();
487   case Operation::Exp:
488     return mpfrInput.exp();
489   case Operation::Exp2:
490     return mpfrInput.exp2();
491   case Operation::Expm1:
492     return mpfrInput.expm1();
493   case Operation::Floor:
494     return mpfrInput.floor();
495   case Operation::Log:
496     return mpfrInput.log();
497   case Operation::Mod2PI:
498     return mpfrInput.mod_2pi();
499   case Operation::ModPIOver2:
500     return mpfrInput.mod_pi_over_2();
501   case Operation::ModPIOver4:
502     return mpfrInput.mod_pi_over_4();
503   case Operation::Round:
504     return mpfrInput.round();
505   case Operation::Sin:
506     return mpfrInput.sin();
507   case Operation::Sqrt:
508     return mpfrInput.sqrt();
509   case Operation::Tan:
510     return mpfrInput.tan();
511   case Operation::Trunc:
512     return mpfrInput.trunc();
513   default:
514     __builtin_unreachable();
515   }
516 }
517 
518 template <typename InputType>
519 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
520 unary_operation_two_outputs(Operation op, InputType input, int &output,
521                             unsigned int precision, RoundingMode rounding) {
522   MPFRNumber mpfrInput(input, precision, rounding);
523   switch (op) {
524   case Operation::Frexp:
525     return mpfrInput.frexp(output);
526   default:
527     __builtin_unreachable();
528   }
529 }
530 
531 template <typename InputType>
532 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
533 binary_operation_one_output(Operation op, InputType x, InputType y,
534                             unsigned int precision, RoundingMode rounding) {
535   MPFRNumber inputX(x, precision, rounding);
536   MPFRNumber inputY(y, precision, rounding);
537   switch (op) {
538   case Operation::Hypot:
539     return inputX.hypot(inputY);
540   default:
541     __builtin_unreachable();
542   }
543 }
544 
545 template <typename InputType>
546 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
547 binary_operation_two_outputs(Operation op, InputType x, InputType y,
548                              int &output, unsigned int precision,
549                              RoundingMode rounding) {
550   MPFRNumber inputX(x, precision, rounding);
551   MPFRNumber inputY(y, precision, rounding);
552   switch (op) {
553   case Operation::RemQuo:
554     return inputX.remquo(inputY, output);
555   default:
556     __builtin_unreachable();
557   }
558 }
559 
560 template <typename InputType>
561 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
562 ternary_operation_one_output(Operation op, InputType x, InputType y,
563                              InputType z, unsigned int precision,
564                              RoundingMode rounding) {
565   // For FMA function, we just need to compare with the mpfr_fma with the same
566   // precision as InputType.  Using higher precision as the intermediate results
567   // to compare might incorrectly fail due to double-rounding errors.
568   MPFRNumber inputX(x, precision, rounding);
569   MPFRNumber inputY(y, precision, rounding);
570   MPFRNumber inputZ(z, precision, rounding);
571   switch (op) {
572   case Operation::Fma:
573     return inputX.fma(inputY, inputZ);
574   default:
575     __builtin_unreachable();
576   }
577 }
578 
579 // Remark: For all the explain_*_error functions, we will use std::stringstream
580 // to build the complete error messages before sending it to the outstream `OS`
581 // once at the end.  This will stop the error messages from interleaving when
582 // the tests are running concurrently.
583 template <typename T>
584 void explain_unary_operation_single_output_error(Operation op, T input,
585                                                  T matchValue,
586                                                  double ulp_tolerance,
587                                                  RoundingMode rounding,
588                                                  testutils::StreamWrapper &OS) {
589   unsigned int precision = get_precision<T>(ulp_tolerance);
590   MPFRNumber mpfrInput(input, precision);
591   MPFRNumber mpfr_result;
592   mpfr_result = unary_operation(op, input, precision, rounding);
593   MPFRNumber mpfrMatchValue(matchValue);
594   std::stringstream ss;
595   ss << "Match value not within tolerance value of MPFR result:\n"
596      << "  Input decimal: " << mpfrInput.str() << '\n';
597   __llvm_libc::fputil::testing::describeValue("     Input bits: ", input, ss);
598   ss << '\n' << "  Match decimal: " << mpfrMatchValue.str() << '\n';
599   __llvm_libc::fputil::testing::describeValue("     Match bits: ", matchValue,
600                                               ss);
601   ss << '\n' << "    MPFR result: " << mpfr_result.str() << '\n';
602   __llvm_libc::fputil::testing::describeValue(
603       "   MPFR rounded: ", mpfr_result.as<T>(), ss);
604   ss << '\n';
605   ss << "      ULP error: " << std::to_string(mpfr_result.ulp(matchValue))
606      << '\n';
607   OS << ss.str();
608 }
609 
610 template void
611 explain_unary_operation_single_output_error<float>(Operation op, float, float,
612                                                    double, RoundingMode,
613                                                    testutils::StreamWrapper &);
614 template void explain_unary_operation_single_output_error<double>(
615     Operation op, double, double, double, RoundingMode,
616     testutils::StreamWrapper &);
617 template void explain_unary_operation_single_output_error<long double>(
618     Operation op, long double, long double, double, RoundingMode,
619     testutils::StreamWrapper &);
620 
621 template <typename T>
622 void explain_unary_operation_two_outputs_error(
623     Operation op, T input, const BinaryOutput<T> &libc_result,
624     double ulp_tolerance, RoundingMode rounding, testutils::StreamWrapper &OS) {
625   unsigned int precision = get_precision<T>(ulp_tolerance);
626   MPFRNumber mpfrInput(input, precision);
627   int mpfrIntResult;
628   MPFRNumber mpfr_result = unary_operation_two_outputs(op, input, mpfrIntResult,
629                                                        precision, rounding);
630   std::stringstream ss;
631 
632   if (mpfrIntResult != libc_result.i) {
633     ss << "MPFR integral result: " << mpfrIntResult << '\n'
634        << "Libc integral result: " << libc_result.i << '\n';
635   } else {
636     ss << "Integral result from libc matches integral result from MPFR.\n";
637   }
638 
639   MPFRNumber mpfrMatchValue(libc_result.f);
640   ss << "Libc floating point result is not within tolerance value of the MPFR "
641      << "result.\n\n";
642 
643   ss << "            Input decimal: " << mpfrInput.str() << "\n\n";
644 
645   ss << "Libc floating point value: " << mpfrMatchValue.str() << '\n';
646   __llvm_libc::fputil::testing::describeValue(
647       " Libc floating point bits: ", libc_result.f, ss);
648   ss << "\n\n";
649 
650   ss << "              MPFR result: " << mpfr_result.str() << '\n';
651   __llvm_libc::fputil::testing::describeValue(
652       "             MPFR rounded: ", mpfr_result.as<T>(), ss);
653   ss << '\n'
654      << "                ULP error: "
655      << std::to_string(mpfr_result.ulp(libc_result.f)) << '\n';
656   OS << ss.str();
657 }
658 
659 template void explain_unary_operation_two_outputs_error<float>(
660     Operation, float, const BinaryOutput<float> &, double, RoundingMode,
661     testutils::StreamWrapper &);
662 template void explain_unary_operation_two_outputs_error<double>(
663     Operation, double, const BinaryOutput<double> &, double, RoundingMode,
664     testutils::StreamWrapper &);
665 template void explain_unary_operation_two_outputs_error<long double>(
666     Operation, long double, const BinaryOutput<long double> &, double,
667     RoundingMode, testutils::StreamWrapper &);
668 
669 template <typename T>
670 void explain_binary_operation_two_outputs_error(
671     Operation op, const BinaryInput<T> &input,
672     const BinaryOutput<T> &libc_result, double ulp_tolerance,
673     RoundingMode rounding, testutils::StreamWrapper &OS) {
674   unsigned int precision = get_precision<T>(ulp_tolerance);
675   MPFRNumber mpfrX(input.x, precision);
676   MPFRNumber mpfrY(input.y, precision);
677   int mpfrIntResult;
678   MPFRNumber mpfr_result = binary_operation_two_outputs(
679       op, input.x, input.y, mpfrIntResult, precision, rounding);
680   MPFRNumber mpfrMatchValue(libc_result.f);
681   std::stringstream ss;
682 
683   ss << "Input decimal: x: " << mpfrX.str() << " y: " << mpfrY.str() << '\n'
684      << "MPFR integral result: " << mpfrIntResult << '\n'
685      << "Libc integral result: " << libc_result.i << '\n'
686      << "Libc floating point result: " << mpfrMatchValue.str() << '\n'
687      << "               MPFR result: " << mpfr_result.str() << '\n';
688   __llvm_libc::fputil::testing::describeValue(
689       "Libc floating point result bits: ", libc_result.f, ss);
690   __llvm_libc::fputil::testing::describeValue(
691       "              MPFR rounded bits: ", mpfr_result.as<T>(), ss);
692   ss << "ULP error: " << std::to_string(mpfr_result.ulp(libc_result.f)) << '\n';
693   OS << ss.str();
694 }
695 
696 template void explain_binary_operation_two_outputs_error<float>(
697     Operation, const BinaryInput<float> &, const BinaryOutput<float> &, double,
698     RoundingMode, testutils::StreamWrapper &);
699 template void explain_binary_operation_two_outputs_error<double>(
700     Operation, const BinaryInput<double> &, const BinaryOutput<double> &,
701     double, RoundingMode, testutils::StreamWrapper &);
702 template void explain_binary_operation_two_outputs_error<long double>(
703     Operation, const BinaryInput<long double> &,
704     const BinaryOutput<long double> &, double, RoundingMode,
705     testutils::StreamWrapper &);
706 
707 template <typename T>
708 void explain_binary_operation_one_output_error(
709     Operation op, const BinaryInput<T> &input, T libc_result,
710     double ulp_tolerance, RoundingMode rounding, testutils::StreamWrapper &OS) {
711   unsigned int precision = get_precision<T>(ulp_tolerance);
712   MPFRNumber mpfrX(input.x, precision);
713   MPFRNumber mpfrY(input.y, precision);
714   FPBits<T> xbits(input.x);
715   FPBits<T> ybits(input.y);
716   MPFRNumber mpfr_result =
717       binary_operation_one_output(op, input.x, input.y, precision, rounding);
718   MPFRNumber mpfrMatchValue(libc_result);
719   std::stringstream ss;
720 
721   ss << "Input decimal: x: " << mpfrX.str() << " y: " << mpfrY.str() << '\n';
722   __llvm_libc::fputil::testing::describeValue("First input bits: ", input.x,
723                                               ss);
724   __llvm_libc::fputil::testing::describeValue("Second input bits: ", input.y,
725                                               ss);
726 
727   ss << "Libc result: " << mpfrMatchValue.str() << '\n'
728      << "MPFR result: " << mpfr_result.str() << '\n';
729   __llvm_libc::fputil::testing::describeValue(
730       "Libc floating point result bits: ", libc_result, ss);
731   __llvm_libc::fputil::testing::describeValue(
732       "              MPFR rounded bits: ", mpfr_result.as<T>(), ss);
733   ss << "ULP error: " << std::to_string(mpfr_result.ulp(libc_result)) << '\n';
734   OS << ss.str();
735 }
736 
737 template void explain_binary_operation_one_output_error<float>(
738     Operation, const BinaryInput<float> &, float, double, RoundingMode,
739     testutils::StreamWrapper &);
740 template void explain_binary_operation_one_output_error<double>(
741     Operation, const BinaryInput<double> &, double, double, RoundingMode,
742     testutils::StreamWrapper &);
743 template void explain_binary_operation_one_output_error<long double>(
744     Operation, const BinaryInput<long double> &, long double, double,
745     RoundingMode, testutils::StreamWrapper &);
746 
747 template <typename T>
748 void explain_ternary_operation_one_output_error(
749     Operation op, const TernaryInput<T> &input, T libc_result,
750     double ulp_tolerance, RoundingMode rounding, testutils::StreamWrapper &OS) {
751   unsigned int precision = get_precision<T>(ulp_tolerance);
752   MPFRNumber mpfrX(input.x, precision);
753   MPFRNumber mpfrY(input.y, precision);
754   MPFRNumber mpfrZ(input.z, precision);
755   FPBits<T> xbits(input.x);
756   FPBits<T> ybits(input.y);
757   FPBits<T> zbits(input.z);
758   MPFRNumber mpfr_result = ternary_operation_one_output(
759       op, input.x, input.y, input.z, precision, rounding);
760   MPFRNumber mpfrMatchValue(libc_result);
761   std::stringstream ss;
762 
763   ss << "Input decimal: x: " << mpfrX.str() << " y: " << mpfrY.str()
764      << " z: " << mpfrZ.str() << '\n';
765   __llvm_libc::fputil::testing::describeValue("First input bits: ", input.x,
766                                               ss);
767   __llvm_libc::fputil::testing::describeValue("Second input bits: ", input.y,
768                                               ss);
769   __llvm_libc::fputil::testing::describeValue("Third input bits: ", input.z,
770                                               ss);
771 
772   ss << "Libc result: " << mpfrMatchValue.str() << '\n'
773      << "MPFR result: " << mpfr_result.str() << '\n';
774   __llvm_libc::fputil::testing::describeValue(
775       "Libc floating point result bits: ", libc_result, ss);
776   __llvm_libc::fputil::testing::describeValue(
777       "              MPFR rounded bits: ", mpfr_result.as<T>(), ss);
778   ss << "ULP error: " << std::to_string(mpfr_result.ulp(libc_result)) << '\n';
779   OS << ss.str();
780 }
781 
782 template void explain_ternary_operation_one_output_error<float>(
783     Operation, const TernaryInput<float> &, float, double, RoundingMode,
784     testutils::StreamWrapper &);
785 template void explain_ternary_operation_one_output_error<double>(
786     Operation, const TernaryInput<double> &, double, double, RoundingMode,
787     testutils::StreamWrapper &);
788 template void explain_ternary_operation_one_output_error<long double>(
789     Operation, const TernaryInput<long double> &, long double, double,
790     RoundingMode, testutils::StreamWrapper &);
791 
792 template <typename T>
793 bool compare_unary_operation_single_output(Operation op, T input, T libc_result,
794                                            double ulp_tolerance,
795                                            RoundingMode rounding) {
796   unsigned int precision = get_precision<T>(ulp_tolerance);
797   MPFRNumber mpfr_result;
798   mpfr_result = unary_operation(op, input, precision, rounding);
799   double ulp = mpfr_result.ulp(libc_result);
800   return (ulp <= ulp_tolerance);
801 }
802 
803 template bool compare_unary_operation_single_output<float>(Operation, float,
804                                                            float, double,
805                                                            RoundingMode);
806 template bool compare_unary_operation_single_output<double>(Operation, double,
807                                                             double, double,
808                                                             RoundingMode);
809 template bool compare_unary_operation_single_output<long double>(
810     Operation, long double, long double, double, RoundingMode);
811 
812 template <typename T>
813 bool compare_unary_operation_two_outputs(Operation op, T input,
814                                          const BinaryOutput<T> &libc_result,
815                                          double ulp_tolerance,
816                                          RoundingMode rounding) {
817   int mpfrIntResult;
818   unsigned int precision = get_precision<T>(ulp_tolerance);
819   MPFRNumber mpfr_result = unary_operation_two_outputs(op, input, mpfrIntResult,
820                                                        precision, rounding);
821   double ulp = mpfr_result.ulp(libc_result.f);
822 
823   if (mpfrIntResult != libc_result.i)
824     return false;
825 
826   return (ulp <= ulp_tolerance);
827 }
828 
829 template bool compare_unary_operation_two_outputs<float>(
830     Operation, float, const BinaryOutput<float> &, double, RoundingMode);
831 template bool compare_unary_operation_two_outputs<double>(
832     Operation, double, const BinaryOutput<double> &, double, RoundingMode);
833 template bool compare_unary_operation_two_outputs<long double>(
834     Operation, long double, const BinaryOutput<long double> &, double,
835     RoundingMode);
836 
837 template <typename T>
838 bool compare_binary_operation_two_outputs(Operation op,
839                                           const BinaryInput<T> &input,
840                                           const BinaryOutput<T> &libc_result,
841                                           double ulp_tolerance,
842                                           RoundingMode rounding) {
843   int mpfrIntResult;
844   unsigned int precision = get_precision<T>(ulp_tolerance);
845   MPFRNumber mpfr_result = binary_operation_two_outputs(
846       op, input.x, input.y, mpfrIntResult, precision, rounding);
847   double ulp = mpfr_result.ulp(libc_result.f);
848 
849   if (mpfrIntResult != libc_result.i) {
850     if (op == Operation::RemQuo) {
851       if ((0x7 & mpfrIntResult) != (0x7 & libc_result.i))
852         return false;
853     } else {
854       return false;
855     }
856   }
857 
858   return (ulp <= ulp_tolerance);
859 }
860 
861 template bool compare_binary_operation_two_outputs<float>(
862     Operation, const BinaryInput<float> &, const BinaryOutput<float> &, double,
863     RoundingMode);
864 template bool compare_binary_operation_two_outputs<double>(
865     Operation, const BinaryInput<double> &, const BinaryOutput<double> &,
866     double, RoundingMode);
867 template bool compare_binary_operation_two_outputs<long double>(
868     Operation, const BinaryInput<long double> &,
869     const BinaryOutput<long double> &, double, RoundingMode);
870 
871 template <typename T>
872 bool compare_binary_operation_one_output(Operation op,
873                                          const BinaryInput<T> &input,
874                                          T libc_result, double ulp_tolerance,
875                                          RoundingMode rounding) {
876   unsigned int precision = get_precision<T>(ulp_tolerance);
877   MPFRNumber mpfr_result =
878       binary_operation_one_output(op, input.x, input.y, precision, rounding);
879   double ulp = mpfr_result.ulp(libc_result);
880 
881   return (ulp <= ulp_tolerance);
882 }
883 
884 template bool compare_binary_operation_one_output<float>(
885     Operation, const BinaryInput<float> &, float, double, RoundingMode);
886 template bool compare_binary_operation_one_output<double>(
887     Operation, const BinaryInput<double> &, double, double, RoundingMode);
888 template bool compare_binary_operation_one_output<long double>(
889     Operation, const BinaryInput<long double> &, long double, double,
890     RoundingMode);
891 
892 template <typename T>
893 bool compare_ternary_operation_one_output(Operation op,
894                                           const TernaryInput<T> &input,
895                                           T libc_result, double ulp_tolerance,
896                                           RoundingMode rounding) {
897   unsigned int precision = get_precision<T>(ulp_tolerance);
898   MPFRNumber mpfr_result = ternary_operation_one_output(
899       op, input.x, input.y, input.z, precision, rounding);
900   double ulp = mpfr_result.ulp(libc_result);
901 
902   return (ulp <= ulp_tolerance);
903 }
904 
905 template bool compare_ternary_operation_one_output<float>(
906     Operation, const TernaryInput<float> &, float, double, RoundingMode);
907 template bool compare_ternary_operation_one_output<double>(
908     Operation, const TernaryInput<double> &, double, double, RoundingMode);
909 template bool compare_ternary_operation_one_output<long double>(
910     Operation, const TernaryInput<long double> &, long double, double,
911     RoundingMode);
912 
913 } // namespace internal
914 
915 template <typename T> bool round_to_long(T x, long &result) {
916   MPFRNumber mpfr(x);
917   return mpfr.round_to_long(result);
918 }
919 
920 template bool round_to_long<float>(float, long &);
921 template bool round_to_long<double>(double, long &);
922 template bool round_to_long<long double>(long double, long &);
923 
924 template <typename T> bool round_to_long(T x, RoundingMode mode, long &result) {
925   MPFRNumber mpfr(x);
926   return mpfr.round_to_long(get_mpfr_rounding_mode(mode), result);
927 }
928 
929 template bool round_to_long<float>(float, RoundingMode, long &);
930 template bool round_to_long<double>(double, RoundingMode, long &);
931 template bool round_to_long<long double>(long double, RoundingMode, long &);
932 
933 template <typename T> T round(T x, RoundingMode mode) {
934   MPFRNumber mpfr(x);
935   MPFRNumber result = mpfr.rint(get_mpfr_rounding_mode(mode));
936   return result.as<T>();
937 }
938 
939 template float round<float>(float, RoundingMode);
940 template double round<double>(double, RoundingMode);
941 template long double round<long double>(long double, RoundingMode);
942 
943 } // namespace mpfr
944 } // namespace testing
945 } // namespace __llvm_libc
946