diff --git a/README.md b/README.md
index b1051d73..f76b5881 100644
--- a/README.md
+++ b/README.md
@@ -26,7 +26,7 @@ requiring extended range and precision.
 
 Multiprecision consists of a generic interface to the mathematics
 of large numbers as well as a selection of big number backends.
-This includes interfaces to GMP, MPFR, MPIR and TomMath
+These include interfaces to GMP, MPFR, MPIR and TomMath
 and also Multiprecision's own collection of Boost-licensed,
 header-only backends for integers, rationals, floats and complex-floats.
 
diff --git a/include/boost/multiprecision/cpp_dec_float.hpp b/include/boost/multiprecision/cpp_dec_float.hpp
index c9b958c0..df1a173b 100644
--- a/include/boost/multiprecision/cpp_dec_float.hpp
+++ b/include/boost/multiprecision/cpp_dec_float.hpp
@@ -427,7 +427,6 @@ class cpp_dec_float
    cpp_dec_float& div_unsigned_long_long(const unsigned long long n);
 
    // Elementary primitives.
-   cpp_dec_float& calculate_inv();
    cpp_dec_float& calculate_sqrt();
 
    void negate()
@@ -568,6 +567,9 @@ class cpp_dec_float
    static bool data_elem_is_non_nine_predicate(const std::uint32_t& d) { return (d != static_cast<std::uint32_t>(cpp_dec_float::cpp_dec_float_elem_mask - 1)); }
    static bool char_is_nonzero_predicate(const char& c) { return (c != static_cast<char>('0')); }
 
+   // Inversion.
+   cpp_dec_float& calculate_inv();
+
    void from_unsigned_long_long(const unsigned long long u);
 
    template <typename InputIteratorTypeLeft,
@@ -1218,20 +1220,16 @@ cpp_dec_float<Digits10, ExponentType, Allocator>& cpp_dec_float<Digits10, Expone
 template <unsigned Digits10, class ExponentType, class Allocator>
 cpp_dec_float<Digits10, ExponentType, Allocator>& cpp_dec_float<Digits10, ExponentType, Allocator>::calculate_inv()
 {
-   // Compute the inverse of *this.
-   const bool b_neg = neg;
-
-   neg = false;
-
-   // Handle special cases like zero, inf and NaN.
+   // Handle the special case of zero.
    if (iszero())
    {
       *this = inf();
-      if (b_neg)
-         negate();
+
       return *this;
    }
 
+   // Handle the special cases of inf and NaN.
+
    if ((isnan)())
    {
       return *this;
@@ -1242,14 +1240,12 @@ cpp_dec_float<Digits10, ExponentType, Allocator>& cpp_dec_float<Digits10, Expone
       return *this = zero();
    }
 
-   if (isone())
-   {
-      if (b_neg)
-         negate();
-      return *this;
-   }
+   // Compute the inverse of *this.
+   const bool b_neg = neg;
 
-   // Save the original *this.
+   neg = false;
+
+   // Save the original (absolute value of) *this.
    cpp_dec_float<Digits10, ExponentType, Allocator> x(*this);
 
    // Generate the initial estimate using division.
diff --git a/include/boost/multiprecision/cpp_double_fp.hpp b/include/boost/multiprecision/cpp_double_fp.hpp
index 4e8245f1..9049e0bc 100644
--- a/include/boost/multiprecision/cpp_double_fp.hpp
+++ b/include/boost/multiprecision/cpp_double_fp.hpp
@@ -163,9 +163,11 @@ struct number_category<backends::cpp_double_fp_backend<FloatingPointType>> : pub
 
 namespace backends {
 
-// A cpp_double_fp_backend is represented by an unevaluated sum of two floating-point
-// units (say a0 and a1) which satisfy |a1| <= (1 / 2) * ulp(a0).
+// A cpp_double_fp_backend is represented by an unevaluated sum of two
+// floating-point units, a0 and a1, which satisfy |a1| <= (1 / 2) * ulp(a0).
 // The type of the floating-point constituents should adhere to IEEE754.
+// This class has been tested with floats having single-precision (4 byte),
+// double-precision (8 byte) and quad precision (16 byte, such as GCC's __float128).
 
 template <typename FloatingPointType>
 class cpp_double_fp_backend
@@ -178,10 +180,10 @@ class cpp_double_fp_backend
          cpp_df_qf_detail::is_floating_point<float_type>::value
       && bool
          {
-               (cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits ==  24)
-            || (cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits ==  53)
-            || (cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits ==  64)
-            || (cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits == 113)
+               ((cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits ==  24) && std::numeric_limits<float_type>::is_specialized && std::numeric_limits<float_type>::is_iec559)
+            || ((cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits ==  53) && std::numeric_limits<float_type>::is_specialized && std::numeric_limits<float_type>::is_iec559)
+            || ((cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits ==  64) && std::numeric_limits<float_type>::is_specialized && std::numeric_limits<float_type>::is_iec559)
+            ||  (cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits == 113)
          }, "Error: float_type does not fulfil the backend requirements of cpp_double_fp_backend"
    );
 
@@ -514,12 +516,10 @@ class cpp_double_fp_backend
          }
       }
 
-      const float_type xlo { data.second };
+      const rep_type thi_tlo { arithmetic::two_sum(data.second, v.data.second) };
 
       data = arithmetic::two_sum(data.first, v.data.first);
 
-      const rep_type thi_tlo { arithmetic::two_sum(xlo, v.data.second) };
-
       data = arithmetic::two_hilo_sum(data.first, data.second + thi_tlo.first);
 
       data = arithmetic::two_hilo_sum(data.first, thi_tlo.second + data.second);
@@ -578,12 +578,10 @@ class cpp_double_fp_backend
          return *this;
       }
 
-      const float_type xlo { data.second };
+      const rep_type thi_tlo { arithmetic::two_diff(data.second, v.data.second) };
 
       data = arithmetic::two_diff(data.first, v.data.first);
 
-      const rep_type thi_tlo { arithmetic::two_diff(xlo, v.data.second) };
-
       data = arithmetic::two_hilo_sum(data.first, data.second + thi_tlo.first);
 
       data = arithmetic::two_hilo_sum(data.first, thi_tlo.second + data.second);
@@ -993,7 +991,8 @@ class cpp_double_fp_backend
       constexpr cpp_double_fp_backend
          my_value_eps_constexpr
          {
-            cpp_df_qf_detail::ccmath::unsafe::ldexp(float_type { 1 }, int { 3 - my_digits })
+              cpp_double_fp_backend(cpp_df_qf_detail::ccmath::numeric_limits<float_type>::epsilon())
+            * cpp_double_fp_backend(cpp_df_qf_detail::ccmath::numeric_limits<float_type>::epsilon())
          };
 
       static_assert
@@ -1613,13 +1612,15 @@ template <typename FloatingPointType,
           typename ::std::enable_if<(cpp_df_qf_detail::is_floating_point<FloatingPointType>::value && ((cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits10 * 2) < 16))>::type const*>
 constexpr auto eval_exp(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void
 {
-   const int fpc { eval_fpclassify(x) };
-
    using double_float_type = cpp_double_fp_backend<FloatingPointType>;
 
+   constexpr double_float_type one { 1 };
+
+   const int fpc { eval_fpclassify(x) };
+
    if (fpc == FP_ZERO)
    {
-      result = double_float_type { 1.0F };
+      result = one;
    }
    else if (fpc != FP_NORMAL)
    {
@@ -1751,7 +1752,7 @@ constexpr auto eval_exp(cpp_double_fp_backend<FloatingPointType>& result, const
 
          if (b_neg)
          {
-            result = double_float_type(1U) / result;
+            result = one / result;
          }
       }
    }
@@ -1761,13 +1762,15 @@ template <typename FloatingPointType,
           typename ::std::enable_if<(cpp_df_qf_detail::is_floating_point<FloatingPointType>::value && (((cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits10 * 2) >= 16) && ((cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits10 * 2) <= 36)))>::type const*>
 constexpr auto eval_exp(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void
 {
-   const int fpc { eval_fpclassify(x) };
-
    using double_float_type = cpp_double_fp_backend<FloatingPointType>;
 
+   constexpr double_float_type one { 1 };
+
+   const int fpc { eval_fpclassify(x) };
+
    if (fpc == FP_ZERO)
    {
-      result = double_float_type(1);
+      result = one;
    }
    else if (fpc != FP_NORMAL)
    {
@@ -1901,7 +1904,7 @@ constexpr auto eval_exp(cpp_double_fp_backend<FloatingPointType>& result, const
 
          if (b_neg)
          {
-            result = double_float_type(1U) / result;
+            result = one / result;
          }
       }
    }
@@ -1911,13 +1914,15 @@ template <typename FloatingPointType,
           typename ::std::enable_if<(cpp_df_qf_detail::is_floating_point<FloatingPointType>::value && ((cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits10 * 2) > 36))>::type const*>
 constexpr auto eval_exp(cpp_double_fp_backend<FloatingPointType>& result, const cpp_double_fp_backend<FloatingPointType>& x) -> void
 {
-   const int fpc { eval_fpclassify(x) };
-
    using double_float_type = cpp_double_fp_backend<FloatingPointType>;
 
+   constexpr double_float_type one { 1 };
+
+   const int fpc { eval_fpclassify(x) };
+
    if (fpc == FP_ZERO)
    {
-      result = double_float_type(1);
+      result = one;
    }
    else if (fpc != FP_NORMAL)
    {
@@ -2012,16 +2017,17 @@ constexpr auto eval_exp(cpp_double_fp_backend<FloatingPointType>& result, const
          // Series expansion of hypergeometric_0f0(; ; x).
          // For this high(er) digit count, a scaled argument with subsequent
          // Taylor series expansion is actually more precise than Pade approximation.
-         for (unsigned n = 2U; n < 64U; ++n)
+         for (unsigned n { 2U }; n < 64U; ++n)
          {
             x_pow_n_div_n_fact *= xh;
+
             x_pow_n_div_n_fact /= typename double_float_type::float_type(n);
 
             int n_tol { };
 
             eval_frexp(dummy, x_pow_n_div_n_fact, &n_tol);
 
-            if ((n > 4U) && (n_tol < -(double_float_type::my_digits - 4)))
+            if ((n > 4U) && (n_tol < -(double_float_type::my_digits - 2)))
             {
                break;
             }
@@ -2051,7 +2057,7 @@ constexpr auto eval_exp(cpp_double_fp_backend<FloatingPointType>& result, const
 
          if (b_neg)
          {
-            result = double_float_type(1U) / result;
+            result = one / result;
          }
       }
    }
diff --git a/include/boost/multiprecision/detail/default_ops.hpp b/include/boost/multiprecision/detail/default_ops.hpp
index d66fe4a1..6e1aeb35 100644
--- a/include/boost/multiprecision/detail/default_ops.hpp
+++ b/include/boost/multiprecision/detail/default_ops.hpp
@@ -499,7 +499,7 @@ inline BOOST_MP_CXX14_CONSTEXPR typename std::enable_if< !std::is_same<T, U>::va
    T t;
    t = number<T>::canonical_value(u);
    return t;
-}
+} // LCOV_EXCL_LINE
 template <class T>
 inline BOOST_MP_CXX14_CONSTEXPR const T& make_T(const T& t)
 {
@@ -1803,6 +1803,9 @@ BOOST_MP_CXX14_CONSTEXPR void eval_karatsuba_sqrt(Backend& result, const Backend
    result = s;
 }
 
+// LCOV_EXCL_START
+// This is known tested in the test-file test_various_edges.cpp.
+
 template <class B>
 void BOOST_MP_CXX14_CONSTEXPR eval_integer_sqrt_bitwise(B& s, B& r, const B& x)
 {
@@ -1861,6 +1864,8 @@ void BOOST_MP_CXX14_CONSTEXPR eval_integer_sqrt_bitwise(B& s, B& r, const B& x)
    } while (g >= 0);
 }
 
+// LCOV_EXCL_STOP
+
 template <class Backend>
 BOOST_MP_CXX14_CONSTEXPR void eval_integer_sqrt(Backend& result, Backend& r, const Backend& x)
 {
@@ -3867,7 +3872,7 @@ conj(const number<Backend, et_off>& arg)
    using default_ops::eval_conj;
    eval_conj(result.backend(), arg.backend());
    return result;
-}
+} // LCOV_EXCL_LINE
 
 template <class tag, class A1, class A2, class A3, class A4>
 inline BOOST_MP_CXX14_CONSTEXPR detail::expression<
diff --git a/test/test_exp.cpp b/test/test_exp.cpp
index 776b6e39..6cc7f1d1 100644
--- a/test/test_exp.cpp
+++ b/test/test_exp.cpp
@@ -14,8 +14,9 @@
 #pragma warning(disable : 4127)
 #endif
 
+#include <test.hpp>
+
 #include <boost/detail/lightweight_test.hpp>
-#include "test.hpp"
 
 #include <array>
 #include <ctime>
diff --git a/test/test_log.cpp b/test/test_log.cpp
index 3860dbdf..87ada463 100644
--- a/test/test_log.cpp
+++ b/test/test_log.cpp
@@ -13,9 +13,11 @@
 #define _SCL_SECURE_NO_WARNINGS
 #endif
 
+#include <test.hpp>
+
 #include <boost/detail/lightweight_test.hpp>
+
 #include <array>
-#include "test.hpp"
 
 #if !defined(TEST_MPF_50) && !defined(TEST_MPF) && !defined(TEST_BACKEND) && !defined(TEST_MPZ) && !defined(TEST_CPP_DEC_FLOAT) && !defined(TEST_MPFR) && !defined(TEST_MPFR_50) && !defined(TEST_MPFI_50) && !defined(TEST_FLOAT128) && !defined(TEST_CPP_BIN_FLOAT) && !defined(TEST_CPP_DOUBLE_FLOAT)
 #define TEST_MPF_50
diff --git a/test/test_various_edges.cpp b/test/test_various_edges.cpp
index 9db4378c..03363efc 100644
--- a/test/test_various_edges.cpp
+++ b/test/test_various_edges.cpp
@@ -28,10 +28,7 @@
 #include <string>
 #include <typeinfo>
 
-#if defined(__clang__)
-#  pragma clang diagnostic push
-#  pragma clang diagnostic ignored "-Wfloat-equal"
-#elif defined(__GNUC__)
+#if (defined(__clang__) || defined(__GNUC__))
 #  pragma GCC diagnostic push
 #  pragma GCC diagnostic ignored "-Wfloat-equal"
 #endif
@@ -97,13 +94,13 @@ namespace local
     {
       delta = fabs(a - b); // LCOV_EXCL_LINE
 
-      result_is_ok = (delta < tol); // LCOV_EXCL_LINE
+      result_is_ok = (delta <= tol); // LCOV_EXCL_LINE
     }
     else
     {
       delta = fabs(1 - (a / b));
 
-      result_is_ok = (delta < tol);
+      result_is_ok = (delta <= tol);
     }
 
     return result_is_ok;
@@ -501,6 +498,103 @@ namespace local
     return result_is_ok;
   }
 
+  template<typename FloatType>
+  auto test_fdim_edge() -> bool
+  {
+    using float_type = FloatType;
+
+    std::mt19937_64 gen { time_point<typename std::mt19937_64::result_type>() };
+
+    std::uniform_real_distribution<float>
+      dist
+      (
+        static_cast<float>(-125.5L),
+        static_cast<float>(+125.5L)
+      );
+
+    auto result_is_ok = true;
+
+    for(auto i = static_cast<unsigned>(UINT8_C(0)); i < static_cast<unsigned>(UINT8_C(4)); ++i)
+    {
+      static_cast<void>(i);
+
+      const float_type val_nrm = ::my_one<float_type>() * static_cast<float_type>(dist(gen));
+      const float_type val_nan = std::numeric_limits<float_type>::quiet_NaN() * static_cast<float_type>(dist(gen));
+      const float_type val_inf = std::numeric_limits<float_type>::infinity() * static_cast<float_type>(dist(gen));
+
+      const double flt_nrm = static_cast<double>(val_nrm);
+      const double flt_nan = std::numeric_limits<double>::quiet_NaN() * static_cast<double>(dist(gen));
+      const double flt_inf = std::numeric_limits<double>::infinity() * static_cast<double>(dist(gen));
+
+      using std::fdim;
+      using std::fpclassify;
+
+      float_type result_fdim { fdim(val_nrm, val_nan) };
+      bool result_fdim_is_ok = (fpclassify(result_fdim) == FP_ZERO);
+      BOOST_TEST(result_fdim_is_ok);
+      result_is_ok = (result_fdim_is_ok && result_is_ok);
+
+      result_fdim = fdim(val_nrm, val_inf);
+      result_fdim_is_ok = (fpclassify(result_fdim) == FP_ZERO);
+      BOOST_TEST(result_fdim_is_ok);
+      result_is_ok = (result_fdim_is_ok && result_is_ok);
+
+      result_fdim = fdim(val_nan, val_nrm);
+      result_fdim_is_ok = (fpclassify(result_fdim) == FP_ZERO);
+      BOOST_TEST(result_fdim_is_ok);
+      result_is_ok = (result_fdim_is_ok && result_is_ok);
+
+      result_fdim = fdim(val_inf, val_nrm);
+      result_fdim_is_ok = (fpclassify(result_fdim) == FP_INFINITE);
+      BOOST_TEST(result_fdim_is_ok);
+      result_is_ok = (result_fdim_is_ok && result_is_ok);
+
+      // MP-type argument-a mixed with built-in argument-b.
+      result_fdim = fdim(val_nrm, flt_nan);
+      result_fdim_is_ok = (fpclassify(result_fdim) == FP_ZERO);
+      BOOST_TEST(result_fdim_is_ok);
+      result_is_ok = (result_fdim_is_ok && result_is_ok);
+
+      result_fdim = fdim(val_nrm, flt_inf);
+      result_fdim_is_ok = (fpclassify(result_fdim) == FP_ZERO);
+      BOOST_TEST(result_fdim_is_ok);
+      result_is_ok = (result_fdim_is_ok && result_is_ok);
+
+      result_fdim = fdim(val_nan, flt_nrm);
+      result_fdim_is_ok = (fpclassify(result_fdim) == FP_ZERO);
+      BOOST_TEST(result_fdim_is_ok);
+      result_is_ok = (result_fdim_is_ok && result_is_ok);
+
+      result_fdim = fdim(val_inf, flt_nrm);
+      result_fdim_is_ok = (fpclassify(result_fdim) == FP_INFINITE);
+      BOOST_TEST(result_fdim_is_ok);
+      result_is_ok = (result_fdim_is_ok && result_is_ok);
+
+      // Built-in argument-a mixed with MP-type argument-b.
+      result_fdim = fdim(flt_nrm, val_nan);
+      result_fdim_is_ok = (fpclassify(result_fdim) == FP_ZERO);
+      BOOST_TEST(result_fdim_is_ok);
+      result_is_ok = (result_fdim_is_ok && result_is_ok);
+
+      result_fdim = fdim(flt_nrm, val_inf);
+      result_fdim_is_ok = (fpclassify(result_fdim) == FP_ZERO);
+      BOOST_TEST(result_fdim_is_ok);
+      result_is_ok = (result_fdim_is_ok && result_is_ok);
+
+      result_fdim = fdim(flt_nan, val_nrm);
+      result_fdim_is_ok = (fpclassify(result_fdim) == FP_ZERO);
+      BOOST_TEST(result_fdim_is_ok);
+      result_is_ok = (result_fdim_is_ok && result_is_ok);
+
+      result_fdim = fdim(flt_inf, val_nrm);
+      result_fdim_is_ok = (fpclassify(result_fdim) == FP_INFINITE);
+      BOOST_TEST(result_fdim_is_ok);
+      result_is_ok = (result_fdim_is_ok && result_is_ok);
+    }
+
+    return result_is_ok;
+  }
+
   template<typename FloatType>
   auto test_sqrt_edge() -> bool
   {
@@ -1523,10 +1617,72 @@ namespace local
 
     return result_is_ok;
   }
+
+  auto test_sqrt_integral_and_constexpr() -> void
+  {
+    {
+      // Select some pseudo-random integers
+
+      // Table[N[Exp[Pi  EulerGamma*m], 36], {m, 1, 301, 100}]
+      // 6.13111418242260482895474318171556018 632831453348
+      // 3.47920348866883608793309754527486715 109524826009*10^79
+      // 1.97433232450131483625758030703105499 465632954349*10^158
+      // 1.12036796360599148200087404494586705 923184966483*10^237
+
+      // N[Sqrt[613111418242260482895474318171556018], 50]
+      // N[Sqrt[347920348866883608793309754527486715], 50]
+      // N[Sqrt[197433232450131483625758030703105499], 50]
+      // N[Sqrt[112036796360599148200087404494586705], 50]
+
+      // 7.8301431547722068635982391629827009581580426009374*10^17
+      // 5.8984773362867438315481181521503243715373594761667*10^17
+      // 4.4433459515339505424085647442953457083310383400346*10^17
+      // 3.3471898117764273004961506045440891118909367547368*10^17
+
+      // 783014315477220686
+      // 589847733628674383
+      // 444334595153395054
+      // 334718981177642730
+
+      constexpr boost::multiprecision::uint128_t un0 = boost::multiprecision::uint128_t(UINT64_C(613111418242260482)) * boost::multiprecision::uint128_t(UINT64_C(1000000000000000000)) + boost::multiprecision::uint128_t(UINT64_C(895474318171556018));
+      constexpr boost::multiprecision::uint128_t un1 = boost::multiprecision::uint128_t(UINT64_C(347920348866883608)) * boost::multiprecision::uint128_t(UINT64_C(1000000000000000000)) + boost::multiprecision::uint128_t(UINT64_C(793309754527486715));
+      constexpr boost::multiprecision::uint128_t un2 = boost::multiprecision::uint128_t(UINT64_C(197433232450131483)) * boost::multiprecision::uint128_t(UINT64_C(1000000000000000000)) + boost::multiprecision::uint128_t(UINT64_C(625758030703105499));
+      constexpr boost::multiprecision::uint128_t un3 = boost::multiprecision::uint128_t(UINT64_C(112036796360599148)) * boost::multiprecision::uint128_t(UINT64_C(1000000000000000000)) + boost::multiprecision::uint128_t(UINT64_C(200087404494586705));
+
+      constexpr boost::multiprecision::uint128_t sqrt_un0(sqrt(un0));
+      constexpr boost::multiprecision::uint128_t sqrt_un1(sqrt(un1));
+      constexpr boost::multiprecision::uint128_t sqrt_un2(sqrt(un2));
+      constexpr boost::multiprecision::uint128_t sqrt_un3(sqrt(un3));
+
+      static_assert(static_cast<std::uint64_t>(sqrt_un0) == UINT64_C(783014315477220686), "Error in constexpr integer square root");
+      static_assert(static_cast<std::uint64_t>(sqrt_un1) == UINT64_C(589847733628674383), "Error in constexpr integer square root");
+      static_assert(static_cast<std::uint64_t>(sqrt_un2) == UINT64_C(444334595153395054), "Error in constexpr integer square root");
+      static_assert(static_cast<std::uint64_t>(sqrt_un3) == UINT64_C(334718981177642730), "Error in constexpr integer square root");
+
+      static_assert(sqrt(boost::multiprecision::uint128_t(0)) == 0, "Error in constexpr integer square root");
+      static_assert(sqrt(boost::multiprecision::uint128_t(1)) == 1, "Error in constexpr integer square root");
+      static_assert(sqrt(boost::multiprecision::uint128_t(2)) == 1, "Error in constexpr integer square root");
+
+      // N[Sqrt[2^128 - 1], 50]
+      // 1.8446744073709551615999999999999999999972894945688*10^19
+      static_assert(sqrt((std::numeric_limits<boost::multiprecision::uint128_t>::max)()) == boost::multiprecision::uint128_t(UINT64_C(18446744073709551615)), "Error in constexpr integer square root");
+
+      BOOST_TEST(sqrt(boost::multiprecision::uint128_t("613111418242260482895474318171556018")) == boost::multiprecision::uint128_t(UINT64_C(783014315477220686)));
+      BOOST_TEST(sqrt(boost::multiprecision::uint128_t("347920348866883608793309754527486715")) == boost::multiprecision::uint128_t(UINT64_C(589847733628674383)));
+      BOOST_TEST(sqrt(boost::multiprecision::uint128_t("197433232450131483625758030703105499")) == boost::multiprecision::uint128_t(UINT64_C(444334595153395054)));
+      BOOST_TEST(sqrt(boost::multiprecision::uint128_t("112036796360599148200087404494586705")) == boost::multiprecision::uint128_t(UINT64_C(334718981177642730)));
+
+      BOOST_TEST(sqrt(boost::multiprecision::uint128_t("0")) == boost::multiprecision::uint128_t(UINT8_C(0)));
+      BOOST_TEST(sqrt(boost::multiprecision::uint128_t("1")) == boost::multiprecision::uint128_t(UINT8_C(1)));
+      BOOST_TEST(sqrt(boost::multiprecision::uint128_t("2")) == boost::multiprecision::uint128_t(UINT8_C(1)));
+    }
+  }
 } // namespace local
 
 auto main() -> int
 {
+  local::test_sqrt_integral_and_constexpr();
+
   #if defined(TEST_CPP_DOUBLE_FLOAT)
   {
     using float_type = boost::multiprecision::cpp_double_float;
@@ -1535,6 +1691,7 @@ auto main() -> int
 
     static_cast<void>(local::test_edges<float_type>());
     static_cast<void>(local::test_edges_extra<float_type>());
+    local::test_fdim_edge<float_type>();
     local::test_sqrt_edge<float_type>();
     local::test_exp_edge<float_type>();
     local::test_log_edge<float_type>();
@@ -1549,6 +1706,7 @@ auto main() -> int
 
     static_cast<void>(local::test_edges<float_type>());
     static_cast<void>(local::test_edges_extra<float_type>());
+    local::test_fdim_edge<float_type>();
     local::test_sqrt_edge<float_type>();
     local::test_exp_edge<float_type>();
     local::test_log_edge<float_type>();
@@ -1563,6 +1721,7 @@ auto main() -> int
 
     static_cast<void>(local::test_edges<float_type>());
     static_cast<void>(local::test_edges_extra<float_type>());
+    local::test_fdim_edge<float_type>();
     local::test_sqrt_edge<float_type>();
     local::test_exp_edge<float_type>();
     local::test_log_edge<float_type>();
@@ -1580,6 +1739,7 @@ auto main() -> int
     std::cout << "Testing type: " << typeid(float_type).name() << std::endl;
 
     static_cast<void>(local::test_edges<float_type>());
+    local::test_fdim_edge<float_type>();
     local::test_sqrt_edge<float_type>();
     local::test_exp_edge<float_type>();
     local::test_log_edge<float_type>();
@@ -1595,6 +1755,7 @@ auto main() -> int
     std::cout << "Testing type: " << typeid(float_type).name() << std::endl;
 
     static_cast<void>(local::test_edges<float_type>());
+    local::test_fdim_edge<float_type>();
     local::test_sqrt_edge<float_type>();
     local::test_exp_edge<float_type>();
     local::test_log_edge<float_type>();
@@ -1610,6 +1771,7 @@ auto main() -> int
     std::cout << "Testing type: " << typeid(float_type).name() << std::endl;
 
     static_cast<void>(local::test_edges<float_type>());
+    local::test_fdim_edge<float_type>();
     local::test_sqrt_edge<float_type>();
     local::test_exp_edge<float_type>();
     local::test_log_edge<float_type>();
@@ -1626,8 +1788,6 @@ template<typename FloatType> auto my_inf () noexcept -> FloatType& { using float
 template<typename FloatType> auto my_nan () noexcept -> FloatType& { using float_type = FloatType; static float_type val_nan  { std::numeric_limits<float_type>::quiet_NaN() }; return val_nan; }
 template<typename FloatType> auto my_exp1() noexcept -> FloatType& { using float_type = FloatType; static float_type val_exp1 { "2.718281828459045235360287471352662497757247093699959574966967627724076630353547594571382178525166427427466391932003059921817413596629043572900334295260595630738132328627943490763233829880753195251019011573834187930702154089149934884167509244761460668082265" }; return val_exp1; }
 
-#if defined(__clang__)
-#  pragma clang diagnostic pop
-#elif defined(__GNUC__)
+#if (defined(__clang__) || defined(__GNUC__))
 #  pragma GCC diagnostic pop
 #endif
diff --git a/test/test_various_edges_more.cpp b/test/test_various_edges_more.cpp
index a0b6257e..e2a0191d 100644
--- a/test/test_various_edges_more.cpp
+++ b/test/test_various_edges_more.cpp
@@ -5,8 +5,9 @@
 
 // Some parts of this test file have been taken from the Boost.Decimal project.
 
-#include <boost/multiprecision/cpp_dec_float.hpp>
 #include <boost/multiprecision/cpp_bin_float.hpp>
+#include <boost/multiprecision/cpp_dec_float.hpp>
+#include <boost/multiprecision/cpp_double_fp.hpp>
 
 #include <test_traits.hpp> // Note: include this AFTER the test-backends are defined
 
@@ -17,10 +18,7 @@
 #include <string>
 #include <vector>
 
-#if defined(__clang__)
-#  pragma clang diagnostic push
-#  pragma clang diagnostic ignored "-Wfloat-equal"
-#elif defined(__GNUC__)
+#if (defined(__clang__) || defined(__GNUC__))
 #  pragma GCC diagnostic push
 #  pragma GCC diagnostic ignored "-Wfloat-equal"
 #endif
@@ -74,13 +72,79 @@ namespace local
     {
       delta = fabs(a - b); // LCOV_EXCL_LINE
 
-      result_is_ok = (delta < tol); // LCOV_EXCL_LINE
+      result_is_ok = (delta <= tol); // LCOV_EXCL_LINE
     }
     else
     {
       delta = fabs(1 - (a / b));
 
-      result_is_ok = (delta < tol);
+      result_is_ok = (delta <= tol);
+    }
+
+    return result_is_ok;
+  }
+
+  template<class AnyFloatType, class OtherFloatType>
+  bool test_convert_and_back(const float epsilon_factor)
+  {
+    using any_float_type = AnyFloatType;
+    using other_float_type = OtherFloatType;
+
+    std::mt19937_64 gen { time_point<typename std::mt19937_64::result_type>() };
+
+    auto dis =
+      std::uniform_real_distribution<float>
+      {
+        static_cast<float>(-1.5F),
+        static_cast<float>(1.5F)
+      };
+
+    bool result_is_ok { true };
+
+    for(int n_loop = static_cast<int>(INT8_C(0)); n_loop < static_cast<int>(INT8_C(64)); ++n_loop)
+    {
+      static_cast<void>(n_loop);
+
+      using string_data_array_type = std::array<std::string, std::size_t { UINT8_C(5) }>;
+
+      // Table[N[Exp[Pi   EulerGamma*m], 100], {m, 1, 41, 10}]
+      // 6.131114182422604828954743181715560166328314533478636289880930665602805209787080979043183175178343525
+      // 4.601860472890328982970020344164597095017995523702194983334509372565826509463730150117899613177592949*10^8
+      // 3.454041008184686240269652058630805927805404685705614002929256804267342980143530922708514273650048389*10^16
+      // 2.592516517287682590319488890844344913780990651072642232646840608443538857001634005575916140667342752*10^24
+      // 1.945877850460679347837790546862626024227287564971475316486097406542338497232351222296834286913627889*10^32
+
+      const string_data_array_type
+        float_number_strings
+        {
+          std::string("6.131114182422604828954743181715560166328314533478636289880930665602805209787080979043183175178343525"),
+          std::string("4.601860472890328982970020344164597095017995523702194983334509372565826509463730150117899613177592949E8"),
+          std::string("3.454041008184686240269652058630805927805404685705614002929256804267342980143530922708514273650048389E16"),
+          std::string("2.592516517287682590319488890844344913780990651072642232646840608443538857001634005575916140667342752E24"),
+          std::string("1.945877850460679347837790546862626024227287564971475316486097406542338497232351222296834286913627889E32")
+        };
+
+      for(std::size_t index { std::size_t { UINT8_C(0) }}; index < std::tuple_size<string_data_array_type>::value; ++index)
+      {
+        const any_float_type start(any_float_type(float_number_strings[index]) * (dis(gen)* dis(gen)));
+        const other_float_type other(start);
+        const any_float_type backto(other);
+
+        const bool
+          result_of_trip_is_ok
+          {
+            local::is_close_fraction
+            (
+              start,
+              backto,
+              any_float_type(std::numeric_limits<any_float_type>::epsilon() * epsilon_factor)
+            )
+          };
+
+        BOOST_TEST(result_of_trip_is_ok);
+
+        result_is_ok = (result_of_trip_is_ok && result_is_ok);
+      }
     }
 
     return result_is_ok;
@@ -194,12 +258,15 @@ namespace local
         const std::string str_ctrl { str128_maker(flt_n128) };
         const std::string str_n128 { str128_maker(static_cast<float_type>(val_n128)) };
 
+        // This test assumes that boost::int128_type is as wide or wider than signed long long.
+        // Consider using a a kind of if-constexpr to verify this "up front".
+
         result_val_n128_is_ok =
           (
                (str_n128 == str_ctrl)
             && (
-                 (!is_neg) ? (val_n128 > static_cast<boost::int128_type>((std::numeric_limits<signed long long>::max)()))
-                           : (val_n128 < static_cast<boost::int128_type>((std::numeric_limits<signed long long>::max)()))
+                 (!is_neg) ? (val_n128 >= static_cast<boost::int128_type>((std::numeric_limits<signed long long>::max)()))
+                           : (val_n128 <= static_cast<boost::int128_type>((std::numeric_limits<signed long long>::max)()))
                )
           );
 
@@ -291,6 +358,106 @@ namespace local
       result_is_ok = (result_funky_strings_is_ok && result_is_ok);
     }
 
+    {
+      for(auto index = static_cast<unsigned>(UINT8_C(0)); index < static_cast<unsigned>(UINT8_C(16)); ++index)
+      {
+        static_cast<void>(index);
+
+        float_type flt_nrm { ::my_one<float_type>() * dis(gen) };
+
+        const bool is_neg { ((index & 1U) != 0U) };
+
+        if(is_neg)
+        {
+          flt_nrm = -flt_nrm;
+        }
+
+        const float_type flt_zer { ::my_zero<float_type>() * dis(gen) };
+
+        flt_nrm /= flt_zer;
+
+        bool result_div_zero_is_ok { (boost::multiprecision::isinf)(flt_nrm) };
+
+        if(is_neg)
+        {
+          result_div_zero_is_ok = ((boost::multiprecision::signbit)(flt_nrm) && result_div_zero_is_ok);
+        }
+
+        BOOST_TEST(result_div_zero_is_ok);
+
+        result_is_ok = (result_div_zero_is_ok && result_is_ok);
+      }
+    }
+
+    {
+      for(auto index = static_cast<unsigned>(UINT8_C(0)); index < static_cast<unsigned>(UINT8_C(16)); ++index)
+      {
+        static_cast<void>(index);
+
+        float_type flt_nrm { ::my_one<float_type>() + 0.25F };
+        float_type flt_one { ::my_one<float_type>() + ::my_zero<float_type>() * dis(gen) };
+
+        const bool is_neg { ((index & 1U) != 0U) };
+
+        if(is_neg)
+        {
+          flt_one = -flt_one;
+        }
+
+        flt_nrm /= flt_one;
+
+        bool result_div_one_is_ok { (is_neg ? (flt_nrm == -1.25F) : (flt_nrm == 1.25F)) };
+
+        BOOST_TEST(result_div_one_is_ok);
+
+        result_is_ok = (result_div_one_is_ok && result_is_ok);
+      }
+    }
+
+    constexpr unsigned max_index { 1024U };
+
+    {
+      using std::ldexp;
+
+      float_type flt_near_min { ldexp((std::numeric_limits<float_type>::min)(), +64) };
+
+      unsigned index { };
+
+      while((index < max_index) && (!((boost::multiprecision::fpclassify)(flt_near_min) == FP_ZERO)))
+      {
+        flt_near_min /= static_cast<unsigned long long>(0xDEADBEEF);
+
+        ++index;
+      }
+
+      const bool result_unf_is_ok { ((index > 1U) && (index < max_index)) };
+
+      BOOST_TEST(result_unf_is_ok);
+
+      result_is_ok = (result_unf_is_ok && result_is_ok);
+    }
+
+    {
+      using std::ldexp;
+
+      float_type flt_near_min { ldexp((std::numeric_limits<float_type>::min)(), +64) };
+
+      unsigned index { };
+
+      while((index < max_index) && (!((boost::multiprecision::fpclassify)(flt_near_min) == FP_ZERO)))
+      {
+        flt_near_min /= static_cast<unsigned long long>(100000000ULL - 3ULL);
+
+        ++index;
+      }
+
+      const bool result_unf_is_ok { ((index > 1U) && (index < max_index)) };
+
+      BOOST_TEST(result_unf_is_ok);
+
+      result_is_ok = (result_unf_is_ok && result_is_ok);
+    }
+
     {
       using std::ldexp;
 
@@ -300,14 +467,14 @@ namespace local
 
       unsigned index { };
 
-      while((index < 1024U) && (!(boost::multiprecision::isinf)(flt_near_max)))
+      while((index < max_index) && (!(boost::multiprecision::isinf)(flt_near_max)))
       {
         flt_near_max += (flt_less_near_max * dis(gen));
 
         ++index;
       }
 
-      const bool result_ovf_is_ok { ((index > 1U) && (index < 1024U)) };
+      const bool result_ovf_is_ok { ((index > 1U) && (index < max_index)) };
 
       BOOST_TEST(result_ovf_is_ok);
 
@@ -323,14 +490,14 @@ namespace local
 
       unsigned index { };
 
-      while((index < 1024U) && (!(boost::multiprecision::isinf)(flt_near_lowest)))
+      while((index < max_index) && (!(boost::multiprecision::isinf)(flt_near_lowest)))
       {
         flt_near_lowest -= (flt_less_near_max * dis(gen));
 
         ++index;
       }
 
-      const bool result_ovf_is_ok { ((index > 1U) && (index < 1024U)) && signbit(flt_near_lowest) };
+      const bool result_ovf_is_ok { ((index > 1U) && (index < max_index)) && signbit(flt_near_lowest) };
 
       BOOST_TEST(result_ovf_is_ok);
 
@@ -346,14 +513,14 @@ namespace local
 
       unsigned index { };
 
-      while((index < 1024U) && (!(boost::multiprecision::isinf)(flt_near_max)))
+      while((index < max_index) && (!(boost::multiprecision::isinf)(flt_near_max)))
       {
         flt_near_max *= static_cast<unsigned long long>(INT32_C(2));
 
         ++index;
       }
 
-      const bool result_ovf_is_ok { ((index > 1U) && (index < 1024U)) };
+      const bool result_ovf_is_ok { ((index > 1U) && (index < max_index)) };
 
       BOOST_TEST(result_ovf_is_ok);
 
@@ -387,21 +554,34 @@ namespace local
       {
         static_cast<void>(index);
 
-        float_type flt_zero_numerator { ::my_zero<float_type>() * dis(gen) };
-        float_type flt_finite_numerator { ::my_one<float_type>() * dis(gen) };
+        float_type flt_zer_numerator { ::my_zero<float_type>() * dis(gen) };
+        float_type flt_nrm_numerator { ::my_one<float_type>() * dis(gen) };
         float_type flt_nan_numerator { ::my_nan<float_type>() * dis(gen) };
+        float_type flt_inf_numerator { ::my_inf<float_type>() * dis(gen) };
 
-        const float_type flt_infinite_result_neg { -flt_finite_numerator /= static_cast<signed long long>(INT8_C(0)) };
+        const float_type flt_inf_result_neg_01 { -flt_nrm_numerator /= static_cast<signed long long>(INT8_C(0)) };
+        const float_type flt_inf_result_neg_02 { -flt_inf_numerator *= static_cast<signed long long>(index + 2U) };
 
-        const bool result_edge_00 { (boost::multiprecision::isnan)(flt_zero_numerator /= static_cast<signed long long>(INT8_C(0))) };
-        const bool result_edge_01 { (boost::multiprecision::isinf)(flt_finite_numerator /= static_cast<signed long long>(INT8_C(0))) };
-        const bool result_edge_02 { (boost::multiprecision::isinf)(flt_infinite_result_neg) && (boost::multiprecision::signbit)(flt_infinite_result_neg) };
+        const bool result_edge_00 { (boost::multiprecision::isnan)(flt_zer_numerator /= static_cast<signed long long>(INT8_C(0))) };
+        const bool result_edge_01 { (boost::multiprecision::isinf)(flt_nrm_numerator /= static_cast<signed long long>(INT8_C(0))) };
+        const bool result_edge_02 { (boost::multiprecision::isinf)(flt_inf_result_neg_01) && (boost::multiprecision::signbit)(flt_inf_result_neg_01) };
         const bool result_edge_03 { (boost::multiprecision::isnan)(flt_nan_numerator /= static_cast<signed long long>(index + 2U)) };
+        const bool result_edge_04 { (boost::multiprecision::isnan)(flt_nan_numerator *= static_cast<signed long long>(index + 2U)) };
+        const bool result_edge_05 { (boost::multiprecision::isnan)(flt_inf_numerator *= static_cast<signed long long>(0)) };
+        const bool result_edge_06 { (boost::multiprecision::isinf)(flt_inf_result_neg_02) && (boost::multiprecision::signbit)(flt_inf_result_neg_02) };
 
         const bool
           result_divs_are_ok
           {
-            (result_edge_00 && result_edge_01 && result_edge_02 && result_edge_03)
+            (
+                 result_edge_00
+              && result_edge_01
+              && result_edge_02
+              && result_edge_03
+              && result_edge_04
+              && result_edge_05
+              && result_edge_06
+            )
           };
 
         BOOST_TEST(result_divs_are_ok);
@@ -413,6 +593,23 @@ namespace local
     return result_is_ok;
   }
 
+  template<class OtherFloatType>
+  auto test_convert_and_back_caller(const float epsilon_factor = 0.0F) -> bool
+  {
+    using other_float_type = OtherFloatType;
+
+    std::cout << "Testing type (in test_convert_and_back): " << typeid(other_float_type).name() << std::endl;
+
+    const bool
+      result_of_trip_is_ok
+      {
+        local::test_convert_and_back<boost::multiprecision::cpp_double_double, other_float_type>(epsilon_factor)
+      };
+
+    BOOST_TEST(result_of_trip_is_ok);
+
+    return result_of_trip_is_ok;
+  }
 } // namespace local
 
 auto main() -> int
@@ -437,6 +634,17 @@ auto main() -> int
     static_cast<void>(local::test_edges<float_type>());
   }
 
+  {
+    using big_bin_float_backend_type = boost::multiprecision::cpp_bin_float<512>;
+    using big_dec_float_backend_type = boost::multiprecision::cpp_dec_float<512>;
+
+    using big_bin_float_type = boost::multiprecision::number<big_bin_float_backend_type, boost::multiprecision::et_off>;
+    using big_dec_float_type = boost::multiprecision::number<big_dec_float_backend_type, boost::multiprecision::et_off>;
+
+    static_cast<void>(local::test_convert_and_back_caller<big_bin_float_type>(0.25F));
+    static_cast<void>(local::test_convert_and_back_caller<big_dec_float_type>(0.50F));
+  }
+
   return boost::report_errors();
 }
 
@@ -445,8 +653,6 @@ template<typename FloatType> auto my_one () noexcept -> FloatType& { using float
 template<typename FloatType> auto my_inf () noexcept -> FloatType& { using float_type = FloatType; static float_type val_inf  { std::numeric_limits<float_type>::infinity() }; return val_inf; }
 template<typename FloatType> auto my_nan () noexcept -> FloatType& { using float_type = FloatType; static float_type val_nan  { std::numeric_limits<float_type>::quiet_NaN() }; return val_nan; }
 
-#if defined(__clang__)
-#  pragma clang diagnostic pop
-#elif defined(__GNUC__)
+#if (defined(__clang__) || defined(__GNUC__))
 #  pragma GCC diagnostic pop
 #endif