Try finishing constexpr-ness

include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_log.hpp

@@ -8,28 +8,162 @@
 #ifndef BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_LOG_2024_12_30_HPP
 #define BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_LOG_2024_12_30_HPP
 
+#include <boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_frexp.hpp>
+#include <boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_limits.hpp>
+
 #include <cmath>
 #include <type_traits>
 
+#if (defined(__GNUC__) && defined(BOOST_MP_CPP_DOUBLE_FP_HAS_FLOAT128))
+//
+// This is the only way we can avoid
+// warning: non-standard suffix on floating constant [-Wpedantic]
+// when building with -Wall -pedantic.  Neither __extension__
+// nor #pragma diagnostic ignored work :(
+//
+#pragma GCC system_header
+#endif
+
 namespace boost { namespace multiprecision { namespace backends { namespace cpp_df_qf_detail { namespace ccmath {
 
 namespace detail {
 
-template <class T>
-constexpr auto log_impl(T x) -> T
+// LCOV_EXCL_START
+template <class Real>
+constexpr auto exp_impl(Real x) noexcept -> Real
 {
-   // Default to the regular log function.
-   using std::log;
+   constexpr int
+     my_digits_10
+     {
+        ::boost::multiprecision::backends::cpp_df_qf_detail::ccmath::numeric_limits<Real>::digits10
+     };
 
-   return log(x);
+   constexpr int
+     loop_count
+     {
+          my_digits_10 <  8 ?  7
+        : my_digits_10 < 16 ? 15
+        : my_digits_10 < 20 ? 19
+        :                     33
+     };
+
+   // Scale the argument with a single factor of 2.
+   // Then square the result upon return.
+   x /= 2;
+
+   // Perform a simple Taylor series of the exponent function.
+   Real term { x };
+   Real sum { Real { 1 } + term };
+
+   for (int loop_index { INT8_C(2) }; loop_index < loop_count; ++loop_index)
+   {
+      term *= x;
+      term /= static_cast<Real>(loop_index);
+      sum += term;
+   }
+
+   // Scale the result.
+   return sum * sum;
 }
 
+template<typename Real>
+constexpr auto log_impl_pade(Real x) noexcept -> typename std::enable_if<(::boost::multiprecision::backends::cpp_df_qf_detail::ccmath::numeric_limits<Real>::digits < 54), Real>::type
+{
+   // PadeApproximant[Log[x], {x, 1, {8, 8}}]
+   // FullSimplify[%]
+
+   const Real
+      top
+      {
+         ((static_cast<Real>(-1.0L) + x) * (static_cast<Real>(1.0L) + x) * (static_cast<Real>(761.0L) + x * (static_cast<Real>(28544.0L) + x * (static_cast<Real>(209305.0L) + x * (static_cast<Real>(423680.0L) + x * (static_cast<Real>(209305.0L) + x * (static_cast<Real>(28544.0L) + static_cast<Real>(761.0L) * x)))))))
+      };
+
+   const Real
+      bot
+      {
+         (static_cast<Real>(140.0L) * (static_cast<Real>(1.0L) + x * (static_cast<Real>(64.0L) + x * (static_cast<Real>(784.0L) + x * (static_cast<Real>(3136.0L) + x * (static_cast<Real>(4900.0L) + x * (static_cast<Real>(3136.0L) + x * (static_cast<Real>(784.0L) + x * (static_cast<Real>(64.0L) + x)))))))))
+      };
+
+   return top / bot;
+}
+
+template<typename Real>
+constexpr auto log_impl_pade(Real x) noexcept -> typename std::enable_if<(!(::boost::multiprecision::backends::cpp_df_qf_detail::ccmath::numeric_limits<Real>::digits < 54)), Real>::type
+{
+   // PadeApproximant[Log[x], {x, 1, {16, 16}}]
+   // FullSimplify[%]
+
+   const Real
+      top
+      {
+         (static_cast<Real>(17.0L) * (static_cast<Real>(-1.0L) + x) * (static_cast<Real>(1.0L) + x) * (static_cast<Real>(143327.0L) + x * (static_cast<Real>(25160192.0L) + x * (static_cast<Real>(1069458527.0L) + x * (static_cast<Real>(17931092992.0L) + x * (static_cast<Real>(144291009727.0L) + x * (static_cast<Real>(613705186816.0L) + x * (static_cast<Real>(1446475477311.0L) + x * (static_cast<Real>(1923749922816.0L) + x * (static_cast<Real>(1446475477311.0L) + x * (static_cast<Real>(613705186816.0L) + x * (static_cast<Real>(144291009727.0L) + x * (static_cast<Real>(17931092992.0L) + x * (static_cast<Real>(1069458527.0L) + x * (static_cast<Real>(25160192.0L) + static_cast<Real>(143327.0L) * x)))))))))))))))
+      };
+
+   const Real
+      bot
+      {
+         (static_cast<Real>(360360.0L) * (static_cast<Real>(1.0L) + x * (static_cast<Real>(256.0L) + x * (static_cast<Real>(14400.0L) + x * (static_cast<Real>(313600.0L) + x * (static_cast<Real>(3312400.0L) + x * (static_cast<Real>(19079424.0L) + x * (static_cast<Real>(64128064.0L) + x * (static_cast<Real>(130873600.0L) + x * (static_cast<Real>(165636900.0L) + x * (static_cast<Real>(130873600.0L) + x * (static_cast<Real>(64128064.0L) + x * (static_cast<Real>(19079424.0L) + x * (static_cast<Real>(3312400.0L) + x * (static_cast<Real>(313600.0L) + x * (static_cast<Real>(14400.0L) + x * (static_cast<Real>(256.0L) + x)))))))))))))))))
+      };
+
+   return top / bot;
+}
+
+// N[Log[2], 101]
+// 0.69314718055994530941723212145817656807550013436025525412068000949339362196969471560586332699641868754
+
+template <typename FloatingPointType> constexpr auto constant_ln_two() noexcept -> typename ::std::enable_if<(::boost::multiprecision::backends::cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits ==  24), FloatingPointType>::type { return static_cast<FloatingPointType>(0.69314718055994530941723212145817656807550013436025525412068L); }
+template <typename FloatingPointType> constexpr auto constant_ln_two() noexcept -> typename ::std::enable_if<(::boost::multiprecision::backends::cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits ==  53), FloatingPointType>::type { return static_cast<FloatingPointType>(0.69314718055994530941723212145817656807550013436025525412068L); }
+template <typename FloatingPointType> constexpr auto constant_ln_two() noexcept -> typename ::std::enable_if<(::boost::multiprecision::backends::cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits ==  64), FloatingPointType>::type { return static_cast<FloatingPointType>(0.69314718055994530941723212145817656807550013436025525412068L); }
+#if defined(BOOST_MP_CPP_DOUBLE_FP_HAS_FLOAT128)
+template <typename FloatingPointType> constexpr auto constant_ln_two() noexcept -> typename ::std::enable_if<(::boost::multiprecision::backends::cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits == 113), FloatingPointType>::type { return static_cast<FloatingPointType>(0.69314718055994530941723212145817656807550013436025525412068Q); }
+#else
+template <typename FloatingPointType> constexpr auto constant_ln_two() noexcept -> typename ::std::enable_if<(::boost::multiprecision::backends::cpp_df_qf_detail::ccmath::numeric_limits<FloatingPointType>::digits == 113), FloatingPointType>::type { return static_cast<FloatingPointType>(0.69314718055994530941723212145817656807550013436025525412068L); }
+#endif
+
+template<typename Real>
+constexpr auto log_impl(Real x) noexcept -> Real
+{
+   int n2 { };
+
+   Real x2 { ::boost::multiprecision::backends::cpp_df_qf_detail::ccmath::detail::frexp_impl(x, &n2) };
+
+   if (x2 > static_cast<Real>(0.875L))
+   {
+     x2 /= 2;
+
+     ++n2;
+   }
+
+   Real s { log_impl_pade(x2) };
+
+   const Real E { exp_impl(s) };
+
+   // Do one single step of Newton-Raphson iteration.
+   s = s + ((x2 - E) / E);
+
+   Real xn2 { static_cast<Real>(n2) * constant_ln_two<Real>() };
+
+   s += xn2;
+
+   return s;
+}
+// LCOV_EXCL_STOP
+
 } // namespace detail
 
 template <typename Real>
 constexpr auto log(Real x) -> Real
 {
-   return cpp_df_qf_detail::ccmath::detail::log_impl<Real>(x);
+   if (BOOST_MP_IS_CONST_EVALUATED(x))
+   {
+      return detail::log_impl<Real>(x); // LCOV_EXCL_LINE
+   }
+   else
+   {
+      using std::log;
+
+      return log(x);
+   }
 }
 
 } } } } } // namespace boost::multiprecision::backends::cpp_df_qf_detail::ccmath

include/boost/multiprecision/cpp_df_qf/cpp_df_qf_detail_ccmath_sqrt.hpp

@@ -52,12 +52,6 @@ constexpr auto sqrt(Real x) noexcept -> Real
    }
 }
 
-template <typename Real>
-constexpr auto sqrt_constexpr(Real x) -> Real
-{
-   return detail::sqrt_impl<Real>(x);
-}
-
 } } } } } // namespace boost::multiprecision::backends::cpp_df_qf_detail::ccmath
 
 #endif // BOOST_MP_CPP_DF_QF_DETAIL_CCMATH_SQRT_2023_01_07_HPP

include/boost/multiprecision/cpp_double_fp.hpp

@@ -943,26 +943,27 @@ class cpp_double_fp_backend
       // Use the non-normalized sum of two maximum values, where the lower
       // value is "shifted" right in the sense of floating-point ldexp.
 
-      return cpp_double_fp_backend
+      return
+         cpp_double_fp_backend
          (
             arithmetic::two_hilo_sum
             (
                float_type
-                  {
-                       (cpp_df_qf_detail::ccmath::numeric_limits<float_type>::max)()
-                     * (
-                            static_cast<float_type>(1.0F)
-                          - static_cast<float_type>(1.5F) * cpp_df_qf_detail::ccmath::sqrt(cpp_df_qf_detail::ccmath::numeric_limits<float_type>::epsilon())
-                       )
-                  },
+               (
+                    (cpp_df_qf_detail::ccmath::numeric_limits<float_type>::max)()
+                  * (
+                         static_cast<float_type>(1.0F)
+                       - static_cast<float_type>(1.5F) * cpp_df_qf_detail::ccmath::sqrt(cpp_df_qf_detail::ccmath::numeric_limits<float_type>::epsilon())
+                    )
+               ),
                float_type
-                  {
-                     cpp_df_qf_detail::ccmath::ldexp
-                     (
-                        (cpp_df_qf_detail::ccmath::numeric_limits<float_type>::max)(),
-                        -cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits
-                     )
-                  }
+               (
+                  cpp_df_qf_detail::ccmath::ldexp
+                  (
+                     (cpp_df_qf_detail::ccmath::numeric_limits<float_type>::max)(),
+                     -cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits
+                  )
+               )
             )
          );
    }
@@ -972,7 +973,8 @@ class cpp_double_fp_backend
       // Use the non-normalized minimum value, where the lower value
       // is "shifted" left in the sense of floating-point ldexp.
 
-      return cpp_double_fp_backend
+      return
+         cpp_double_fp_backend
          (
             float_type
             (
@@ -987,7 +989,8 @@ class cpp_double_fp_backend
 
    static constexpr auto my_value_eps() noexcept -> cpp_double_fp_backend
    {
-      return cpp_double_fp_backend
+      return
+         cpp_double_fp_backend
          (
             float_type(cpp_df_qf_detail::ccmath::ldexp(float_type { 1 }, int { 3 - my_digits }))
          );
@@ -1005,12 +1008,40 @@ class cpp_double_fp_backend
 
    static constexpr auto my_value_logmax() -> cpp_double_fp_backend
    {
-      return float_type(cpp_df_qf_detail::ccmath::log(my_value_max().my_first()));
+      return
+         cpp_double_fp_backend
+         (
+            cpp_df_qf_detail::ccmath::log
+            (
+               float_type
+               (
+                    (cpp_df_qf_detail::ccmath::numeric_limits<float_type>::max)()
+                  * (
+                         static_cast<float_type>(1.0F)
+                       - static_cast<float_type>(1.5F) * cpp_df_qf_detail::ccmath::sqrt(cpp_df_qf_detail::ccmath::numeric_limits<float_type>::epsilon())
+                    )
+               )
+            )
+         );
    }
 
    static constexpr auto my_value_logmin() -> cpp_double_fp_backend
    {
-      return float_type(cpp_df_qf_detail::ccmath::log(my_value_min().my_first()));
+      return
+         cpp_double_fp_backend
+         (
+            cpp_df_qf_detail::ccmath::log
+            (
+               float_type
+               (
+                  cpp_df_qf_detail::ccmath::ldexp
+                  (
+                     (cpp_df_qf_detail::ccmath::numeric_limits<float_type>::min)(),
+                      cpp_df_qf_detail::ccmath::numeric_limits<float_type>::digits
+                  )
+               )
+            )
+         );
    }
 
  private:

test/test_arithmetic_df.cpp

@@ -30,14 +30,17 @@ auto test_constexpr_ness() noexcept -> void
 
    constexpr local_float_type dd_s { sqrt(dd_a) };
    constexpr local_float_type sqrt_hundred { sqrt(local_float_type(100)) };
+   constexpr local_float_type log_hundred { log(local_float_type(100)) };
 
    static_assert(dd_c > 5, "Error in constexpr multiplication");
    static_assert(dd_s > 1, "Error in constexpr square root");
    static_assert(sqrt_hundred == 10, "Error in constexpr square root");
+   static_assert(log_hundred > local_float_type(4.605L), "Error in constexpr logarithm");
 
    BOOST_CHECK(dd_c > 5);
    BOOST_CHECK(dd_s > 1);
    BOOST_CHECK(sqrt_hundred == 10);
+   BOOST_CHECK(log_hundred > local_float_type(4.605L));
 }
 
 int main()