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