diff --git a/doc/changelog.rst b/doc/changelog.rst index 2247902fb..a46cc651c 100644 --- a/doc/changelog.rst +++ b/doc/changelog.rst @@ -7,6 +7,9 @@ Changelog New ~~~ +- Implement additional three-way comparison functions + for :cpp:class:`~mppp::real` + (`#243 `__). - Add the :cpp:func:`~mppp::set_ui_2exp()` and :cpp:func:`~mppp::set_si_2exp()` functions for :cpp:class:`~mppp::real`, and implement constructors from an integral multiple of a power of 2 diff --git a/doc/real.rst b/doc/real.rst index 34418f294..2eeb864cd 100644 --- a/doc/real.rst +++ b/doc/real.rst @@ -667,6 +667,14 @@ The real class :return: a reference to ``this``. + .. cpp:function:: real &sqr() + + .. versionadded:: 0.19 + + Square ``this`` in place. + + :return: a reference to ``this``. + .. cpp:function:: real &sqrt() .. cpp:function:: real &rec_sqrt() .. cpp:function:: real &sqrt1pm1() @@ -702,14 +710,6 @@ The real class :exception std\:\:invalid_argument: if the conversion between Arb and MPFR types fails because of (unlikely) overflow conditions. - .. cpp:function:: real &sqr() - - .. versionadded:: 0.19 - - Square ``this`` in place. - - :return: a reference to ``this``. - .. cpp:function:: real &sin() .. cpp:function:: real &cos() .. cpp:function:: real &tan() @@ -954,6 +954,11 @@ Types Concepts -------- +.. cpp:concept:: template mppp::cvr_real + + This concept is satisfied if the type ``T``, after the removal of reference and cv qualifiers, + is the same as :cpp:class:`mppp::real`. + .. cpp:concept:: template mppp::real_interoperable This concept is satisfied if the type ``T`` can interoperate with :cpp:class:`~mppp::real`. @@ -964,11 +969,6 @@ Concepts * a :cpp:class:`~mppp::rational`, or * :cpp:class:`~mppp::real128`. -.. cpp:concept:: template mppp::cvr_real - - This concept is satisfied if the type ``T``, after the removal of reference and cv qualifiers, - is the same as :cpp:class:`mppp::real`. - .. cpp:concept:: template mppp::real_set_args This concept is satisfied if the types in the parameter pack ``Args`` @@ -1310,6 +1310,34 @@ Arithmetic :return: *x* multiplied/divided by :math:`2^n`. +.. cpp:function:: template mppp::real &mppp::sqr(mppp::real &rop, T &&op) + + .. versionadded:: 0.19 + + Binary :cpp:class:`~mppp::real` squaring. + + This function will compute the square of *op* and store it + into *rop*. The precision of the result will be equal to the precision + of *op*. + + :param rop: the return value. + :param op: the operand. + + :return: a reference to *rop*. + +.. cpp:function:: template mppp::real mppp::sqr(T &&r) + + .. versionadded:: 0.19 + + Unary :cpp:class:`~mppp::real` squaring. + + This function will compute and return the square of *r*. + The precision of the result will be equal to the precision of *r*. + + :param r: the operand. + + :return: the square of *r*. + .. _real_comparison: Comparison @@ -1380,6 +1408,50 @@ Comparison :exception std\:\:domain_error: if at least one of the operands is NaN. +.. cpp:function:: int mppp::cmpabs(const mppp::real &a, const mppp::real &b) + + .. versionadded:: 0.20 + + Three-way comparison of absolute values. + + This function will compare *a* and *b*, returning: + + * zero if :math:`\left|a\right|=\left|b\right|`, + * a negative value if :math:`\left|a\right|<\left|b\right|`, + * a positive value if :math:`\left|a\right|>\left|b\right|`. + + If at least one NaN value is involved in the comparison, an error will be raised. + + :param a: the first operand. + :param b: the second operand. + + :return: an integral value expressing how the absolute values of *a* and *b* compare. + + :exception std\:\:domain_error: if at least one of the operands is NaN. + +.. cpp:function:: int mppp::cmp_ui_2exp(const mppp::real &a, unsigned long n, mpfr_exp_t e) +.. cpp:function:: int mppp::cmp_si_2exp(const mppp::real &a, long n, mpfr_exp_t e) + + .. versionadded:: 0.20 + + Comparison with integral multiples of powers of 2. + + This function will compare *a* to :math:`n\times 2^e`, returning: + + * zero if :math:`a=n\times 2^e`, + * a negative value if :math:`an\times 2^e`. + + If *a* is NaN, an error will be raised. + + :param a: the first operand. + :param n: the integral multiplier. + :param e: the power of 2. + + :return: an integral value expressing how *a* compares to :math:`n\times 2^e`. + + :exception std\:\:domain_error: if *a* is NaN. + .. cpp:function:: bool mppp::real_equal_to(const mppp::real &a, const mppp::real &b) Equality predicate with special handling for NaN. @@ -1634,34 +1706,6 @@ Exponentiation :return: *op1* raised to the power of *op2*. -.. cpp:function:: template mppp::real &mppp::sqr(mppp::real &rop, T &&op) - - .. versionadded:: 0.19 - - Binary :cpp:class:`~mppp::real` squaring. - - This function will compute the square of *op* and store it - into *rop*. The precision of the result will be equal to the precision - of *op*. - - :param rop: the return value. - :param op: the operand. - - :return: a reference to *rop*. - -.. cpp:function:: template mppp::real mppp::sqr(T &&r) - - .. versionadded:: 0.19 - - Unary :cpp:class:`~mppp::real` squaring. - - This function will compute and return the square of *r*. - The precision of the result will be equal to the precision of *r*. - - :param r: the operand. - - :return: the square of *r*. - .. _real_trig: Trigonometry @@ -2503,7 +2547,7 @@ Input/Output :return: a reference to *os*. - :exception unspecified: any exception thrown by :cpp:func`mppp::real::to_string()`. + :exception unspecified: any exception thrown by :cpp:func:`mppp::real::to_string()`. .. _real_operators: diff --git a/include/mp++/real.hpp b/include/mp++/real.hpp index 5685e1a72..1ab050d7c 100644 --- a/include/mp++/real.hpp +++ b/include/mp++/real.hpp @@ -1553,8 +1553,8 @@ template = 0> #endif inline real mul_2ui(T &&x, unsigned long n) { - auto mul_2ui_wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_mul_2ui(r, o, n, MPFR_RNDN); }; - return detail::mpfr_nary_op_return_impl(0, mul_2ui_wrapper, std::forward(x)); + auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_mul_2ui(r, o, n, MPFR_RNDN); }; + return detail::mpfr_nary_op_return_impl(0, wrapper, std::forward(x)); } #if defined(MPPP_HAVE_CONCEPTS) @@ -1564,8 +1564,8 @@ template = 0> #endif inline real &mul_2ui(real &rop, T &&x, unsigned long n) { - auto mul_2ui_wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_mul_2ui(r, o, n, MPFR_RNDN); }; - return detail::mpfr_nary_op_impl(0, mul_2ui_wrapper, rop, std::forward(x)); + auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_mul_2ui(r, o, n, MPFR_RNDN); }; + return detail::mpfr_nary_op_impl(0, wrapper, rop, std::forward(x)); } #if defined(MPPP_HAVE_CONCEPTS) @@ -1575,8 +1575,8 @@ template = 0> #endif inline real mul_2si(T &&x, long n) { - auto mul_2si_wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_mul_2si(r, o, n, MPFR_RNDN); }; - return detail::mpfr_nary_op_return_impl(0, mul_2si_wrapper, std::forward(x)); + auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_mul_2si(r, o, n, MPFR_RNDN); }; + return detail::mpfr_nary_op_return_impl(0, wrapper, std::forward(x)); } #if defined(MPPP_HAVE_CONCEPTS) @@ -1586,8 +1586,8 @@ template = 0> #endif inline real &mul_2si(real &rop, T &&x, long n) { - auto mul_2si_wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_mul_2si(r, o, n, MPFR_RNDN); }; - return detail::mpfr_nary_op_impl(0, mul_2si_wrapper, rop, std::forward(x)); + auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_mul_2si(r, o, n, MPFR_RNDN); }; + return detail::mpfr_nary_op_impl(0, wrapper, rop, std::forward(x)); } #if defined(MPPP_HAVE_CONCEPTS) @@ -1597,8 +1597,8 @@ template = 0> #endif inline real div_2ui(T &&x, unsigned long n) { - auto div_2ui_wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_div_2ui(r, o, n, MPFR_RNDN); }; - return detail::mpfr_nary_op_return_impl(0, div_2ui_wrapper, std::forward(x)); + auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_div_2ui(r, o, n, MPFR_RNDN); }; + return detail::mpfr_nary_op_return_impl(0, wrapper, std::forward(x)); } #if defined(MPPP_HAVE_CONCEPTS) @@ -1608,8 +1608,8 @@ template = 0> #endif inline real &div_2ui(real &rop, T &&x, unsigned long n) { - auto div_2ui_wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_div_2ui(r, o, n, MPFR_RNDN); }; - return detail::mpfr_nary_op_impl(0, div_2ui_wrapper, rop, std::forward(x)); + auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_div_2ui(r, o, n, MPFR_RNDN); }; + return detail::mpfr_nary_op_impl(0, wrapper, rop, std::forward(x)); } #if defined(MPPP_HAVE_CONCEPTS) @@ -1619,8 +1619,8 @@ template = 0> #endif inline real div_2si(T &&x, long n) { - auto div_2si_wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_div_2si(r, o, n, MPFR_RNDN); }; - return detail::mpfr_nary_op_return_impl(0, div_2si_wrapper, std::forward(x)); + auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_div_2si(r, o, n, MPFR_RNDN); }; + return detail::mpfr_nary_op_return_impl(0, wrapper, std::forward(x)); } #if defined(MPPP_HAVE_CONCEPTS) @@ -1630,8 +1630,8 @@ template = 0> #endif inline real &div_2si(real &rop, T &&x, long n) { - auto div_2si_wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_div_2si(r, o, n, MPFR_RNDN); }; - return detail::mpfr_nary_op_impl(0, div_2si_wrapper, rop, std::forward(x)); + auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_div_2si(r, o, n, MPFR_RNDN); }; + return detail::mpfr_nary_op_impl(0, wrapper, rop, std::forward(x)); } // Detect NaN. @@ -1691,6 +1691,13 @@ inline bool signbit(const real &r) // Comparison. MPPP_DLL_PUBLIC int cmp(const real &, const real &); +// Comparison of absolute values. +MPPP_DLL_PUBLIC int cmpabs(const real &, const real &); + +// Comparison with integral multiples of powers of 2. +MPPP_DLL_PUBLIC int cmp_ui_2exp(const real &, unsigned long, ::mpfr_exp_t); +MPPP_DLL_PUBLIC int cmp_si_2exp(const real &, long, ::mpfr_exp_t); + // Equality predicate with special NaN handling. MPPP_DLL_PUBLIC bool real_equal_to(const real &, const real &); @@ -1892,6 +1899,15 @@ inline real dispatch_real_pow(T &&a, const U &n) } } +// Special casing for bool +template ::value, int> = 0> +inline real dispatch_real_pow(T &&a, const bool &n) +{ + auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_pow_ui(r, o, static_cast(n), MPFR_RNDN); }; + + return mpfr_nary_op_return_impl(real_deduce_precision(n), wrapper, std::forward(a)); +} + // real-signed integral. template , is_cpp_signed_integral>::value, int> = 0> inline real dispatch_real_pow(T &&a, const U &n) @@ -1951,6 +1967,15 @@ inline real dispatch_real_pow(const T &n, U &&a) } } +// Special casing for bool. +template ::value, int> = 0> +inline real dispatch_real_pow(const bool &n, T &&a) +{ + auto wrapper = [n](::mpfr_t r, const ::mpfr_t o) { ::mpfr_ui_pow(r, static_cast(n), o, MPFR_RNDN); }; + + return mpfr_nary_op_return_impl(real_deduce_precision(n), wrapper, std::forward(a)); +} + } // namespace detail // Binary exponentiation. diff --git a/src/real.cpp b/src/real.cpp index 0b37833eb..f2ef57b1b 100644 --- a/src/real.cpp +++ b/src/real.cpp @@ -2566,6 +2566,41 @@ int cmp(const real &a, const real &b) return retval; } +// Three-way comparison of absolute values. +int cmpabs(const real &a, const real &b) +{ + ::mpfr_clear_erangeflag(); + auto retval = ::mpfr_cmpabs(a.get_mpfr_t(), b.get_mpfr_t()); + if (mppp_unlikely(::mpfr_erangeflag_p())) { + ::mpfr_clear_erangeflag(); + throw std::domain_error("Cannot compare the absolute values of two reals if at least one of them is NaN"); + } + return retval; +} + +// Comparison with n*2**e. +int cmp_ui_2exp(const real &x, unsigned long n, ::mpfr_exp_t e) +{ + ::mpfr_clear_erangeflag(); + auto retval = ::mpfr_cmp_ui_2exp(x.get_mpfr_t(), n, e); + if (mppp_unlikely(::mpfr_erangeflag_p())) { + ::mpfr_clear_erangeflag(); + throw std::domain_error("Cannot compare a real NaN to an integral multiple of a power of 2"); + } + return retval; +} + +int cmp_si_2exp(const real &x, long n, ::mpfr_exp_t e) +{ + ::mpfr_clear_erangeflag(); + auto retval = ::mpfr_cmp_si_2exp(x.get_mpfr_t(), n, e); + if (mppp_unlikely(::mpfr_erangeflag_p())) { + ::mpfr_clear_erangeflag(); + throw std::domain_error("Cannot compare a real NaN to an integral multiple of a power of 2"); + } + return retval; +} + // Equality predicate with special NaN handling. bool real_equal_to(const real &a, const real &b) { diff --git a/test/real_arith.cpp b/test/real_arith.cpp index c71013d42..50d814b37 100644 --- a/test/real_arith.cpp +++ b/test/real_arith.cpp @@ -6,6 +6,7 @@ // Public License v. 2.0. If a copy of the MPL was not distributed // with this file, You can obtain one at http://mozilla.org/MPL/2.0/. +#include #include #include @@ -295,6 +296,7 @@ TEST_CASE("real fma") { real r1, r2, r3, r4; fma_wrap(r1, r2, r3, r4); + REQUIRE(std::is_same::value); REQUIRE(r1.zero_p()); REQUIRE(r1.get_prec() == r3.get_prec()); fma_wrap(r1, real{2, 12}, real{3, 7}, real{14, 128}); @@ -309,6 +311,7 @@ TEST_CASE("real fma") REQUIRE(::mpfr_equal_p(r1.get_mpfr_t(), real{44}.get_mpfr_t())); REQUIRE(r1.get_prec() == 128); r1 = fma_wrap(real{14, 128}, real{3, 7}, real{2, 12}); + REQUIRE(std::is_same::value); REQUIRE(::mpfr_equal_p(r1.get_mpfr_t(), real{44}.get_mpfr_t())); REQUIRE(r1.get_prec() == 128); r1 = fma_wrap(static_cast(real{14, 128}), real{3, 7}, real{2, 12}); @@ -348,6 +351,7 @@ TEST_CASE("real fms") { real r1, r2, r3, r4; fms_wrap(r1, r2, r3, r4); + REQUIRE(std::is_same::value); REQUIRE(r1.zero_p()); REQUIRE(r1.get_prec() == r3.get_prec()); fms_wrap(r1, real{2, 12}, real{3, 7}, real{14, 128}); @@ -362,6 +366,7 @@ TEST_CASE("real fms") REQUIRE(::mpfr_equal_p(r1.get_mpfr_t(), real{40}.get_mpfr_t())); REQUIRE(r1.get_prec() == 128); r1 = fms_wrap(real{14, 128}, real{3, 7}, real{2, 12}); + REQUIRE(std::is_same::value); REQUIRE(::mpfr_equal_p(r1.get_mpfr_t(), real{40}.get_mpfr_t())); REQUIRE(r1.get_prec() == 128); r1 = fms_wrap(static_cast(real{14, 128}), real{3, 7}, real{2, 12}); diff --git a/test/real_cmp.cpp b/test/real_cmp.cpp index fc4b1cfc7..eb9eac6f9 100644 --- a/test/real_cmp.cpp +++ b/test/real_cmp.cpp @@ -164,3 +164,54 @@ TEST_CASE("real is_one") REQUIRE(!real{"-nan", 5}.is_one()); REQUIRE(!::mpfr_erangeflag_p()); } + +TEST_CASE("real cmpabs") +{ + REQUIRE(cmpabs(real{}, real{}) == 0); + REQUIRE(cmpabs(real{1}, real{1}) == 0); + REQUIRE(cmpabs(real{1}, real{0}) > 0); + REQUIRE(cmpabs(real{-1}, real{0}) > 0); + REQUIRE(cmpabs(real{0}, real{1}) < 0); + REQUIRE(cmpabs(real{0}, real{-1}) < 0); + REQUIRE(cmpabs(real{"inf", 64}, real{45}) > 0); + REQUIRE(cmpabs(-real{"inf", 64}, real{45}) > 0); + REQUIRE(cmpabs(real{45}, real{"inf", 64}) < 0); + REQUIRE(cmpabs(real{45}, -real{"inf", 64}) < 0); + REQUIRE(cmpabs(-real{"inf", 64}, -real{"inf", 4}) == 0); + REQUIRE(cmpabs(real{"inf", 64}, -real{"inf", 4}) == 0); + REQUIRE_THROWS_PREDICATE(cmpabs(real{"nan", 5}, real{6}), std::domain_error, [](const std::domain_error &ex) { + return ex.what() + == std::string("Cannot compare the absolute values of two reals if at least one of them is NaN"); + }); + REQUIRE_THROWS_PREDICATE(cmpabs(real{6}, real{"nan", 5}), std::domain_error, [](const std::domain_error &ex) { + return ex.what() + == std::string("Cannot compare the absolute values of two reals if at least one of them is NaN"); + }); + REQUIRE_THROWS_PREDICATE( + cmpabs(real{"nan", 5}, real{"nan", 5}), std::domain_error, [](const std::domain_error &ex) { + return ex.what() + == std::string("Cannot compare the absolute values of two reals if at least one of them is NaN"); + }); +} + +TEST_CASE("real cmp_2_exp") +{ + REQUIRE(cmp_ui_2exp(real{}, 0, 0) == 0); + REQUIRE(cmp_si_2exp(real{}, 0, 0) == 0); + + REQUIRE(cmp_ui_2exp(real{16}, 4, 2) == 0); + REQUIRE(cmp_si_2exp(real{16}, 4, 2) == 0); + + REQUIRE(cmp_ui_2exp(real{2}, 0, 0) > 0); + REQUIRE(cmp_si_2exp(real{2}, 0, 0) > 0); + + REQUIRE(cmp_ui_2exp(real{-2}, 4, 2) < 0); + REQUIRE(cmp_si_2exp(real{-32}, -4, 2) < 0); + + REQUIRE_THROWS_PREDICATE(cmp_ui_2exp(real{"nan", 5}, 4, 2), std::domain_error, [](const std::domain_error &ex) { + return ex.what() == std::string("Cannot compare a real NaN to an integral multiple of a power of 2"); + }); + REQUIRE_THROWS_PREDICATE(cmp_si_2exp(real{"nan", 5}, 4, 2), std::domain_error, [](const std::domain_error &ex) { + return ex.what() == std::string("Cannot compare a real NaN to an integral multiple of a power of 2"); + }); +} diff --git a/test/real_mul_div_2.cpp b/test/real_mul_div_2.cpp index e0cc929b6..2e841ab10 100644 --- a/test/real_mul_div_2.cpp +++ b/test/real_mul_div_2.cpp @@ -19,6 +19,7 @@ TEST_CASE("real mul_2ui") { // The return form. REQUIRE(mul_2ui(2_r128, 0) == 2); + REQUIRE(std::is_same::value); REQUIRE(mul_2ui(2_r128, 1) == 4); REQUIRE(mul_2ui(2_r128, 2) == 8); REQUIRE(mul_2ui(2_r128, 2).get_prec() == 128); @@ -61,6 +62,7 @@ TEST_CASE("real mul_2si") { // The return form. REQUIRE(mul_2si(2_r128, 0) == 2); + REQUIRE(std::is_same::value); REQUIRE(mul_2si(2_r128, -1) == 1); REQUIRE(mul_2si(2_r128, -2) == 2_r128 / 4); REQUIRE(mul_2si(2_r128, 2).get_prec() == 128); @@ -103,6 +105,7 @@ TEST_CASE("real div_2ui") { // The return form. REQUIRE(div_2ui(2_r128, 0) == 2); + REQUIRE(std::is_same::value); REQUIRE(div_2ui(2_r128, 1) == 1); REQUIRE(div_2ui(2_r128, 2) == 2_r128 / 4); REQUIRE(div_2ui(2_r128, 2).get_prec() == 128); @@ -145,6 +148,7 @@ TEST_CASE("real div_2si") { // The return form. REQUIRE(div_2si(2_r128, 0) == 2); + REQUIRE(std::is_same::value); REQUIRE(div_2si(2_r128, 1) == 1); REQUIRE(div_2si(2_r128, -1) == 4); REQUIRE(div_2si(2_r128, -2) == 8); diff --git a/test/real_pow.cpp b/test/real_pow.cpp index d728893ea..05e8648ce 100644 --- a/test/real_pow.cpp +++ b/test/real_pow.cpp @@ -159,6 +159,21 @@ TEST_CASE("real pow") // Ensure that x**(1/3) is almost identical to cbrt(1.1). REQUIRE(abs(pow(1.1_r512, 1_q1 / 3) - cbrt(1.1_r512)) <= pow(2_r512, -500)); + + // Special casing for bool. + REQUIRE(mppp::pow(real{123}, false) == 1); + REQUIRE(mppp::pow(real{123}, false).get_prec() + == std::max(detail::real_deduce_precision(123), detail::real_deduce_precision(false))); + REQUIRE(mppp::pow(real{123}, true) == 123); + REQUIRE(mppp::pow(real{123}, true).get_prec() + == std::max(detail::real_deduce_precision(123), detail::real_deduce_precision(false))); + + REQUIRE(mppp::pow(false, real{123}) == 0); + REQUIRE(mppp::pow(false, real{123}).get_prec() + == std::max(detail::real_deduce_precision(123), detail::real_deduce_precision(false))); + REQUIRE(mppp::pow(true, real{123}) == 1); + REQUIRE(mppp::pow(true, real{123}).get_prec() + == std::max(detail::real_deduce_precision(123), detail::real_deduce_precision(false))); } TEST_CASE("real sqr") diff --git a/test/real_roots.cpp b/test/real_roots.cpp index 2ff95ce24..3404072be 100644 --- a/test/real_roots.cpp +++ b/test/real_roots.cpp @@ -21,12 +21,15 @@ TEST_CASE("real sqrt") { real r0{0}; r0.sqrt(); + REQUIRE(std::is_same::value); REQUIRE(r0.get_prec() == detail::real_deduce_precision(0)); REQUIRE(r0.zero_p()); real rop; REQUIRE(sqrt(rop, r0).zero_p()); + REQUIRE(std::is_same::value); REQUIRE(rop.get_prec() == detail::real_deduce_precision(0)); REQUIRE(sqrt(r0).zero_p()); + REQUIRE(std::is_same::value); REQUIRE(sqrt(std::move(r0)).zero_p()); REQUIRE(!r0.get_mpfr_t()->_mpfr_d); r0 = real{16, 128};