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