How to set rounding mode in Boost Multiprecision when using MPFR

Question

I am trying to figure out how I can format mpfr_float numbers using a rounding mode in Boost Multiprecision. In the below example, I expect 1.55 to round to either 1.5 or 1.6 depending on which rounding mode is used, but instead for all cases it outputs 1.5. How can I achieve this simple functionality in Boost with MPFR?

#include <iostream>
#include <boost/multiprecision/mpfr.hpp>

void setRoundingMode(boost::multiprecision::mpfr_float m, mpfr_rnd_t r)
{
    mpfr_t tmp;
    mpfr_init(tmp);
    mpfr_set(tmp, m.backend().data(), r);
    mpfr_clear(tmp);
}

int main()
{
    using namespace boost::multiprecision;
    using std::cout;
    using std::endl;
    using std::setprecision;

    mpfr_float::default_precision(50);
    mpfr_float a("1.55");

    setRoundingMode(a, MPFR_RNDN); /* round to nearest, with ties to even */
    cout << setprecision(2) << a << endl;

    setRoundingMode(a, MPFR_RNDZ); /* round toward zero */
    cout << setprecision(2) << a << endl;

    setRoundingMode(a, MPFR_RNDU); /* round toward +Inf */
    cout << setprecision(2) << a << endl;

    setRoundingMode(a, MPFR_RNDD); /* round toward -Inf */
    cout << setprecision(2) << a << endl;

    setRoundingMode(a, MPFR_RNDA); /* round away from zero */
    cout << setprecision(2) << a << endl;

    setRoundingMode(a, MPFR_RNDF); /* faithful rounding */
    cout << setprecision(2) << a << endl;

    setRoundingMode(a, MPFR_RNDNA); /* round to nearest, with ties away from zero (mpfr_round) */
    cout << setprecision(2) << a << endl;

    return 0;
}

Answer 1

The docs say:

All conversions are performed by the underlying MPFR library

However the doc for mpfr float backend states:

Things you should know when using this type:

• A default constructed mpfr_float_backend is set to zero (Note that this is not the default MPFR behavior).
• All operations use round to nearest.

(emphasis mine)

I've found that that MPFR doesn't have a global default rounding override, so it's only specific operations (assignment and precision change) that take a mpfr_rnd_t.

You realized this, hence your:

void setRoundingMode(boost::multiprecision::mpfr_float m, mpfr_rnd_t r)
{
    mpfr_t tmp;
    mpfr_init(tmp);
    mpfr_set(tmp, m.backend().data(), r);
    mpfr_clear(tmp);
}

However, this doesn't "set" a "RoundingMode". Instead, it copies a value to a temporary, using that rounding mode if required and then forgets about the temporary.

That's... not doing anything effectively.

So...

How Does IO work?

operator<< invokes:

template <class Backend, expression_template_option ExpressionTemplates>
inline std::ostream& operator<<(std::ostream& os, const number<Backend, ExpressionTemplates>& r)
{
   std::streamsize d  = os.precision();
   std::string     s  = r.str(d, os.flags());
   std::streamsize ss = os.width();
   if (ss > static_cast<std::streamsize>(s.size()))
   {
      char fill = os.fill();
      if ((os.flags() & std::ios_base::left) == std::ios_base::left)
         s.append(static_cast<std::string::size_type>(ss - s.size()), fill);
      else
         s.insert(static_cast<std::string::size_type>(0), static_cast<std::string::size_type>(ss - s.size()), fill);
   }
   return os << s;
}

So far, so good. The meat is in

std::string     s  = r.str(d, os.flags());

The str(...) implementation relays from the frontend (number<>) to the backend, and ends up doing (amongst a ton of different things):

 char* ps = mpfr_get_str(0, &e, 10, static_cast<std::size_t>(digits), m_data, GMP_RNDN);

So there we have it. It's just hardcoded.

If you are willing to make some assumptions, you can write your own simplified output routine (see below), but in reality I doubt that you care so much about rounding if it's just for display.

Assuming that performance is not the prime concern (since it's string IO anyways, which is rarely fast or required to be), I'm going to assume it's more valuable to be able to actually get the rounded number, so you can print that using the existing library facilities.

What Is Rounding, Really?

Rounding in MPFR isn't what you think it is: it does not round to decimal positions. It rounds to binary digits in the representation.

Lets demonstrate this by implementing rounding the MPFR way:

Live On Coliru

#include <iostream>
#include <boost/multiprecision/mpfr.hpp>

namespace bmp = boost::multiprecision;

namespace detail {
    template <
        unsigned srcDigits10, bmp::mpfr_allocation_type srcAlloc,
        unsigned dstDigits10, bmp::mpfr_allocation_type dstAlloc
    >
    void round(
            bmp::mpfr_float_backend<srcDigits10, srcAlloc> const& src, 
            mpfr_rnd_t r,
            bmp::mpfr_float_backend<dstDigits10, dstAlloc>& dst)
    {
        mpfr_set(dst.data(), src.data(), r);
    }

    template <unsigned dstDigits10, unsigned srcDigits10, bmp::mpfr_allocation_type alloc>
    auto round(bmp::mpfr_float_backend<srcDigits10, alloc> const& src, mpfr_rnd_t r) {
        bmp::mpfr_float_backend<dstDigits10, alloc> dst;
        round(src, r, dst);
        return dst;
    }
}

template <unsigned dstDigits10, typename Number>
auto round(Number const& src, mpfr_rnd_t r) {
    auto dst = detail::round<dstDigits10>(src.backend(), r);
    return bmp::number<decltype(dst)>(dst);
}

int main() {
    using bmp::mpfr_float;
    mpfr_float::default_precision(50);
    mpfr_float const a("1.55");

    std::cout << std::setprecision(20) << std::fixed;

    for (mpfr_rnd_t r : {
             MPFR_RNDN, /* round to nearest, with ties to even */
             MPFR_RNDZ, /* round toward zero */
             MPFR_RNDU, /* round toward +Inf */
             MPFR_RNDD, /* round toward -Inf */
             MPFR_RNDA, /* round away from zero */
             MPFR_RNDF, /* faithful rounding */
             MPFR_RNDNA, /* round to nearest, with ties away from zero (mpfr_round) */
         })
    {
        std::cout << round<2>(a, r) << std::endl;
    }
}

Prints

1.54687500000000000000
1.54687500000000000000
1.55468750000000000000
1.54687500000000000000
1.55468750000000000000
1.55468750000000000000
1.55468750000000000000

So, what do we expect?

Let's re-implement the stringification:

Live On Coliru

#include <iostream>
#include <boost/multiprecision/mpfr.hpp>

namespace bmp = boost::multiprecision;

template <unsigned srcDigits10, bmp::mpfr_allocation_type alloc>
auto to_string(bmp::mpfr_float_backend<srcDigits10, alloc> const& src,
               unsigned digits, mpfr_rnd_t r, std::ios::fmtflags fmtflags) {
    std::streamsize org_digits(digits);
    std::string result;

    mpfr_exp_t e = 0;
    char* ps = mpfr_get_str(0, &e, 10, static_cast<std::size_t>(digits),
                            src.data(), r);
    --e; // To match with what our formatter expects.
    if (e != -1) {
        // Oops we actually need a different number of digits to what we asked
        // for:
        mpfr_free_str(ps);
        digits += e + 1;
        if (digits == 0) {
            // We need to get *all* the digits and then possibly round up,
            // we end up with either "0" or "1" as the result.
            ps = mpfr_get_str(0, &e, 10, 0, src.data(), r);
            --e;
            unsigned offset = *ps == '-' ? 1 : 0;
            if (ps[offset] > '5') {
                ++e;
                ps[offset] = '1';
                ps[offset + 1] = 0;
            } else if (ps[offset] == '5') {
                unsigned i = offset + 1;
                bool round_up = false;
                while (ps[i] != 0) {
                    if (ps[i] != '0') {
                        round_up = true;
                        break;
                    }
                    ++i;
                }
                if (round_up) {
                    ++e;
                    ps[offset] = '1';
                    ps[offset + 1] = 0;
                } else {
                    ps[offset] = '0';
                    ps[offset + 1] = 0;
                }
            } else {
                ps[offset] = '0';
                ps[offset + 1] = 0;
            }
        } else if (digits > 0) {
            mp_exp_t old_e = e;
            ps = mpfr_get_str(0, &e, 10, static_cast<std::size_t>(digits),
                              src.data(), r);
            --e; // To match with what our formatter expects.
            if (old_e > e) {
                // in some cases, when we ask for more digits of precision, it
                // will change the number of digits to the left of the decimal,
                // if that happens, account for it here. example: cout << fixed
                // << setprecision(3) << mpf_float_50("99.9809")
                digits -= old_e - e;
                ps = mpfr_get_str(0, &e, 10, static_cast<std::size_t>(digits),
                                  src.data(), r);
                --e; // To match with what our formatter expects.
            }
        } else {
            ps = mpfr_get_str(0, &e, 10, 1, src.data(), r);
            --e;
            unsigned offset = *ps == '-' ? 1 : 0;
            ps[offset] = '0';
            ps[offset + 1] = 0;
        }
    }
    result = ps ? ps : "0";
    if (ps)
        mpfr_free_str(ps);
    bmp::detail::format_float_string(result, e, org_digits, fmtflags,
                                     0 != mpfr_zero_p(src.data()));
    return result;
}

template <unsigned srcDigits10, bmp::mpfr_allocation_type alloc>
auto to_string(
    bmp::number<bmp::mpfr_float_backend<srcDigits10, alloc>> const& src,
    unsigned digits, mpfr_rnd_t r,
    std::ios::fmtflags fmtflags = std::ios::fixed) {
    return to_string(src.backend(), digits, r, fmtflags);
}

int main() {
    using bmp::mpfr_float;
    mpfr_float::default_precision(50);
    mpfr_float const a("1.55");

    std::cout << std::setprecision(20) << std::fixed;

    for (mpfr_rnd_t r : {
             MPFR_RNDN, /* round to nearest, with ties to even */
             MPFR_RNDZ, /* round toward zero */
             MPFR_RNDU, /* round toward +Inf */
             MPFR_RNDD, /* round toward -Inf */
             MPFR_RNDA, /* round away from zero */
             MPFR_RNDF, /* faithful rounding */
             MPFR_RNDNA, /* round to nearest, with ties away from zero (mpfr_round) */
         })
    {
        std::cout
            << " -- " << to_string(a, 2, r)
            << ", " << to_string(a, 1, r)
            << " -- "  <<  to_string(a, 2, r, std::ios::scientific)
            << ", "  <<  to_string(a, 1, r, std::ios::scientific) << std::endl;
    }
}

Prints

 -- 1.55, 1.5 -- 1.55e+00, 1.5e+00
 -- 1.54, 1.5 -- 1.54e+00, 1.5e+00
 -- 1.55, 1.6 -- 1.55e+00, 1.6e+00
 -- 1.54, 1.5 -- 1.54e+00, 1.5e+00
 -- 1.55, 1.6 -- 1.55e+00, 1.6e+00
 -- 1.55, 1.5 -- 1.55e+00, 1.5e+00
 -- 1.55, 1.5 -- 1.55e+00, 1.5e+00

Disclaimer: I initially dropped some code to drop scientific notation support, so things might not be 100% up to par. Also, not tested with subnormal, infinity, nan. YMMV

Closing Thoughts

If indeed it wasn't about presentation but rounding the numbers in memory, you might construct a new mpfr_float from string representation.

My expectiation in that case would be really you rather wanted a decimal float (cpp_dec_float e.g.) in the first place.