diff --git a/libc/utils/MPFRWrapper/MPFRUtils.h b/libc/utils/MPFRWrapper/MPFRUtils.h
--- a/libc/utils/MPFRWrapper/MPFRUtils.h
+++ b/libc/utils/MPFRWrapper/MPFRUtils.h
@@ -71,6 +71,8 @@
   EndTernaryOperationsSingleOutput,
 };
 
+enum class RoundingMode : uint8_t { Upward, Downward, TowardZero, Nearest };
+
 template <typename T> struct BinaryInput {
   static_assert(
       __llvm_libc::cpp::IsFloatingPointType<T>::Value,
@@ -108,50 +110,56 @@
 
 template <typename T>
 bool compare_unary_operation_single_output(Operation op, T input, T libc_output,
-                                           double t);
+                                           double t, RoundingMode rounding);
 template <typename T>
 bool compare_unary_operation_two_outputs(Operation op, T input,
                                          const BinaryOutput<T> &libc_output,
-                                         double t);
+                                         double t, RoundingMode rounding);
 template <typename T>
 bool compare_binary_operation_two_outputs(Operation op,
                                           const BinaryInput<T> &input,
                                           const BinaryOutput<T> &libc_output,
-                                          double t);
+                                          double t, RoundingMode rounding);
 
 template <typename T>
 bool compare_binary_operation_one_output(Operation op,
                                          const BinaryInput<T> &input,
-                                         T libc_output, double t);
+                                         T libc_output, double t,
+                                         RoundingMode rounding);
 
 template <typename T>
 bool compare_ternary_operation_one_output(Operation op,
                                           const TernaryInput<T> &input,
-                                          T libc_output, double t);
+                                          T libc_output, double t,
+                                          RoundingMode rounding);
 
 template <typename T>
 void explain_unary_operation_single_output_error(Operation op, T input,
                                                  T match_value,
+                                                 RoundingMode rounding,
                                                  testutils::StreamWrapper &OS);
 template <typename T>
 void explain_unary_operation_two_outputs_error(
     Operation op, T input, const BinaryOutput<T> &match_value,
-    testutils::StreamWrapper &OS);
+    RoundingMode rounding, testutils::StreamWrapper &OS);
 template <typename T>
 void explain_binary_operation_two_outputs_error(
     Operation op, const BinaryInput<T> &input,
-    const BinaryOutput<T> &match_value, testutils::StreamWrapper &OS);
+    const BinaryOutput<T> &match_value, RoundingMode rounding,
+    testutils::StreamWrapper &OS);
 
 template <typename T>
 void explain_binary_operation_one_output_error(Operation op,
                                                const BinaryInput<T> &input,
                                                T match_value,
+                                               RoundingMode rounding,
                                                testutils::StreamWrapper &OS);
 
 template <typename T>
 void explain_ternary_operation_one_output_error(Operation op,
                                                 const TernaryInput<T> &input,
                                                 T match_value,
+                                                RoundingMode rounding,
                                                 testutils::StreamWrapper &OS);
 
 template <Operation op, typename InputType, typename OutputType>
@@ -159,14 +167,15 @@
   InputType input;
   OutputType match_value;
   double ulp_tolerance;
+  RoundingMode rounding;
 
 public:
-  MPFRMatcher(InputType testInput, double ulp_tolerance)
-      : input(testInput), ulp_tolerance(ulp_tolerance) {}
+  MPFRMatcher(InputType testInput, double ulp_tolerance, RoundingMode rounding)
+      : input(testInput), ulp_tolerance(ulp_tolerance), rounding(rounding) {}
 
   bool match(OutputType libcResult) {
     match_value = libcResult;
-    return match(input, match_value, ulp_tolerance);
+    return match(input, match_value);
   }
 
   // This method is marked with NOLINT because it the name `explainError`
@@ -176,59 +185,59 @@
   }
 
 private:
-  template <typename T> static bool match(T in, T out, double tolerance) {
-    return compare_unary_operation_single_output(op, in, out, tolerance);
+  template <typename T> bool match(T in, T out) {
+    return compare_unary_operation_single_output(op, in, out, ulp_tolerance,
+                                                 rounding);
   }
 
-  template <typename T>
-  static bool match(T in, const BinaryOutput<T> &out, double tolerance) {
-    return compare_unary_operation_two_outputs(op, in, out, tolerance);
+  template <typename T> bool match(T in, const BinaryOutput<T> &out) {
+    return compare_unary_operation_two_outputs(op, in, out, ulp_tolerance,
+                                               rounding);
   }
 
-  template <typename T>
-  static bool match(const BinaryInput<T> &in, T out, double tolerance) {
-    return compare_binary_operation_one_output(op, in, out, tolerance);
+  template <typename T> bool match(const BinaryInput<T> &in, T out) {
+    return compare_binary_operation_one_output(op, in, out, ulp_tolerance,
+                                               rounding);
   }
 
   template <typename T>
-  static bool match(BinaryInput<T> in, const BinaryOutput<T> &out,
-                    double tolerance) {
-    return compare_binary_operation_two_outputs(op, in, out, tolerance);
+  bool match(BinaryInput<T> in, const BinaryOutput<T> &out) {
+    return compare_binary_operation_two_outputs(op, in, out, ulp_tolerance,
+                                                rounding);
   }
 
-  template <typename T>
-  static bool match(const TernaryInput<T> &in, T out, double tolerance) {
-    return compare_ternary_operation_one_output(op, in, out, tolerance);
+  template <typename T> bool match(const TernaryInput<T> &in, T out) {
+    return compare_ternary_operation_one_output(op, in, out, ulp_tolerance,
+                                                rounding);
   }
 
   template <typename T>
-  static void explain_error(T in, T out, testutils::StreamWrapper &OS) {
-    explain_unary_operation_single_output_error(op, in, out, OS);
+  void explain_error(T in, T out, testutils::StreamWrapper &OS) {
+    explain_unary_operation_single_output_error(op, in, out, rounding, OS);
   }
 
   template <typename T>
-  static void explain_error(T in, const BinaryOutput<T> &out,
-                            testutils::StreamWrapper &OS) {
-    explain_unary_operation_two_outputs_error(op, in, out, OS);
+  void explain_error(T in, const BinaryOutput<T> &out,
+                     testutils::StreamWrapper &OS) {
+    explain_unary_operation_two_outputs_error(op, in, out, rounding, OS);
   }
 
   template <typename T>
-  static void explain_error(const BinaryInput<T> &in,
-                            const BinaryOutput<T> &out,
-                            testutils::StreamWrapper &OS) {
-    explain_binary_operation_two_outputs_error(op, in, out, OS);
+  void explain_error(const BinaryInput<T> &in, const BinaryOutput<T> &out,
+                     testutils::StreamWrapper &OS) {
+    explain_binary_operation_two_outputs_error(op, in, out, rounding, OS);
   }
 
   template <typename T>
-  static void explain_error(const BinaryInput<T> &in, T out,
-                            testutils::StreamWrapper &OS) {
-    explain_binary_operation_one_output_error(op, in, out, OS);
+  void explain_error(const BinaryInput<T> &in, T out,
+                     testutils::StreamWrapper &OS) {
+    explain_binary_operation_one_output_error(op, in, out, rounding, OS);
   }
 
   template <typename T>
-  static void explain_error(const TernaryInput<T> &in, T out,
-                            testutils::StreamWrapper &OS) {
-    explain_ternary_operation_one_output_error(op, in, out, OS);
+  void explain_error(const TernaryInput<T> &in, T out,
+                     testutils::StreamWrapper &OS) {
+    explain_ternary_operation_one_output_error(op, in, out, rounding, OS);
   }
 };
 
@@ -264,12 +273,12 @@
 __attribute__((no_sanitize("address")))
 cpp::EnableIfType<is_valid_operation<op, InputType, OutputType>(),
                   internal::MPFRMatcher<op, InputType, OutputType>>
-get_mpfr_matcher(InputType input, OutputType output_unused, double t) {
-  return internal::MPFRMatcher<op, InputType, OutputType>(input, t);
+get_mpfr_matcher(InputType input, OutputType output_unused, double tolerance,
+                 RoundingMode rounding) {
+  return internal::MPFRMatcher<op, InputType, OutputType>(input, tolerance,
+                                                          rounding);
 }
 
-enum class RoundingMode : uint8_t { Upward, Downward, TowardZero, Nearest };
-
 template <typename T> T round(T x, RoundingMode mode);
 
 template <typename T> bool round_to_long(T x, long &result);
@@ -279,12 +288,40 @@
 } // namespace testing
 } // namespace __llvm_libc
 
-#define EXPECT_MPFR_MATCH(op, input, match_value, tolerance)                   \
+#define MISSING_ARGS(...) static_assert(true, "Missing arguments");
+
+#define GET_MPFR_MACRO(_1, _2, _3, _4, _5, NAME, ...) NAME
+
+#define EXPECT_MPFR_MATCH_DEFAULT(op, input, match_value, tolerance)           \
+  EXPECT_THAT(match_value,                                                     \
+              __llvm_libc::testing::mpfr::get_mpfr_matcher<op>(                \
+                  input, match_value, tolerance,                               \
+                  __llvm_libc::testing::mpfr::RoundingMode::Nearest))
+
+#define EXPECT_MPFR_MATCH_ROUNDING(op, input, match_value, tolerance,          \
+                                   rounding)                                   \
   EXPECT_THAT(match_value, __llvm_libc::testing::mpfr::get_mpfr_matcher<op>(   \
-                               input, match_value, tolerance))
+                               input, match_value, tolerance, rounding))
 
-#define ASSERT_MPFR_MATCH(op, input, match_value, tolerance)                   \
+#define EXPECT_MPFR_MATCH(...)                                                 \
+  GET_MPFR_MACRO(__VA_ARGS__, EXPECT_MPFR_MATCH_ROUNDING,                      \
+                 EXPECT_MPFR_MATCH_DEFAULT, MISSING_ARGS)                      \
+  (__VA_ARGS__)
+
+#define ASSERT_MPFR_MATCH_DEFAULT(op, input, match_value, tolerance)           \
+  ASSERT_THAT(match_value,                                                     \
+              __llvm_libc::testing::mpfr::get_mpfr_matcher<op>(                \
+                  input, match_value, tolerance,                               \
+                  __llvm_libc::testing::mpfr::RoundingMode::Nearest))
+
+#define ASSERT_MPFR_MATCH_ROUNDING(op, input, match_value, tolerance,          \
+                                   rounding)                                   \
   ASSERT_THAT(match_value, __llvm_libc::testing::mpfr::get_mpfr_matcher<op>(   \
-                               input, match_value, tolerance))
+                               input, match_value, tolerance, rounding))
+
+#define ASSERT_MPFR_MATCH(...)                                                 \
+  GET_MPFR_MACRO(__VA_ARGS__, ASSERT_MPFR_MATCH_ROUNDING,                      \
+                 ASSERT_MPFR_MATCH_DEFAULT, MISSING_ARGS)                      \
+  (__VA_ARGS__)
 
 #endif // LLVM_LIBC_UTILS_TESTUTILS_MPFRUTILS_H
diff --git a/libc/utils/MPFRWrapper/MPFRUtils.cpp b/libc/utils/MPFRWrapper/MPFRUtils.cpp
--- a/libc/utils/MPFRWrapper/MPFRUtils.cpp
+++ b/libc/utils/MPFRWrapper/MPFRUtils.cpp
@@ -55,67 +55,113 @@
 };
 #endif
 
+template <typename T> struct ExtraPrecision;
+
+template <> struct ExtraPrecision<float> {
+  static constexpr unsigned int VALUE = 128;
+};
+
+template <> struct ExtraPrecision<double> {
+  static constexpr unsigned int VALUE = 256;
+};
+
+template <> struct ExtraPrecision<long double> {
+  static constexpr unsigned int VALUE = 256;
+};
+
+static inline mpfr_rnd_t get_mpfr_rounding_mode(RoundingMode mode) {
+  switch (mode) {
+  case RoundingMode::Upward:
+    return MPFR_RNDU;
+    break;
+  case RoundingMode::Downward:
+    return MPFR_RNDD;
+    break;
+  case RoundingMode::TowardZero:
+    return MPFR_RNDZ;
+    break;
+  case RoundingMode::Nearest:
+    return MPFR_RNDN;
+    break;
+  }
+}
+
 class MPFRNumber {
   // A precision value which allows sufficiently large additional
   // precision even compared to quad-precision floating point values.
   unsigned int mpfr_precision;
+  mpfr_rnd_t mpfr_rounding;
 
   mpfr_t value;
 
 public:
-  MPFRNumber() : mpfr_precision(256) { mpfr_init2(value, mpfr_precision); }
+  MPFRNumber() : mpfr_precision(256), mpfr_rounding(MPFR_RNDN) {
+    mpfr_init2(value, mpfr_precision);
+  }
 
   // We use explicit EnableIf specializations to disallow implicit
   // conversions. Implicit conversions can potentially lead to loss of
   // precision.
   template <typename XType,
             cpp::EnableIfType<cpp::IsSame<float, XType>::Value, int> = 0>
-  explicit MPFRNumber(XType x, int precision = 128)
-      : mpfr_precision(precision) {
+  explicit MPFRNumber(XType x, int precision = ExtraPrecision<float>::VALUE,
+                      RoundingMode rounding = RoundingMode::Nearest)
+      : mpfr_precision(precision),
+        mpfr_rounding(get_mpfr_rounding_mode(rounding)) {
     mpfr_init2(value, mpfr_precision);
-    mpfr_set_flt(value, x, MPFR_RNDN);
+    mpfr_set_flt(value, x, mpfr_rounding);
   }
 
   template <typename XType,
             cpp::EnableIfType<cpp::IsSame<double, XType>::Value, int> = 0>
-  explicit MPFRNumber(XType x, int precision = 128)
-      : mpfr_precision(precision) {
+  explicit MPFRNumber(XType x, int precision = ExtraPrecision<double>::VALUE,
+                      RoundingMode rounding = RoundingMode::Nearest)
+      : mpfr_precision(precision),
+        mpfr_rounding(get_mpfr_rounding_mode(rounding)) {
     mpfr_init2(value, mpfr_precision);
-    mpfr_set_d(value, x, MPFR_RNDN);
+    mpfr_set_d(value, x, mpfr_rounding);
   }
 
   template <typename XType,
             cpp::EnableIfType<cpp::IsSame<long double, XType>::Value, int> = 0>
-  explicit MPFRNumber(XType x, int precision = 128)
-      : mpfr_precision(precision) {
+  explicit MPFRNumber(XType x,
+                      int precision = ExtraPrecision<long double>::VALUE,
+                      RoundingMode rounding = RoundingMode::Nearest)
+      : mpfr_precision(precision),
+        mpfr_rounding(get_mpfr_rounding_mode(rounding)) {
     mpfr_init2(value, mpfr_precision);
-    mpfr_set_ld(value, x, MPFR_RNDN);
+    mpfr_set_ld(value, x, mpfr_rounding);
   }
 
   template <typename XType,
             cpp::EnableIfType<cpp::IsIntegral<XType>::Value, int> = 0>
-  explicit MPFRNumber(XType x, int precision = 128)
-      : mpfr_precision(precision) {
+  explicit MPFRNumber(XType x, int precision = ExtraPrecision<float>::VALUE,
+                      RoundingMode rounding = RoundingMode::Nearest)
+      : mpfr_precision(precision),
+        mpfr_rounding(get_mpfr_rounding_mode(rounding)) {
     mpfr_init2(value, mpfr_precision);
-    mpfr_set_sj(value, x, MPFR_RNDN);
+    mpfr_set_sj(value, x, mpfr_rounding);
   }
 
-  MPFRNumber(const MPFRNumber &other) : mpfr_precision(other.mpfr_precision) {
+  MPFRNumber(const MPFRNumber &other)
+      : mpfr_precision(other.mpfr_precision),
+        mpfr_rounding(other.mpfr_rounding) {
     mpfr_init2(value, mpfr_precision);
-    mpfr_set(value, other.value, MPFR_RNDN);
+    mpfr_set(value, other.value, mpfr_rounding);
   }
 
   ~MPFRNumber() { mpfr_clear(value); }
 
   MPFRNumber &operator=(const MPFRNumber &rhs) {
     mpfr_precision = rhs.mpfr_precision;
-    mpfr_set(value, rhs.value, MPFR_RNDN);
+    mpfr_rounding = rhs.mpfr_rounding;
+    mpfr_set(value, rhs.value, mpfr_rounding);
     return *this;
   }
 
   MPFRNumber abs() const {
     MPFRNumber result;
-    mpfr_abs(result.value, value, MPFR_RNDN);
+    mpfr_abs(result.value, value, mpfr_rounding);
     return result;
   }
 
@@ -127,25 +173,25 @@
 
   MPFRNumber cos() const {
     MPFRNumber result;
-    mpfr_cos(result.value, value, MPFR_RNDN);
+    mpfr_cos(result.value, value, mpfr_rounding);
     return result;
   }
 
   MPFRNumber exp() const {
     MPFRNumber result;
-    mpfr_exp(result.value, value, MPFR_RNDN);
+    mpfr_exp(result.value, value, mpfr_rounding);
     return result;
   }
 
   MPFRNumber exp2() const {
     MPFRNumber result;
-    mpfr_exp2(result.value, value, MPFR_RNDN);
+    mpfr_exp2(result.value, value, mpfr_rounding);
     return result;
   }
 
   MPFRNumber expm1() const {
     MPFRNumber result;
-    mpfr_expm1(result.value, value, MPFR_RNDN);
+    mpfr_expm1(result.value, value, mpfr_rounding);
     return result;
   }
 
@@ -158,27 +204,27 @@
   MPFRNumber frexp(int &exp) {
     MPFRNumber result;
     mpfr_exp_t resultExp;
-    mpfr_frexp(&resultExp, result.value, value, MPFR_RNDN);
+    mpfr_frexp(&resultExp, result.value, value, mpfr_rounding);
     exp = resultExp;
     return result;
   }
 
   MPFRNumber hypot(const MPFRNumber &b) {
     MPFRNumber result;
-    mpfr_hypot(result.value, value, b.value, MPFR_RNDN);
+    mpfr_hypot(result.value, value, b.value, mpfr_rounding);
     return result;
   }
 
   MPFRNumber log() const {
     MPFRNumber result;
-    mpfr_log(result.value, value, MPFR_RNDN);
+    mpfr_log(result.value, value, mpfr_rounding);
     return result;
   }
 
   MPFRNumber remquo(const MPFRNumber &divisor, int &quotient) {
     MPFRNumber remainder;
     long q;
-    mpfr_remquo(remainder.value, &q, value, divisor.value, MPFR_RNDN);
+    mpfr_remquo(remainder.value, &q, value, divisor.value, mpfr_rounding);
     quotient = q;
     return remainder;
   }
@@ -240,19 +286,19 @@
 
   MPFRNumber sin() const {
     MPFRNumber result;
-    mpfr_sin(result.value, value, MPFR_RNDN);
+    mpfr_sin(result.value, value, mpfr_rounding);
     return result;
   }
 
   MPFRNumber sqrt() const {
     MPFRNumber result;
-    mpfr_sqrt(result.value, value, MPFR_RNDN);
+    mpfr_sqrt(result.value, value, mpfr_rounding);
     return result;
   }
 
   MPFRNumber tan() const {
     MPFRNumber result;
-    mpfr_tan(result.value, value, MPFR_RNDN);
+    mpfr_tan(result.value, value, mpfr_rounding);
     return result;
   }
 
@@ -264,7 +310,7 @@
 
   MPFRNumber fma(const MPFRNumber &b, const MPFRNumber &c) {
     MPFRNumber result(*this);
-    mpfr_fma(result.value, value, b.value, c.value, MPFR_RNDN);
+    mpfr_fma(result.value, value, b.value, c.value, mpfr_rounding);
     return result;
   }
 
@@ -282,10 +328,14 @@
   // These functions are useful for debugging.
   template <typename T> T as() const;
 
-  template <> float as<float>() const { return mpfr_get_flt(value, MPFR_RNDN); }
-  template <> double as<double>() const { return mpfr_get_d(value, MPFR_RNDN); }
+  template <> float as<float>() const {
+    return mpfr_get_flt(value, mpfr_rounding);
+  }
+  template <> double as<double>() const {
+    return mpfr_get_d(value, mpfr_rounding);
+  }
   template <> long double as<long double>() const {
-    return mpfr_get_ld(value, MPFR_RNDN);
+    return mpfr_get_ld(value, mpfr_rounding);
   }
 
   void dump(const char *msg) const { mpfr_printf("%s%.128Rf\n", msg, value); }
@@ -378,8 +428,8 @@
 
 template <typename InputType>
 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
-unary_operation(Operation op, InputType input) {
-  MPFRNumber mpfrInput(input);
+unary_operation(Operation op, InputType input, RoundingMode rounding) {
+  MPFRNumber mpfrInput(input, ExtraPrecision<InputType>::VALUE, rounding);
   switch (op) {
   case Operation::Abs:
     return mpfrInput.abs();
@@ -420,8 +470,9 @@
 
 template <typename InputType>
 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
-unary_operation_two_outputs(Operation op, InputType input, int &output) {
-  MPFRNumber mpfrInput(input);
+unary_operation_two_outputs(Operation op, InputType input, int &output,
+                            RoundingMode rounding) {
+  MPFRNumber mpfrInput(input, ExtraPrecision<InputType>::VALUE, rounding);
   switch (op) {
   case Operation::Frexp:
     return mpfrInput.frexp(output);
@@ -432,8 +483,10 @@
 
 template <typename InputType>
 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
-binary_operation_one_output(Operation op, InputType x, InputType y) {
-  MPFRNumber inputX(x), inputY(y);
+binary_operation_one_output(Operation op, InputType x, InputType y,
+                            RoundingMode rounding) {
+  MPFRNumber inputX(x, ExtraPrecision<InputType>::VALUE, rounding);
+  MPFRNumber inputY(y, ExtraPrecision<InputType>::VALUE, rounding);
   switch (op) {
   case Operation::Hypot:
     return inputX.hypot(inputY);
@@ -445,8 +498,9 @@
 template <typename InputType>
 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
 binary_operation_two_outputs(Operation op, InputType x, InputType y,
-                             int &output) {
-  MPFRNumber inputX(x), inputY(y);
+                             int &output, RoundingMode rounding) {
+  MPFRNumber inputX(x, ExtraPrecision<InputType>::VALUE, rounding);
+  MPFRNumber inputY(y, ExtraPrecision<InputType>::VALUE, rounding);
   switch (op) {
   case Operation::RemQuo:
     return inputX.remquo(inputY, output);
@@ -458,12 +512,14 @@
 template <typename InputType>
 cpp::EnableIfType<cpp::IsFloatingPointType<InputType>::Value, MPFRNumber>
 ternary_operation_one_output(Operation op, InputType x, InputType y,
-                             InputType z) {
+                             InputType z, RoundingMode rounding) {
   // For FMA function, we just need to compare with the mpfr_fma with the same
   // precision as InputType.  Using higher precision as the intermediate results
   // to compare might incorrectly fail due to double-rounding errors.
   constexpr unsigned int prec = Precision<InputType>::VALUE;
-  MPFRNumber inputX(x, prec), inputY(y, prec), inputZ(z, prec);
+  MPFRNumber inputX(x, prec, rounding);
+  MPFRNumber inputY(y, prec, rounding);
+  MPFRNumber inputZ(z, prec, rounding);
   switch (op) {
   case Operation::Fma:
     return inputX.fma(inputY, inputZ);
@@ -475,13 +531,11 @@
 template <typename T>
 void explain_unary_operation_single_output_error(Operation op, T input,
                                                  T matchValue,
+                                                 RoundingMode rounding,
                                                  testutils::StreamWrapper &OS) {
   MPFRNumber mpfrInput(input);
-  MPFRNumber mpfr_result = unary_operation(op, input);
+  MPFRNumber mpfr_result = unary_operation(op, input, rounding);
   MPFRNumber mpfrMatchValue(matchValue);
-  FPBits<T> inputBits(input);
-  FPBits<T> matchBits(matchValue);
-  FPBits<T> mpfr_resultBits(mpfr_result.as<T>());
   OS << "Match value not within tolerance value of MPFR result:\n"
      << "  Input decimal: " << mpfrInput.str() << '\n';
   __llvm_libc::fputil::testing::describeValue("     Input bits: ", input, OS);
@@ -496,23 +550,22 @@
      << '\n';
 }
 
-template void
-explain_unary_operation_single_output_error<float>(Operation op, float, float,
-                                                   testutils::StreamWrapper &);
+template void explain_unary_operation_single_output_error<float>(
+    Operation op, float, float, RoundingMode, testutils::StreamWrapper &);
 template void explain_unary_operation_single_output_error<double>(
-    Operation op, double, double, testutils::StreamWrapper &);
+    Operation op, double, double, RoundingMode, testutils::StreamWrapper &);
 template void explain_unary_operation_single_output_error<long double>(
-    Operation op, long double, long double, testutils::StreamWrapper &);
+    Operation op, long double, long double, RoundingMode,
+    testutils::StreamWrapper &);
 
 template <typename T>
 void explain_unary_operation_two_outputs_error(
     Operation op, T input, const BinaryOutput<T> &libc_result,
-    testutils::StreamWrapper &OS) {
+    RoundingMode rounding, testutils::StreamWrapper &OS) {
   MPFRNumber mpfrInput(input);
-  FPBits<T> inputBits(input);
   int mpfrIntResult;
   MPFRNumber mpfr_result =
-      unary_operation_two_outputs(op, input, mpfrIntResult);
+      unary_operation_two_outputs(op, input, mpfrIntResult, rounding);
 
   if (mpfrIntResult != libc_result.i) {
     OS << "MPFR integral result: " << mpfrIntResult << '\n'
@@ -541,26 +594,25 @@
 }
 
 template void explain_unary_operation_two_outputs_error<float>(
-    Operation, float, const BinaryOutput<float> &, testutils::StreamWrapper &);
-template void
-explain_unary_operation_two_outputs_error<double>(Operation, double,
-                                                  const BinaryOutput<double> &,
-                                                  testutils::StreamWrapper &);
+    Operation, float, const BinaryOutput<float> &, RoundingMode,
+    testutils::StreamWrapper &);
+template void explain_unary_operation_two_outputs_error<double>(
+    Operation, double, const BinaryOutput<double> &, RoundingMode,
+    testutils::StreamWrapper &);
 template void explain_unary_operation_two_outputs_error<long double>(
-    Operation, long double, const BinaryOutput<long double> &,
+    Operation, long double, const BinaryOutput<long double> &, RoundingMode,
     testutils::StreamWrapper &);
 
 template <typename T>
 void explain_binary_operation_two_outputs_error(
     Operation op, const BinaryInput<T> &input,
-    const BinaryOutput<T> &libc_result, testutils::StreamWrapper &OS) {
+    const BinaryOutput<T> &libc_result, RoundingMode rounding,
+    testutils::StreamWrapper &OS) {
   MPFRNumber mpfrX(input.x);
   MPFRNumber mpfrY(input.y);
-  FPBits<T> xbits(input.x);
-  FPBits<T> ybits(input.y);
   int mpfrIntResult;
-  MPFRNumber mpfr_result =
-      binary_operation_two_outputs(op, input.x, input.y, mpfrIntResult);
+  MPFRNumber mpfr_result = binary_operation_two_outputs(
+      op, input.x, input.y, mpfrIntResult, rounding);
   MPFRNumber mpfrMatchValue(libc_result.f);
 
   OS << "Input decimal: x: " << mpfrX.str() << " y: " << mpfrY.str() << '\n'
@@ -577,24 +629,27 @@
 
 template void explain_binary_operation_two_outputs_error<float>(
     Operation, const BinaryInput<float> &, const BinaryOutput<float> &,
-    testutils::StreamWrapper &);
+    RoundingMode, testutils::StreamWrapper &);
 template void explain_binary_operation_two_outputs_error<double>(
     Operation, const BinaryInput<double> &, const BinaryOutput<double> &,
-    testutils::StreamWrapper &);
+    RoundingMode, testutils::StreamWrapper &);
 template void explain_binary_operation_two_outputs_error<long double>(
     Operation, const BinaryInput<long double> &,
-    const BinaryOutput<long double> &, testutils::StreamWrapper &);
+    const BinaryOutput<long double> &, RoundingMode,
+    testutils::StreamWrapper &);
 
 template <typename T>
 void explain_binary_operation_one_output_error(Operation op,
                                                const BinaryInput<T> &input,
                                                T libc_result,
+                                               RoundingMode rounding,
                                                testutils::StreamWrapper &OS) {
   MPFRNumber mpfrX(input.x);
   MPFRNumber mpfrY(input.y);
   FPBits<T> xbits(input.x);
   FPBits<T> ybits(input.y);
-  MPFRNumber mpfr_result = binary_operation_one_output(op, input.x, input.y);
+  MPFRNumber mpfr_result =
+      binary_operation_one_output(op, input.x, input.y, rounding);
   MPFRNumber mpfrMatchValue(libc_result);
 
   OS << "Input decimal: x: " << mpfrX.str() << " y: " << mpfrY.str() << '\n';
@@ -613,17 +668,20 @@
 }
 
 template void explain_binary_operation_one_output_error<float>(
-    Operation, const BinaryInput<float> &, float, testutils::StreamWrapper &);
+    Operation, const BinaryInput<float> &, float, RoundingMode,
+    testutils::StreamWrapper &);
 template void explain_binary_operation_one_output_error<double>(
-    Operation, const BinaryInput<double> &, double, testutils::StreamWrapper &);
+    Operation, const BinaryInput<double> &, double, RoundingMode,
+    testutils::StreamWrapper &);
 template void explain_binary_operation_one_output_error<long double>(
-    Operation, const BinaryInput<long double> &, long double,
+    Operation, const BinaryInput<long double> &, long double, RoundingMode,
     testutils::StreamWrapper &);
 
 template <typename T>
 void explain_ternary_operation_one_output_error(Operation op,
                                                 const TernaryInput<T> &input,
                                                 T libc_result,
+                                                RoundingMode rounding,
                                                 testutils::StreamWrapper &OS) {
   MPFRNumber mpfrX(input.x, Precision<T>::VALUE);
   MPFRNumber mpfrY(input.y, Precision<T>::VALUE);
@@ -632,7 +690,7 @@
   FPBits<T> ybits(input.y);
   FPBits<T> zbits(input.z);
   MPFRNumber mpfr_result =
-      ternary_operation_one_output(op, input.x, input.y, input.z);
+      ternary_operation_one_output(op, input.x, input.y, input.z, rounding);
   MPFRNumber mpfrMatchValue(libc_result);
 
   OS << "Input decimal: x: " << mpfrX.str() << " y: " << mpfrY.str()
@@ -654,20 +712,22 @@
 }
 
 template void explain_ternary_operation_one_output_error<float>(
-    Operation, const TernaryInput<float> &, float, testutils::StreamWrapper &);
+    Operation, const TernaryInput<float> &, float, RoundingMode,
+    testutils::StreamWrapper &);
 template void explain_ternary_operation_one_output_error<double>(
-    Operation, const TernaryInput<double> &, double,
+    Operation, const TernaryInput<double> &, double, RoundingMode,
     testutils::StreamWrapper &);
 template void explain_ternary_operation_one_output_error<long double>(
-    Operation, const TernaryInput<long double> &, long double,
+    Operation, const TernaryInput<long double> &, long double, RoundingMode,
     testutils::StreamWrapper &);
 
 template <typename T>
 bool compare_unary_operation_single_output(Operation op, T input, T libc_result,
-                                           double ulp_error) {
+                                           double ulp_error,
+                                           RoundingMode rounding) {
   // If the ulp error is exactly 0.5 (i.e a tie), we would check that the result
   // is rounded to the nearest even.
-  MPFRNumber mpfr_result = unary_operation(op, input);
+  MPFRNumber mpfr_result = unary_operation(op, input, rounding);
   double ulp = mpfr_result.ulp(libc_result);
   bool bits_are_even = ((FPBits<T>(libc_result).uintval() & 1) == 0);
   return (ulp < ulp_error) ||
@@ -675,21 +735,22 @@
 }
 
 template bool compare_unary_operation_single_output<float>(Operation, float,
-                                                           float, double);
+                                                           float, double,
+                                                           RoundingMode);
 template bool compare_unary_operation_single_output<double>(Operation, double,
-                                                            double, double);
-template bool compare_unary_operation_single_output<long double>(Operation,
-                                                                 long double,
-                                                                 long double,
-                                                                 double);
+                                                            double, double,
+                                                            RoundingMode);
+template bool compare_unary_operation_single_output<long double>(
+    Operation, long double, long double, double, RoundingMode);
 
 template <typename T>
 bool compare_unary_operation_two_outputs(Operation op, T input,
                                          const BinaryOutput<T> &libc_result,
-                                         double ulp_error) {
+                                         double ulp_error,
+                                         RoundingMode rounding) {
   int mpfrIntResult;
   MPFRNumber mpfr_result =
-      unary_operation_two_outputs(op, input, mpfrIntResult);
+      unary_operation_two_outputs(op, input, mpfrIntResult, rounding);
   double ulp = mpfr_result.ulp(libc_result.f);
 
   if (mpfrIntResult != libc_result.i)
@@ -700,22 +761,23 @@
          ((ulp == ulp_error) && ((ulp != 0.5) || bits_are_even));
 }
 
-template bool
-compare_unary_operation_two_outputs<float>(Operation, float,
-                                           const BinaryOutput<float> &, double);
+template bool compare_unary_operation_two_outputs<float>(
+    Operation, float, const BinaryOutput<float> &, double, RoundingMode);
 template bool compare_unary_operation_two_outputs<double>(
-    Operation, double, const BinaryOutput<double> &, double);
+    Operation, double, const BinaryOutput<double> &, double, RoundingMode);
 template bool compare_unary_operation_two_outputs<long double>(
-    Operation, long double, const BinaryOutput<long double> &, double);
+    Operation, long double, const BinaryOutput<long double> &, double,
+    RoundingMode);
 
 template <typename T>
 bool compare_binary_operation_two_outputs(Operation op,
                                           const BinaryInput<T> &input,
                                           const BinaryOutput<T> &libc_result,
-                                          double ulp_error) {
+                                          double ulp_error,
+                                          RoundingMode rounding) {
   int mpfrIntResult;
-  MPFRNumber mpfr_result =
-      binary_operation_two_outputs(op, input.x, input.y, mpfrIntResult);
+  MPFRNumber mpfr_result = binary_operation_two_outputs(
+      op, input.x, input.y, mpfrIntResult, rounding);
   double ulp = mpfr_result.ulp(libc_result.f);
 
   if (mpfrIntResult != libc_result.i) {
@@ -733,19 +795,22 @@
 }
 
 template bool compare_binary_operation_two_outputs<float>(
-    Operation, const BinaryInput<float> &, const BinaryOutput<float> &, double);
+    Operation, const BinaryInput<float> &, const BinaryOutput<float> &, double,
+    RoundingMode);
 template bool compare_binary_operation_two_outputs<double>(
     Operation, const BinaryInput<double> &, const BinaryOutput<double> &,
-    double);
+    double, RoundingMode);
 template bool compare_binary_operation_two_outputs<long double>(
     Operation, const BinaryInput<long double> &,
-    const BinaryOutput<long double> &, double);
+    const BinaryOutput<long double> &, double, RoundingMode);
 
 template <typename T>
 bool compare_binary_operation_one_output(Operation op,
                                          const BinaryInput<T> &input,
-                                         T libc_result, double ulp_error) {
-  MPFRNumber mpfr_result = binary_operation_one_output(op, input.x, input.y);
+                                         T libc_result, double ulp_error,
+                                         RoundingMode rounding) {
+  MPFRNumber mpfr_result =
+      binary_operation_one_output(op, input.x, input.y, rounding);
   double ulp = mpfr_result.ulp(libc_result);
 
   bool bits_are_even = ((FPBits<T>(libc_result).uintval() & 1) == 0);
@@ -754,18 +819,20 @@
 }
 
 template bool compare_binary_operation_one_output<float>(
-    Operation, const BinaryInput<float> &, float, double);
+    Operation, const BinaryInput<float> &, float, double, RoundingMode);
 template bool compare_binary_operation_one_output<double>(
-    Operation, const BinaryInput<double> &, double, double);
+    Operation, const BinaryInput<double> &, double, double, RoundingMode);
 template bool compare_binary_operation_one_output<long double>(
-    Operation, const BinaryInput<long double> &, long double, double);
+    Operation, const BinaryInput<long double> &, long double, double,
+    RoundingMode);
 
 template <typename T>
 bool compare_ternary_operation_one_output(Operation op,
                                           const TernaryInput<T> &input,
-                                          T libc_result, double ulp_error) {
+                                          T libc_result, double ulp_error,
+                                          RoundingMode rounding) {
   MPFRNumber mpfr_result =
-      ternary_operation_one_output(op, input.x, input.y, input.z);
+      ternary_operation_one_output(op, input.x, input.y, input.z, rounding);
   double ulp = mpfr_result.ulp(libc_result);
 
   bool bits_are_even = ((FPBits<T>(libc_result).uintval() & 1) == 0);
@@ -774,28 +841,12 @@
 }
 
 template bool compare_ternary_operation_one_output<float>(
-    Operation, const TernaryInput<float> &, float, double);
+    Operation, const TernaryInput<float> &, float, double, RoundingMode);
 template bool compare_ternary_operation_one_output<double>(
-    Operation, const TernaryInput<double> &, double, double);
+    Operation, const TernaryInput<double> &, double, double, RoundingMode);
 template bool compare_ternary_operation_one_output<long double>(
-    Operation, const TernaryInput<long double> &, long double, double);
-
-static mpfr_rnd_t get_mpfr_rounding_mode(RoundingMode mode) {
-  switch (mode) {
-  case RoundingMode::Upward:
-    return MPFR_RNDU;
-    break;
-  case RoundingMode::Downward:
-    return MPFR_RNDD;
-    break;
-  case RoundingMode::TowardZero:
-    return MPFR_RNDZ;
-    break;
-  case RoundingMode::Nearest:
-    return MPFR_RNDN;
-    break;
-  }
-}
+    Operation, const TernaryInput<long double> &, long double, double,
+    RoundingMode);
 
 } // namespace internal
 
@@ -810,7 +861,7 @@
 
 template <typename T> bool round_to_long(T x, RoundingMode mode, long &result) {
   MPFRNumber mpfr(x);
-  return mpfr.roung_to_long(internal::get_mpfr_rounding_mode(mode), result);
+  return mpfr.roung_to_long(get_mpfr_rounding_mode(mode), result);
 }
 
 template bool round_to_long<float>(float, RoundingMode, long &);
@@ -819,7 +870,7 @@
 
 template <typename T> T round(T x, RoundingMode mode) {
   MPFRNumber mpfr(x);
-  MPFRNumber result = mpfr.rint(internal::get_mpfr_rounding_mode(mode));
+  MPFRNumber result = mpfr.rint(get_mpfr_rounding_mode(mode));
   return result.as<T>();
 }