real.hpp

/* Copyright 2009-2017 Francesco Biscani (bluescarni@gmail.com)

This file is part of the Piranha library.

The Piranha library is free software; you can redistribute it and/or modify
it under the terms of either:

  * the GNU Lesser General Public License as published by the Free
    Software Foundation; either version 3 of the License, or (at your
    option) any later version.

or

  * the GNU General Public License as published by the Free Software
    Foundation; either version 3 of the License, or (at your option) any
    later version.

or both in parallel, as here.

The Piranha library is distributed in the hope that it will be useful, but
WITHOUT ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
for more details.

You should have received copies of the GNU General Public License and the
GNU Lesser General Public License along with the Piranha library.  If not,
see https://www.gnu.org/licenses/. */

#ifndef PIRANHA_REAL_HPP
#define PIRANHA_REAL_HPP

#include <algorithm>
#include <array>
#include <boost/lexical_cast.hpp>
#include <cmath>
#include <cstddef>
#include <cstdint>
#include <iostream>
#include <limits>
#include <memory>
#include <sstream>
#include <stdexcept>
#include <string>
#include <type_traits>

#include <piranha/binomial.hpp>
#include <piranha/config.hpp>
#include <piranha/detail/demangle.hpp>
#include <piranha/detail/mpfr.hpp>
#include <piranha/detail/real_fwd.hpp>
#include <piranha/exceptions.hpp>
#include <piranha/is_cf.hpp>
#include <piranha/math.hpp>
#include <piranha/mp_integer.hpp>
#include <piranha/mp_rational.hpp>
#include <piranha/pow.hpp>
#include <piranha/s11n.hpp>
#include <piranha/safe_cast.hpp>
#include <piranha/type_traits.hpp>

namespace piranha
{

inline namespace impl
{

template <typename = int>
struct real_base {
    // Default rounding mode.
    // All operations will use the MPFR_RNDN (round to nearest) rounding mode.
    static const ::mpfr_rnd_t default_rnd = MPFR_RNDN;
    // Default significand precision.
    // The precision is the number of bits used to represent the significand of a floating-point number.
    // This default value is equivalent to the IEEE 754 quadruple-precision binary floating-point format.
    static const ::mpfr_prec_t default_prec = 113;
};

template <typename T>
const ::mpfr_rnd_t real_base<T>::default_rnd;

template <typename T>
const ::mpfr_prec_t real_base<T>::default_prec;

// Types interoperable with real.
template <typename T>
struct is_real_interoperable_type
    : disjunction<mppp::mppp_impl::is_supported_interop<T>, is_mp_integer<T>, is_mp_rational<T>> {
};
}


// NOTES:
// - if we overhaul the tests, put random precision values as well.
// - maybe we should have a setter as well for the global default precision. It would need to be an atomic
//   variable, and we need perf measures to understand the performance impact of this.
// - For series evaluation, we need to be careful performance-wise with the possible conversions that might go
//   on when mixing real with other types. E.g., pow(real,int) when evaluating polynomials. We need to make sure
//   the conversions are as fast as possible.
// - At the moment, this class is technically not sortable because moved-from reals cannot be compared. For use in
//   std::sort, we should add special casing for moved-from objects. See:
//   http://stackoverflow.com/questions/26579132/what-is-the-post-condition-of-a-move-constructor
class real : public real_base<>
{
    // Shortcut for interop type detector.
    template <typename T>
    using is_interoperable_type = is_real_interoperable_type<T>;
    // Enabler for generic ctor.
    template <typename T>
    using generic_ctor_enabler = typename std::enable_if<is_interoperable_type<T>::value, int>::type;
    // Enabler for conversion operator.
    template <typename T>
    using cast_enabler = generic_ctor_enabler<T>;
    // Enabler for in-place arithmetic operations with interop on the left.
    template <typename T>
    using generic_in_place_enabler
        = enable_if_t<conjunction<is_interoperable_type<T>, negation<std::is_const<T>>>::value, int>;
    // Precision check.
    static void prec_check(const ::mpfr_prec_t &prec)
    {
        if (unlikely(prec < MPFR_PREC_MIN || prec > MPFR_PREC_MAX)) {
            piranha_throw(std::invalid_argument, "invalid significand precision requested");
        }
    }
    static bool is_digit(char c)
    {
        // NOTE: check this answer:
        // http://stackoverflow.com/questions/13827180/char-ascii-relation
        // """
        // The mapping of integer values for characters does have one guarantee given
        // by the Standard: the values of the decimal digits are continguous.
        // (i.e., '1' - '0' == 1, ... '9' - '0' == 9)
        // """
        // It should be possible to implement this with a binary search.
        const char digits[] = "0123456789";
        return std::find(digits, digits + 10, c) != (digits + 10);
    }
    // Construction.
    void construct_from_string(const char *str, const ::mpfr_prec_t &prec)
    {
        prec_check(prec);
        ::mpfr_init2(m_value, prec);
        const int retval = ::mpfr_set_str(m_value, str, 10, default_rnd);
        if (unlikely(retval != 0)) {
            ::mpfr_clear(m_value);
            piranha_throw(std::invalid_argument, "invalid string input for real");
        }
    }
    // The idea here is that we use the largest integral and fp types supported by the MPFR api for construction,
    // and down-cast as needed.
    template <typename T, typename std::enable_if<std::is_floating_point<T>::value, int>::type = 0>
    void construct_from_generic(const T &x)
    {
        ::mpfr_set_ld(m_value, static_cast<long double>(x), default_rnd);
    }
    template <typename T, enable_if_t<conjunction<std::is_integral<T>, std::is_signed<T>>::value, int> = 0>
    void construct_from_generic(const T &si)
    {
        ::mpfr_set_sj(m_value, static_cast<std::intmax_t>(si), default_rnd);
    }
    template <typename T, enable_if_t<conjunction<std::is_integral<T>, std::is_unsigned<T>>::value, int> = 0>
    void construct_from_generic(const T &ui)
    {
        ::mpfr_set_uj(m_value, static_cast<std::uintmax_t>(ui), default_rnd);
    }
    template <typename T, typename std::enable_if<is_mp_integer<T>::value, int>::type = 0>
    void construct_from_generic(const T &n)
    {
        auto v = n.get_mpz_view();
        ::mpfr_set_z(m_value, v, default_rnd);
    }
    template <typename T, typename std::enable_if<is_mp_rational<T>::value, int>::type = 0>
    void construct_from_generic(const T &q)
    {
        auto v = q.get_mpq_view();
        ::mpfr_set_q(m_value, v, default_rnd);
    }
    // Assignment.
    void assign_from_string(const char *str)
    {
        piranha_assert(m_value->_mpfr_d);
        const int retval = ::mpfr_set_str(m_value, str, 10, default_rnd);
        if (retval != 0) {
            // Reset the internal value, as it might have been changed by ::mpfr_set_str().
            ::mpfr_set_zero(m_value, 0);
            piranha_throw(std::invalid_argument, "invalid string input for real");
        }
    }
    // Conversion.
    template <typename T>
    typename std::enable_if<std::is_same<T, bool>::value, T>::type convert_to_impl() const
    {
        return (sign() != 0);
    }
    template <typename T>
    typename std::enable_if<is_mp_integer<T>::value, T>::type convert_to_impl() const
    {
        if (is_nan() || is_inf()) {
            piranha_throw(std::overflow_error, "cannot convert non-finite real to an integral value");
        }
        T retval;
        retval.promote();
        // Explicitly request rounding to zero in this case.
        ::mpfr_get_z(&retval._get_union().g_dy(), m_value, MPFR_RNDZ);
        // NOTE: demote candidate.
        return retval;
    }
    template <typename T>
    typename std::enable_if<std::is_integral<T>::value && !std::is_same<T, bool>::value, T>::type
    convert_to_impl() const
    {
        // NOTE: of course, this can be optimised by avoiding going through the integer conversion and
        // using directly the MPFR functions.
        return static_cast<T>(static_cast<integer>(*this));
    }
    template <typename T>
    typename std::enable_if<std::is_floating_point<T>::value, T>::type convert_to_impl() const
    {
        if (is_nan()) {
            if (std::numeric_limits<T>::has_quiet_NaN) {
                return std::numeric_limits<T>::quiet_NaN();
            } else {
                piranha_throw(std::overflow_error, "cannot convert NaN to floating-point type");
            }
        }
        if (is_inf()) {
            piranha_assert(sign() != 0);
            if (std::numeric_limits<T>::has_infinity && sign() > 0) {
                return std::numeric_limits<T>::infinity();
            } else if (std::numeric_limits<T>::has_infinity && sign() < 0) {
                return std::copysign(std::numeric_limits<T>::infinity(), std::numeric_limits<T>::lowest());
            } else {
                piranha_throw(std::overflow_error, "cannot convert infinity to floating-point type");
            }
        }
        if (std::is_same<T, long double>::value) {
            return static_cast<T>(::mpfr_get_ld(m_value, default_rnd));
        }
        if (std::is_same<T, double>::value) {
            return static_cast<T>(::mpfr_get_d(m_value, default_rnd));
        }
        return static_cast<T>(::mpfr_get_flt(m_value, default_rnd));
    }
    // Smart pointer to handle the string output from mpfr.
    typedef std::unique_ptr<char, void (*)(char *)> smart_mpfr_str;
    template <typename T>
    typename std::enable_if<is_mp_rational<T>::value, T>::type convert_to_impl() const
    {
        if (is_nan()) {
            piranha_throw(std::overflow_error, "cannot convert NaN to rational");
        }
        if (is_inf()) {
            piranha_throw(std::overflow_error, "cannot convert infinity to rational");
        }
        if (sign() == 0) {
            return T{};
        }
        // Get string representation.
        ::mpfr_exp_t exp(0);
        char *cptr = ::mpfr_get_str(nullptr, &exp, 10, 0, m_value, default_rnd);
        if (unlikely(!cptr)) {
            piranha_throw(std::overflow_error,
                          "error in conversion of real to rational: the call to the MPFR function failed");
        }
        smart_mpfr_str str_ptr(cptr, ::mpfr_free_str);
        // Transform into fraction.
        std::size_t digits = 0u;
        for (; *cptr != '\0'; ++cptr) {
            if (is_digit(*cptr)) {
                ++digits;
            }
        }
        if (!digits) {
            piranha_throw(std::overflow_error, "error in conversion of real to rational: invalid number of digits");
        }
        // NOTE: here the only exception that can be thrown is when raising to a power
        // that cannot be represented by unsigned long.
        try {
            T retval(str_ptr.get());
            // NOTE: possible optimizations here include going through direct GMP routines.
            retval *= T(1, 10).pow(digits);
            retval *= T(10).pow(exp);
            return retval;
        } catch (...) {
            piranha_throw(std::overflow_error, "error in conversion of real to rational: exponent is too large");
        }
    }
    // In-place addition.
    // NOTE: all sorts of optimisations, here and in binary add, are possible (e.g., steal from rvalue ref,
    // avoid setting precision twice in binary operators, etc.). For the moment we keep it basic.
    real &in_place_add(const real &r)
    {
        if (r.get_prec() > get_prec()) {
            // Re-init this with the prec of r.
            *this = real{*this, r.get_prec()};
        }
        ::mpfr_add(m_value, m_value, r.m_value, default_rnd);
        return *this;
    }
    template <typename T, typename std::enable_if<is_mp_rational<T>::value, int>::type = 0>
    real &in_place_add(const T &q)
    {
        auto v = q.get_mpq_view();
        ::mpfr_add_q(m_value, m_value, v, default_rnd);
        return *this;
    }
    template <typename T, typename std::enable_if<is_mp_integer<T>::value, int>::type = 0>
    real &in_place_add(const T &n)
    {
        auto v = n.get_mpz_view();
        ::mpfr_add_z(m_value, m_value, v, default_rnd);
        return *this;
    }
    // NOTE: possible optimisations here.
    template <typename T, typename std::enable_if<std::is_integral<T>::value, int>::type = 0>
    real &in_place_add(const T &n)
    {
        return in_place_add(integer(n));
    }
    // NOTE: possible optimisations here as well.
    template <typename T, typename std::enable_if<std::is_floating_point<T>::value, int>::type = 0>
    real &in_place_add(const T &x)
    {
        // Construct real with the same precision as this, then add.
        return in_place_add(real{x, get_prec()});
    }
    // Binary add.
    static real binary_add(const real &a, const real &b)
    {
        real retval{a};
        retval += b;
        return retval;
    }
    // Single implementation for all interoperable types.
    template <typename T, typename std::enable_if<is_interoperable_type<T>::value, int>::type = 0>
    static real binary_add(const real &a, const T &b)
    {
        real retval{a};
        retval += b;
        return retval;
    }
    template <typename T, typename std::enable_if<is_interoperable_type<T>::value, int>::type = 0>
    static real binary_add(const T &a, const real &b)
    {
        return binary_add(b, a);
    }
    // In-place subtraction.
    real &in_place_sub(const real &r)
    {
        if (r.get_prec() > get_prec()) {
            *this = real{*this, r.get_prec()};
        }
        ::mpfr_sub(m_value, m_value, r.m_value, default_rnd);
        return *this;
    }
    template <typename T, typename std::enable_if<is_mp_rational<T>::value, int>::type = 0>
    real &in_place_sub(const T &q)
    {
        auto v = q.get_mpq_view();
        ::mpfr_sub_q(m_value, m_value, v, default_rnd);
        return *this;
    }
    template <typename T, typename std::enable_if<is_mp_integer<T>::value, int>::type = 0>
    real &in_place_sub(const T &n)
    {
        auto v = n.get_mpz_view();
        ::mpfr_sub_z(m_value, m_value, v, default_rnd);
        return *this;
    }
    template <typename T, typename std::enable_if<std::is_integral<T>::value, int>::type = 0>
    real &in_place_sub(const T &n)
    {
        return in_place_sub(integer(n));
    }
    template <typename T, typename std::enable_if<std::is_floating_point<T>::value, int>::type = 0>
    real &in_place_sub(const T &x)
    {
        return in_place_sub(real{x, get_prec()});
    }
    // Binary sub.
    static real binary_sub(const real &a, const real &b)
    {
        real retval{a};
        retval -= b;
        return retval;
    }
    template <typename T, typename std::enable_if<is_interoperable_type<T>::value, int>::type = 0>
    static real binary_sub(const real &a, const T &b)
    {
        real retval{a};
        retval -= b;
        return retval;
    }
    template <typename T, typename std::enable_if<is_interoperable_type<T>::value, int>::type = 0>
    static real binary_sub(const T &a, const real &b)
    {
        auto retval = binary_sub(b, a);
        retval.negate();
        return retval;
    }
    // In-place multiplication.
    real &in_place_mul(const real &r)
    {
        if (r.get_prec() > get_prec()) {
            *this = real{*this, r.get_prec()};
        }
        ::mpfr_mul(m_value, m_value, r.m_value, default_rnd);
        return *this;
    }
    template <typename T, typename std::enable_if<is_mp_rational<T>::value, int>::type = 0>
    real &in_place_mul(const T &q)
    {
        auto v = q.get_mpq_view();
        ::mpfr_mul_q(m_value, m_value, v, default_rnd);
        return *this;
    }
    template <typename T, typename std::enable_if<is_mp_integer<T>::value, int>::type = 0>
    real &in_place_mul(const T &n)
    {
        auto v = n.get_mpz_view();
        ::mpfr_mul_z(m_value, m_value, v, default_rnd);
        return *this;
    }
    template <typename T, typename std::enable_if<std::is_integral<T>::value, int>::type = 0>
    real &in_place_mul(const T &n)
    {
        return in_place_mul(integer(n));
    }
    template <typename T, typename std::enable_if<std::is_floating_point<T>::value, int>::type = 0>
    real &in_place_mul(const T &x)
    {
        return in_place_mul(real{x, get_prec()});
    }
    // Binary mul.
    static real binary_mul(const real &a, const real &b)
    {
        real retval{a};
        retval *= b;
        return retval;
    }
    template <typename T, typename std::enable_if<is_interoperable_type<T>::value, int>::type = 0>
    static real binary_mul(const real &a, const T &b)
    {
        real retval{a};
        retval *= b;
        return retval;
    }
    template <typename T, typename std::enable_if<is_interoperable_type<T>::value, int>::type = 0>
    static real binary_mul(const T &a, const real &b)
    {
        return binary_mul(b, a);
    }
    // In-place division.
    real &in_place_div(const real &r)
    {
        if (r.get_prec() > get_prec()) {
            *this = real{*this, r.get_prec()};
        }
        ::mpfr_div(m_value, m_value, r.m_value, default_rnd);
        return *this;
    }
    template <typename T, typename std::enable_if<is_mp_rational<T>::value, int>::type = 0>
    real &in_place_div(const T &q)
    {
        auto v = q.get_mpq_view();
        ::mpfr_div_q(m_value, m_value, v, default_rnd);
        return *this;
    }
    template <typename T, typename std::enable_if<is_mp_integer<T>::value, int>::type = 0>
    real &in_place_div(const T &n)
    {
        auto v = n.get_mpz_view();
        ::mpfr_div_z(m_value, m_value, v, default_rnd);
        return *this;
    }
    template <typename T, typename std::enable_if<std::is_integral<T>::value, int>::type = 0>
    real &in_place_div(const T &n)
    {
        return in_place_div(integer(n));
    }
    template <typename T, typename std::enable_if<std::is_floating_point<T>::value, int>::type = 0>
    real &in_place_div(const T &x)
    {
        return in_place_div(real{x, get_prec()});
    }
    // Binary div.
    static real binary_div(const real &a, const real &b)
    {
        real retval{a};
        retval /= b;
        return retval;
    }
    template <typename T, typename std::enable_if<is_interoperable_type<T>::value, int>::type = 0>
    static real binary_div(const real &a, const T &b)
    {
        real retval{a};
        retval /= b;
        return retval;
    }
    template <typename T, typename std::enable_if<is_interoperable_type<T>::value, int>::type = 0>
    static real binary_div(const T &a, const real &b)
    {
        // Create with same precision as b.
        real retval{a, b.get_prec()};
        retval /= b;
        return retval;
    }
    // Equality.
    static bool binary_equality(const real &r1, const real &r2)
    {
        return (::mpfr_equal_p(r1.m_value, r2.m_value) != 0);
    }
    template <typename T, typename std::enable_if<is_mp_integer<T>::value, int>::type = 0>
    static bool binary_equality(const real &r, const T &n)
    {
        if (r.is_nan()) {
            return false;
        }
        auto v = n.get_mpz_view();
        return (::mpfr_cmp_z(r.m_value, v) == 0);
    }
    template <typename T, typename std::enable_if<is_mp_rational<T>::value, int>::type = 0>
    static bool binary_equality(const real &r, const T &q)
    {
        if (r.is_nan()) {
            return false;
        }
        auto v = q.get_mpq_view();
        return (::mpfr_cmp_q(r.m_value, v) == 0);
    }
    template <typename T, typename std::enable_if<std::is_integral<T>::value, int>::type = 0>
    static bool binary_equality(const real &r, const T &n)
    {
        if (r.is_nan()) {
            return false;
        }
        return r == integer(n);
    }
    template <typename T, typename std::enable_if<std::is_floating_point<T>::value, int>::type = 0>
    static bool binary_equality(const real &r, const T &x)
    {
        if (r.is_nan() || std::isnan(x)) {
            return false;
        }
        return r == real{x, r.get_prec()};
    }
    // NOTE: this is the reverse of above.
    template <typename T, typename std::enable_if<is_interoperable_type<T>::value, int>::type = 0>
    static bool binary_equality(const T &x, const real &r)
    {
        return binary_equality(r, x);
    }
    // Binary less-than.
    static bool binary_less_than(const real &r1, const real &r2)
    {
        return (::mpfr_less_p(r1.m_value, r2.m_value) != 0);
    }
    template <typename T, typename std::enable_if<is_mp_rational<T>::value, int>::type = 0>
    static bool binary_less_than(const real &r, const T &q)
    {
        auto v = q.get_mpq_view();
        return (::mpfr_cmp_q(r.m_value, v) < 0);
    }
    template <typename T, typename std::enable_if<is_mp_integer<T>::value, int>::type = 0>
    static bool binary_less_than(const real &r, const T &n)
    {
        auto v = n.get_mpz_view();
        return (::mpfr_cmp_z(r.m_value, v) < 0);
    }
    template <typename T, typename std::enable_if<std::is_integral<T>::value, int>::type = 0>
    static bool binary_less_than(const real &r, const T &n)
    {
        return r < integer(n);
    }
    template <typename T, typename std::enable_if<std::is_floating_point<T>::value, int>::type = 0>
    static bool binary_less_than(const real &r, const T &x)
    {
        return r < real{x, r.get_prec()};
    }
    // Binary less-than or equal.
    static bool binary_leq(const real &r1, const real &r2)
    {
        return (::mpfr_lessequal_p(r1.m_value, r2.m_value) != 0);
    }
    template <typename T, typename std::enable_if<is_mp_rational<T>::value, int>::type = 0>
    static bool binary_leq(const real &r, const T &q)
    {
        auto v = q.get_mpq_view();
        return (::mpfr_cmp_q(r.m_value, v) <= 0);
    }
    template <typename T, typename std::enable_if<is_mp_integer<T>::value, int>::type = 0>
    static bool binary_leq(const real &r, const T &n)
    {
        auto v = n.get_mpz_view();
        return (::mpfr_cmp_z(r.m_value, v) <= 0);
    }
    template <typename T, typename std::enable_if<std::is_integral<T>::value, int>::type = 0>
    static bool binary_leq(const real &r, const T &n)
    {
        return r <= integer(n);
    }
    template <typename T, typename std::enable_if<std::is_floating_point<T>::value, int>::type = 0>
    static bool binary_leq(const real &r, const T &x)
    {
        return r <= real{x, r.get_prec()};
    }
    // Inverse forms of less-than and leq.
    template <typename T, typename std::enable_if<is_interoperable_type<T>::value, int>::type = 0>
    static bool binary_less_than(const T &x, const real &r)
    {
        return !binary_leq(r, x);
    }
    template <typename T, typename std::enable_if<is_interoperable_type<T>::value, int>::type = 0>
    static bool binary_leq(const T &x, const real &r)
    {
        return !binary_less_than(r, x);
    }
    // NOTE: we need to handle separately the NaNs as we cannot resort to the inversion of the comparison operators for
    // them.
    static bool check_nan(const real &r)
    {
        return r.is_nan();
    }
    template <typename T>
    static bool check_nan(const T &x, typename std::enable_if<std::is_floating_point<T>::value>::type * = nullptr)
    {
        return std::isnan(x);
    }
    template <typename T>
    static bool check_nan(const T &, typename std::enable_if<!std::is_floating_point<T>::value>::type * = nullptr)
    {
        return false;
    }
    template <typename T, typename U>
    static bool is_nan_comparison(const T &a, const U &b)
    {
        return (check_nan(a) || check_nan(b));
    }
    // Serialization support.
    friend class boost::serialization::access;
    // Utility function to infer the size (in number of limbs) from the precision.
    static ::mpfr_prec_t size_from_prec(::mpfr_prec_t prec)
    {
        const ::mpfr_prec_t q = prec / ::mp_bits_per_limb, r = prec % ::mp_bits_per_limb;
        return q + (r != 0);
    }
    // Portable s11n.
    template <class Archive, enable_if_t<!std::is_same<Archive, boost::archive::binary_oarchive>::value, int> = 0>
    void save(Archive &ar, unsigned) const
    {
        // NOTE: like in mp_integer, the performance here can be improved significantly.
        std::ostringstream oss;
        oss << *this;
        auto prec = get_prec();
        auto s = oss.str();
        piranha::boost_save(ar, prec);
        piranha::boost_save(ar, s);
    }
    template <class Archive, enable_if_t<!std::is_same<Archive, boost::archive::binary_iarchive>::value, int> = 0>
    void load(Archive &ar, unsigned)
    {
        ::mpfr_prec_t prec;
        PIRANHA_MAYBE_TLS std::string s;
        piranha::boost_load(ar, prec);
        piranha::boost_load(ar, s);
        *this = real(s, prec);
    }
    // Binary s11n.
    template <class Archive, enable_if_t<std::is_same<Archive, boost::archive::binary_oarchive>::value, int> = 0>
    void save(Archive &ar, unsigned) const
    {
        piranha::boost_save(ar, m_value->_mpfr_prec);
        piranha::boost_save(ar, m_value->_mpfr_sign);
        piranha::boost_save(ar, m_value->_mpfr_exp);
        const ::mpfr_prec_t s = size_from_prec(m_value->_mpfr_prec);
        // NOTE: no need to save the size, as it can be recovered from the prec.
        for (::mpfr_prec_t i = 0; i < s; ++i) {
            piranha::boost_save(ar, m_value->_mpfr_d[i]);
        }
    }
    template <class Archive, enable_if_t<std::is_same<Archive, boost::archive::binary_iarchive>::value, int> = 0>
    void load(Archive &ar, unsigned)
    {
        // First we recover the non-limb members.
        ::mpfr_prec_t prec;
        decltype(m_value->_mpfr_sign) sign;
        decltype(m_value->_mpfr_exp) exp;
        piranha::boost_load(ar, prec);
        piranha::boost_load(ar, sign);
        piranha::boost_load(ar, exp);
        // Recover the size in limbs from prec.
        const ::mpfr_prec_t s = size_from_prec(prec);
        // Set the precision.
        set_prec(prec);
        piranha_assert(m_value->_mpfr_prec == prec);
        m_value->_mpfr_sign = sign;
        m_value->_mpfr_exp = exp;
        try {
            // NOTE: protect in try/catch as in theory boost_load() could throw even
            // in case of valid archive (e.g., memory errors maybe?) and we want
            // to deal with this case.
            for (::mpfr_prec_t i = 0; i < s; ++i) {
                piranha::boost_load(ar, *(m_value->_mpfr_d + i));
            }
        } catch (...) {
            // Set to zero before re-throwing.
            ::mpfr_set_ui(m_value, 0u, default_rnd);
            throw;
        }
    }
    BOOST_SERIALIZATION_SPLIT_MEMBER()
public:

    real()
    {
        ::mpfr_init2(m_value, default_prec);
        ::mpfr_set_zero(m_value, 0);
    }

    real(const real &other)
    {
        // Init with the same precision as other, and then set.
        ::mpfr_init2(m_value, other.get_prec());
        ::mpfr_set(m_value, other.m_value, default_rnd);
    }

    real(real &&other) noexcept
    {
        m_value->_mpfr_prec = other.m_value->_mpfr_prec;
        m_value->_mpfr_sign = other.m_value->_mpfr_sign;
        m_value->_mpfr_exp = other.m_value->_mpfr_exp;
        m_value->_mpfr_d = other.m_value->_mpfr_d;
        // Erase other.
        other.m_value->_mpfr_prec = 0;
        other.m_value->_mpfr_sign = 0;
        other.m_value->_mpfr_exp = 0;
        other.m_value->_mpfr_d = nullptr;
    }

    explicit real(const char *str, const ::mpfr_prec_t &prec = default_prec)
    {
        construct_from_string(str, prec);
    }

    explicit real(const std::string &str, const ::mpfr_prec_t &prec = default_prec)
    {
        construct_from_string(str.c_str(), prec);
    }

    explicit real(const real &other, const ::mpfr_prec_t &prec)
    {
        prec_check(prec);
        ::mpfr_init2(m_value, prec);
        ::mpfr_set(m_value, other.m_value, default_rnd);
    }

    template <typename T, typename = generic_ctor_enabler<T>>
    explicit real(const T &x, const ::mpfr_prec_t &prec = default_prec)
    {
        prec_check(prec);
        ::mpfr_init2(m_value, prec);
        construct_from_generic(x);
    }

    ~real();

    real &operator=(const real &other)
    {
        if (this != &other) {
            // Handle assignment to moved-from objects.
            if (m_value->_mpfr_d) {
                // Copy the precision. This will also reset the internal value.
                set_prec(other.get_prec());
            } else {
                piranha_assert(!m_value->_mpfr_prec && !m_value->_mpfr_sign && !m_value->_mpfr_exp);
                // Reinit before setting.
                ::mpfr_init2(m_value, other.get_prec());
            }
            ::mpfr_set(m_value, other.m_value, default_rnd);
        }
        return *this;
    }

    real &operator=(real &&other) noexcept
    {
        // NOTE: swap() already has the check for this.
        swap(other);
        return *this;
    }

    real &operator=(const std::string &str)
    {
        return operator=(str.c_str());
    }

    real &operator=(const char *str)
    {
        // Handle moved-from objects.
        if (m_value->_mpfr_d) {
            assign_from_string(str);
        } else {
            piranha_assert(!m_value->_mpfr_prec && !m_value->_mpfr_sign && !m_value->_mpfr_exp);
            construct_from_string(str, default_prec);
        }
        return *this;
    }

    template <typename T, typename = generic_ctor_enabler<T>>
    real &operator=(const T &x)
    {
        if (!m_value->_mpfr_d) {
            piranha_assert(!m_value->_mpfr_prec && !m_value->_mpfr_sign && !m_value->_mpfr_exp);
            // Re-init with default prec if it was moved-from.
            ::mpfr_init2(m_value, default_prec);
        }
        // NOTE: all construct_from_generic() methods here are really assignments.
        construct_from_generic(x);
        return *this;
    }

    template <typename T, typename = cast_enabler<T>>
    explicit operator T() const
    {
        return convert_to_impl<T>();
    }

    void swap(real &other)
    {
        if (this == &other) {
            return;
        }
        ::mpfr_swap(m_value, other.m_value);
    }

    int sign() const
    {
        if (is_nan()) {
            return 0;
        } else {
            return mpfr_sgn(m_value);
        }
    }

    bool is_zero() const
    {
        return (mpfr_zero_p(m_value) != 0);
    }

    bool is_nan() const
    {
        return mpfr_nan_p(m_value) != 0;
    }

    bool is_inf() const
    {
        return mpfr_inf_p(m_value) != 0;
    }

    ::mpfr_prec_t get_prec() const
    {
        return mpfr_get_prec(m_value);
    }

    void set_prec(const ::mpfr_prec_t &prec)
    {
        prec_check(prec);
        ::mpfr_set_prec(m_value, prec);
    }

    void negate()
    {
        ::mpfr_neg(m_value, m_value, default_rnd);
    }

    void truncate()
    {
        if (is_inf() || is_nan()) {
            return;
        }
        ::mpfr_trunc(m_value, m_value);
    }

    real truncated() const
    {
        real retval{*this};
        retval.truncate();
        return retval;
    }

    template <typename T>
    auto operator+=(const T &x) -> decltype(this->in_place_add(x))
    {
        return in_place_add(x);
    }

    template <typename T, generic_in_place_enabler<T> = 0>
    friend T &operator+=(T &x, const real &r)
    {
        // NOTE: for the supported types, move assignment can never throw.
        return x = static_cast<T>(r + x);
    }

    template <typename T, typename U>
    friend auto operator+(const T &x, const U &y) -> decltype(real::binary_add(x, y))
    {
        return binary_add(x, y);
    }

    real operator+() const
    {
        return *this;
    }

    real &operator++()
    {
        return operator+=(1);
    }

    real operator++(int)
    {
        const real retval(*this);
        ++(*this);
        return retval;
    }

    template <typename T>
    auto operator-=(const T &x) -> decltype(this->in_place_sub(x))
    {
        return in_place_sub(x);
    }

    template <typename T, generic_in_place_enabler<T> = 0>
    friend T &operator-=(T &x, const real &r)
    {
        return x = static_cast<T>(x - r);
    }

    template <typename T, typename U>
    friend auto operator-(const T &x, const U &y) -> decltype(real::binary_sub(x, y))
    {
        return binary_sub(x, y);
    }

    real operator-() const
    {
        real retval(*this);
        retval.negate();
        return retval;
    }

    real &operator--()
    {
        return operator-=(1);
    }

    real operator--(int)
    {
        const real retval(*this);
        --(*this);
        return retval;
    }

    template <typename T>
    auto operator*=(const T &x) -> decltype(this->in_place_mul(x))
    {
        return in_place_mul(x);
    }

    template <typename T, generic_in_place_enabler<T> = 0>
    friend T &operator*=(T &x, const real &r)
    {
        return x = static_cast<T>(x * r);
    }

    template <typename T, typename U>
    friend auto operator*(const T &x, const U &y) -> decltype(real::binary_mul(x, y))
    {
        return binary_mul(x, y);
    }

    template <typename T>
    auto operator/=(const T &x) -> decltype(this->in_place_div(x))
    {
        return in_place_div(x);
    }

    template <typename T, generic_in_place_enabler<T> = 0>
    friend T &operator/=(T &x, const real &r)
    {
        return x = static_cast<T>(x / r);
    }

    template <typename T, typename U>
    friend auto operator/(const T &x, const U &y) -> decltype(real::binary_div(x, y))
    {
        return binary_div(x, y);
    }

    real &multiply_accumulate(const real &r1, const real &r2)
    {
        const ::mpfr_prec_t prec1 = std::max<::mpfr_prec_t>(r1.get_prec(), r2.get_prec());
        if (prec1 > get_prec()) {
            *this = real{*this, prec1};
        }
// So the story here is that mpfr_fma has been reported to be slower than the two separate
// operations. Benchmarks on fateman1 indicate this is indeed the case (3.6 vs 2.7 secs
// on 4 threads). Hopefully it will be fixed in the future, for now adopt the workaround.
// http://www.loria.fr/~zimmerma/mpfr-mpc-2014.html
// NOTE: this optimisation requires the thread_local keyword.
#if defined(PIRANHA_HAVE_THREAD_LOCAL)
        static thread_local real tmp;
        // NOTE: set the same precision as this, which is now the max precision of the 3 operands.
        // If we do not do this, then tmp has an undeterminate precision. Use the raw MPFR function
        // in order to avoid the checks in get_prec(), as we know the precision has a sane value.
        ::mpfr_set_prec(tmp.m_value, mpfr_get_prec(m_value));
        ::mpfr_mul(tmp.m_value, r1.m_value, r2.m_value, MPFR_RNDN);
        ::mpfr_add(m_value, m_value, tmp.m_value, MPFR_RNDN);
#else
        ::mpfr_fma(m_value, r1.m_value, r2.m_value, m_value, default_rnd);
#endif
        return *this;
    }

    template <typename T, typename U>
    friend auto operator==(const T &x, const U &y) -> decltype(real::binary_equality(x, y))
    {
        return binary_equality(x, y);
    }

    template <typename T, typename U>
    friend auto operator!=(const T &x, const U &y) -> decltype(!real::binary_equality(x, y))
    {
        return !binary_equality(x, y);
    }

    template <typename T, typename U>
    friend auto operator<(const T &x, const U &y) -> decltype(real::binary_less_than(x, y))
    {
        if (is_nan_comparison(x, y)) {
            return false;
        }
        return binary_less_than(x, y);
    }

    template <typename T, typename U>
    friend auto operator<=(const T &x, const U &y) -> decltype(real::binary_leq(x, y))
    {
        if (is_nan_comparison(x, y)) {
            return false;
        }
        return binary_leq(x, y);
    }

    template <typename T, typename U>
    friend auto operator>(const T &x, const U &y) -> decltype(real::binary_less_than(y, x))
    {
        if (is_nan_comparison(x, y)) {
            return false;
        }
        return y < x;
    }

    template <typename T, typename U>
    friend auto operator>=(const T &x, const U &y) -> decltype(real::binary_leq(y, x))
    {
        if (is_nan_comparison(x, y)) {
            return false;
        }
        return y <= x;
    }

    real pow(const real &exp) const
    {
        real retval{0, get_prec()};
        if (exp.get_prec() > get_prec()) {
            retval.set_prec(exp.get_prec());
        }
        ::mpfr_pow(retval.m_value, m_value, exp.m_value, default_rnd);
        return retval;
    }

    real gamma() const
    {
        real retval{0, get_prec()};
        ::mpfr_gamma(retval.m_value, m_value, default_rnd);
        return retval;
    }

    real lgamma() const
    {
        real retval{0, get_prec()};
        // This is the sign of gamma(*this). We don't use this.
        int sign;
        ::mpfr_lgamma(retval.m_value, &sign, m_value, default_rnd);
        return retval;
    }

    real exp() const
    {
        real retval{0, get_prec()};
        ::mpfr_exp(retval.m_value, m_value, default_rnd);
        return retval;
    }
    real binomial(const real &) const;

    real abs() const
    {
        real retval(*this);
        ::mpfr_abs(retval.m_value, retval.m_value, default_rnd);
        return retval;
    }

    real sin() const
    {
        real retval(0, get_prec());
        ::mpfr_sin(retval.m_value, m_value, default_rnd);
        return retval;
    }

    real cos() const
    {
        real retval(0, get_prec());
        ::mpfr_cos(retval.m_value, m_value, default_rnd);
        return retval;
    }

    real pi() const
    {
        real retval(0, get_prec());
        ::mpfr_const_pi(retval.m_value, default_rnd);
        return retval;
    }

    friend std::ostream &operator<<(std::ostream &os, const real &r)
    {
        if (r.is_nan()) {
            os << "nan";
            return os;
        }
        if (r.is_inf()) {
            if (r.sign() > 0) {
                os << "inf";
            } else {
                os << "-inf";
            }
            return os;
        }
        ::mpfr_exp_t exp(0);
        char *cptr = ::mpfr_get_str(nullptr, &exp, 10, 0, r.m_value, default_rnd);
        if (unlikely(!cptr)) {
            piranha_throw(std::overflow_error,
                          "error in conversion of real to rational: the call to the MPFR function failed");
        }
        smart_mpfr_str str(cptr, ::mpfr_free_str);
        // Copy into C++ string.
        std::string cpp_str(str.get());
        // Insert the radix point.
        auto it = std::find_if(cpp_str.begin(), cpp_str.end(), is_digit);
        if (it != cpp_str.end()) {
            ++it;
            cpp_str.insert(it, '.');
            if (exp == std::numeric_limits<::mpfr_exp_t>::min()) {
                piranha_throw(std::overflow_error, "overflow in conversion of real to string");
            }
            --exp;
            if (exp != ::mpfr_exp_t(0) && r.sign() != 0) {
                cpp_str.append("e" + boost::lexical_cast<std::string>(exp));
            }
        }
        os << cpp_str;
        return os;
    }

    friend std::istream &operator>>(std::istream &is, real &r)
    {
        std::string tmp_str;
        std::getline(is, tmp_str);
        r = tmp_str;
        return is;
    }
    std::remove_extent<::mpfr_t>::type *get_mpfr_t()
    {
        return &m_value[0u];
    }
    const std::remove_extent<::mpfr_t>::type *get_mpfr_t() const
    {
        return &m_value[0u];
    }

#if defined(PIRANHA_WITH_MSGPACK)
private:
    // msgpack enabler.
    template <typename Stream>
    using msgpack_pack_enabler
        = enable_if_t<conjunction<is_msgpack_stream<Stream>, has_msgpack_pack<Stream, ::mpfr_prec_t>,
                                  has_msgpack_pack<Stream, std::string>,
                                  has_msgpack_pack<Stream, decltype(std::declval<const ::mpfr_t &>()->_mpfr_sign)>,
                                  has_msgpack_pack<Stream, decltype(std::declval<const ::mpfr_t &>()->_mpfr_exp)>,
                                  has_msgpack_pack<Stream, ::mp_limb_t>,
                                  has_safe_cast<std::uint32_t, ::mpfr_prec_t>>::value,
                      int>;
    template <typename U>
    using msgpack_convert_enabler
        = enable_if_t<conjunction<std::is_same<U, U>, // For SFINAE.
                                  has_msgpack_convert<::mpfr_prec_t>, has_msgpack_convert<std::string>,
                                  has_msgpack_convert<decltype(std::declval<const ::mpfr_t &>()->_mpfr_sign)>,
                                  has_msgpack_convert<decltype(std::declval<const ::mpfr_t &>()->_mpfr_exp)>,
                                  has_msgpack_convert<::mp_limb_t>>::value,
                      int>;

public:

    template <typename Stream, msgpack_pack_enabler<Stream> = 0>
    void msgpack_pack(msgpack::packer<Stream> &p, msgpack_format f) const
    {
        if (f == msgpack_format::portable) {
            std::ostringstream oss;
            oss << *this;
            auto prec = get_prec();
            auto s = oss.str();
            p.pack_array(2);
            piranha::msgpack_pack(p, prec, f);
            piranha::msgpack_pack(p, s, f);
        } else {
            // NOTE: storing both prec and the number of limbs is a bit redundant: it is possible to
            // infer the number of limbs from prec but not viceversa (only an upper/lower bound). So let's
            // store them both.
            p.pack_array(4);
            piranha::msgpack_pack(p, m_value->_mpfr_prec, f);
            piranha::msgpack_pack(p, m_value->_mpfr_sign, f);
            piranha::msgpack_pack(p, m_value->_mpfr_exp, f);
            const auto s = safe_cast<std::uint32_t>(size_from_prec(m_value->_mpfr_prec));
            p.pack_array(s);
            // NOTE: no need to save the size, as it can be recovered from the prec.
            for (std::uint32_t i = 0; i < s; ++i) {
                piranha::msgpack_pack(p, m_value->_mpfr_d[i], f);
            }
        }
    }

    template <typename U = real, msgpack_convert_enabler<U> = 0>
    void msgpack_convert(const msgpack::object &o, msgpack_format f)
    {
        if (f == msgpack_format::portable) {
            PIRANHA_MAYBE_TLS std::array<msgpack::object, 2> vobj;
            o.convert(vobj);
            ::mpfr_prec_t prec;
            PIRANHA_MAYBE_TLS std::string s;
            piranha::msgpack_convert(prec, vobj[0], f);
            piranha::msgpack_convert(s, vobj[1], f);
            *this = real(s, prec);
        } else {
            PIRANHA_MAYBE_TLS std::array<msgpack::object, 4> vobj;
            o.convert(vobj);
            // First let's handle the non-limbs members.
            ::mpfr_prec_t prec;
            decltype(m_value->_mpfr_sign) sign;
            decltype(m_value->_mpfr_exp) exp;
            piranha::msgpack_convert(prec, vobj[0], f);
            piranha::msgpack_convert(sign, vobj[1], f);
            piranha::msgpack_convert(exp, vobj[2], f);
            set_prec(prec);
            piranha_assert(m_value->_mpfr_prec == prec);
            m_value->_mpfr_sign = sign;
            m_value->_mpfr_exp = exp;
            // Next the limbs. Protect in try/catch so if anything goes wrong we can fix this in the
            // catch block before re-throwing.
            try {
                PIRANHA_MAYBE_TLS std::vector<msgpack::object> vlimbs;
                vobj[3].convert(vlimbs);
                const auto s = safe_cast<std::vector<msgpack::object>::size_type>(size_from_prec(prec));
                if (unlikely(s != vlimbs.size())) {
                    piranha_throw(std::invalid_argument,
                                  std::string("error in the msgpack deserialization of a real: the number "
                                              "of serialized limbs (")
                                      + std::to_string(vlimbs.size())
                                      + ") is not consistent with the number of limbs inferred from the precision ("
                                      + std::to_string(s) + ")");
                }
                for (decltype(vlimbs.size()) i = 0; i < s; ++i) {
                    piranha::msgpack_convert(m_value->_mpfr_d[i], vlimbs[i], f);
                }
            } catch (...) {
                // Set to zero before re-throwing.
                ::mpfr_set_ui(m_value, 0u, default_rnd);
                throw;
            }
        }
    }
#endif

private:
    ::mpfr_t m_value;
};

namespace math
{

template <typename T>
struct negate_impl<T, typename std::enable_if<std::is_same<T, real>::value>::type> {

    void operator()(real &x) const
    {
        x.negate();
    }
};

template <typename T>
struct is_zero_impl<T, typename std::enable_if<std::is_same<T, real>::value>::type> {

    bool operator()(const T &r) const
    {
        return r.is_zero();
    }
};
}

namespace detail
{

// Enabler for real pow.
template <typename T, typename U>
using real_pow_enabler =
    typename std::enable_if<(std::is_same<real, T>::value && is_real_interoperable_type<U>::value)
                            || (std::is_same<real, U>::value && is_real_interoperable_type<T>::value)
                            || (std::is_same<real, T>::value && std::is_same<real, U>::value)>::type;

// Binomial follows the same rules as pow.
template <typename T, typename U>
using real_binomial_enabler = real_pow_enabler<T, U>;
}

namespace math
{


template <typename T, typename U>
struct pow_impl<T, U, detail::real_pow_enabler<T, U>> {

    real operator()(const real &r, const real &x) const
    {
        return r.pow(x);
    }

    template <typename T2>
    real operator()(const real &r, const T2 &x) const
    {
        // NOTE: init with the same precision as r in order
        // to maintain the same precision in the result.
        return r.pow(real{x, r.get_prec()});
    }

    template <typename T2>
    real operator()(const T2 &r, const real &x) const
    {
        return real{r, x.get_prec()}.pow(x);
    }
};

template <typename T>
struct sin_impl<T, typename std::enable_if<std::is_same<T, real>::value>::type> {

    real operator()(const T &r) const
    {
        return r.sin();
    }
};

template <typename T>
struct cos_impl<T, typename std::enable_if<std::is_same<T, real>::value>::type> {

    real operator()(const T &r) const
    {
        return r.cos();
    }
};

template <typename T>
struct abs_impl<T, typename std::enable_if<std::is_same<T, real>::value>::type> {

    T operator()(const T &x) const
    {
        return x.abs();
    }
};

template <typename T>
struct partial_impl<T, typename std::enable_if<std::is_same<T, real>::value>::type> {

    real operator()(const real &, const std::string &) const
    {
        return real(0);
    }
};


template <typename T, typename U>
struct binomial_impl<T, U, detail::real_binomial_enabler<T, U>> {

    real operator()(const real &x, const real &y) const
    {
        return x.binomial(y);
    }

    template <typename T2>
    real operator()(const real &x, const T2 &y) const
    {
        // NOTE: init with the same precision as r in order
        // to maintain the same precision in the result.
        return x.binomial(real{y, x.get_prec()});
    }

    template <typename T2>
    real operator()(const T2 &x, const real &y) const
    {
        return real{x, y.get_prec()}.binomial(y);
    }
};

template <>
struct multiply_accumulate_impl<real> {

    void operator()(real &x, const real &y, const real &z) const
    {
        x.multiply_accumulate(y, z);
    }
};
}

inline real::~real()
{
    PIRANHA_TT_CHECK(is_cf, real);
    static_assert(default_prec >= MPFR_PREC_MIN && default_prec <= MPFR_PREC_MAX,
                  "Invalid value for default precision.");
    if (m_value->_mpfr_d) {
        ::mpfr_clear(m_value);
    } else {
// NOTE: the story here is that ICC has a weird behaviour when the thread_local
// storage. Essentially, the thread-local static variable in the fma() function
// upon destruction has the _mpfr_d set to zero for some reason but the other members
// are not zeroed out. This results in the asserts below firing, and probably a memory
// leak as well as the variable is not cleared. We just disable the asserts for now.
#if !defined(PIRANHA_COMPILER_IS_INTEL)
        piranha_assert(!m_value->_mpfr_prec);
        piranha_assert(!m_value->_mpfr_sign);
        piranha_assert(!m_value->_mpfr_exp);
#endif
    }
}

namespace detail
{

// Compute gamma(a)/(gamma(b) * gamma(c)), assuming a, b and c are not negative ints. The logarithm
// of the gamma function is used internally.
inline real real_compute_3_gamma(const real &a, const real &b, const real &c, const ::mpfr_prec_t &prec)
{
    // Here we should never enter with negative ints.
    piranha_assert(a.sign() >= 0 || a.truncated() != a);
    piranha_assert(b.sign() >= 0 || b.truncated() != b);
    piranha_assert(c.sign() >= 0 || c.truncated() != c);
    const real pi = real{0, prec}.pi();
    real tmp0(0), tmp1(1);
    if (a.sign() < 0) {
        tmp0 -= (1 - a).lgamma();
        tmp1 *= pi / (a * pi).sin();
    } else {
        tmp0 += a.lgamma();
    }
    if (b.sign() < 0) {
        tmp0 += (1 - b).lgamma();
        tmp1 *= (b * pi).sin() / pi;
    } else {
        tmp0 -= b.lgamma();
    }
    if (c.sign() < 0) {
        tmp0 += (1 - c).lgamma();
        tmp1 *= (c * pi).sin() / pi;
    } else {
        tmp0 -= c.lgamma();
    }
    return tmp0.exp() * tmp1;
}
}


inline real real::binomial(const real &y) const
{
    if (unlikely(is_nan() || is_inf() || y.is_nan() || y.is_inf())) {
        piranha_throw(std::invalid_argument, "cannot compute binomial coefficient with non-finite real argument(s)");
    }
    // Work with the max precision.
    const ::mpfr_prec_t max_prec = std::max<::mpfr_prec_t>(get_prec(), y.get_prec());
    const bool neg_int_x = truncated() == (*this) && sign() < 0, neg_int_y = y.truncated() == y && y.sign() < 0,
               neg_int_x_y = ((*this) - y).truncated() == ((*this) - y) && ((*this) - y).sign() < 0;
    const unsigned mask = unsigned(neg_int_x) + (unsigned(neg_int_y) << 1u) + (unsigned(neg_int_x_y) << 2u);
    switch (mask) {
        case 0u:
            // Case 0 is the non-special one, use the default implementation.
            return detail::real_compute_3_gamma((*this) + 1, y + 1, (*this) - y + 1, max_prec);
        // NOTE: case 1 is not possible: x < 0, y > 0 implies x - y < 0 always.
        case 2u:
        case 4u:
            // These are finite numerators with infinite denominators.
            return real{0, max_prec};
        // NOTE: case 6 is not possible: x > 0, y < 0 implies x - y > 0 always.
        case 3u: {
            // 3 and 5 are the cases with 1 inf in num and 1 inf in den. Use the transformation
            // formula to make them finite.
            // NOTE: the phase here is really just a sign, but it seems tricky to compute this exactly
            // due to potential rounding errors. We are attempting to err on the safe side by using pow()
            // here.
            const auto phase = math::pow(-1, (*this) + 1) / math::pow(-1, y + 1);
            return detail::real_compute_3_gamma(-y, -(*this), (*this) - y + 1, max_prec) * phase;
        }
        case 5u: {
            const auto phase = math::pow(-1, (*this) - y + 1) / math::pow(-1, (*this) + 1);
            return detail::real_compute_3_gamma(-((*this) - y), y + 1, -(*this), max_prec) * phase;
        }
    }
    // Case 7 returns zero -> from inf / (inf * inf) it becomes a / (b * inf) after the transform.
    // NOTE: put it here so the compiler does not complain about missing return statement in the switch block.
    return real{0, max_prec};
}

inline namespace impl
{

template <typename To, typename From>
using sc_real_enabler
    = enable_if_t<conjunction<disjunction<std::is_integral<To>, is_mp_integer<To>, is_mp_rational<To>>,
                              std::is_same<From, real>>::value>;
}


template <typename To, typename From>
struct safe_cast_impl<To, From, sc_real_enabler<To, From>> {
private:
    template <typename T>
    using integral_enabler = enable_if_t<disjunction<std::is_integral<T>, is_mp_integer<T>>::value, int>;
    template <typename T>
    using rational_enabler = enable_if_t<is_mp_rational<T>::value, int>;

public:

    template <typename T = To, integral_enabler<T> = 0>
    T operator()(const real &r) const
    {
        // NOTE: the finiteness check here is repeated in the cast below, but we need it here in order
        // to make sure that we can compute the truncated r.
        if (unlikely(r.is_inf() || r.is_nan() || r.truncated() != r)) {
            piranha_throw(safe_cast_failure, "cannot convert the input real " + boost::lexical_cast<std::string>(r)
                                                 + " to the integral type '" + detail::demangle<To>()
                                                 + "', as the input real does not represent a finite integral value");
        }
        try {
            return static_cast<T>(r);
        } catch (const std::overflow_error &) {
            piranha_throw(safe_cast_failure, "cannot convert the input real " + boost::lexical_cast<std::string>(r)
                                                 + " to the integral type '" + detail::demangle<To>()
                                                 + "', as the conversion would not preserve the value");
        }
    }

    template <typename T = To, rational_enabler<T> = 0>
    T operator()(const real &r) const
    {
        try {
            return static_cast<T>(r);
        } catch (const std::overflow_error &) {
            piranha_throw(safe_cast_failure, "cannot convert the input real " + boost::lexical_cast<std::string>(r)
                                                 + " to the rational type '" + detail::demangle<To>()
                                                 + "', as the conversion would not preserve the value");
        }
    }
};

inline namespace literals
{


inline real operator"" _r(const char *s)
{
    return real(s);
}
}

inline namespace impl
{

template <typename Archive>
using real_boost_save_enabler
    = enable_if_t<conjunction<has_boost_save<Archive, ::mpfr_prec_t>, has_boost_save<Archive, std::string>,
                              has_boost_save<Archive, decltype(std::declval<const ::mpfr_t &>()->_mpfr_sign)>,
                              has_boost_save<Archive, decltype(std::declval<const ::mpfr_t &>()->_mpfr_exp)>,
                              has_boost_save<Archive, ::mp_limb_t>>::value>;

template <typename Archive>
using real_boost_load_enabler
    = enable_if_t<conjunction<has_boost_load<Archive, ::mpfr_prec_t>, has_boost_load<Archive, std::string>,
                              has_boost_load<Archive, decltype(std::declval<const ::mpfr_t &>()->_mpfr_sign)>,
                              has_boost_load<Archive, decltype(std::declval<const ::mpfr_t &>()->_mpfr_exp)>,
                              has_boost_load<Archive, ::mp_limb_t>>::value>;
}


template <typename Archive>
struct boost_save_impl<Archive, real, real_boost_save_enabler<Archive>> : boost_save_via_boost_api<Archive, real> {
};


template <typename Archive>
struct boost_load_impl<Archive, real, real_boost_load_enabler<Archive>> : boost_load_via_boost_api<Archive, real> {
};

#if defined(PIRANHA_WITH_MSGPACK)

inline namespace impl
{

// Enablers for msgpack serialization.
template <typename Stream>
using real_msgpack_pack_enabler = enable_if_t<is_detected<msgpack_pack_member_t, Stream, real>::value>;

template <typename T>
using real_msgpack_convert_enabler
    = enable_if_t<conjunction<std::is_same<real, T>, is_detected<msgpack_convert_member_t, T>>::value>;
}


template <typename Stream>
struct msgpack_pack_impl<Stream, real, real_msgpack_pack_enabler<Stream>> {

    void operator()(msgpack::packer<Stream> &p, const real &x, msgpack_format f) const
    {
        x.msgpack_pack(p, f);
    }
};


template <typename T>
struct msgpack_convert_impl<T, real_msgpack_convert_enabler<T>> {

    void operator()(T &x, const msgpack::object &o, msgpack_format f) const
    {
        x.msgpack_convert(o, f);
    }
};

#endif

inline namespace impl
{

template <typename T>
using real_zero_is_absorbing_enabler = enable_if_t<std::is_same<uncvref_t<T>, real>::value>;
}


template <typename T>
struct zero_is_absorbing<T, real_zero_is_absorbing_enabler<T>> {
    static const bool value = false;
};

template <typename T>
const bool zero_is_absorbing<T, real_zero_is_absorbing_enabler<T>>::value;
}

#endif