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