Initial commit log10 but more tests needed

Author: ckormanyos

include/boost/decimal/detail/cmath/log10.hpp

@@ -7,9 +7,10 @@
 #define BOOST_DECIMAL_DETAIL_CMATH_LOG10_HPP
 
 #include <boost/decimal/fwd.hpp> // NOLINT(llvm-include-order)
-#include <boost/decimal/detail/type_traits.hpp>
+#include <boost/decimal/detail/cmath/impl/log_impl.hpp>
 #include <boost/decimal/detail/concepts.hpp>
 #include <boost/decimal/detail/config.hpp>
+#include <boost/decimal/detail/type_traits.hpp>
 #include <boost/decimal/numbers.hpp>
 
 #ifndef BOOST_DECIMAL_BUILD_MODULE
@@ -26,21 +27,102 @@ template <typename T>
 constexpr auto log10_impl(T x) noexcept
     BOOST_DECIMAL_REQUIRES(detail::is_decimal_floating_point_v, T)
 {
-    // TODO(ckormanyos) Actually this is eintirely preliminary. The implementation
-    // of log10 will/should-be entirely re-worked with specific base-10-relevant details.
-
-    // TODO(ckormanyos) Handle non-normal varguments.
-
-    // TODO(ckormanyos) Put in a basic check for pure powers of 10, resulting
-    // in an exact result.
-
-    // Starting point: See implementation of log() and adapt to the following series.
-    // Series[Log[10] Log[10, (1 + (z/2))/(1 - (z/2))], {z, 0, 17}]
-
-    return log(x) / numbers::ln10_v<T>;
+    constexpr T zero { 0, 0 };
+    constexpr T one  { 1, 0 };
+
+    T result { };
+
+    if (isnan(x))
+    {
+        result = x;
+    }
+    else if (isinf(x))
+    {
+        result = (!signbit(x)) ? x: std::numeric_limits<T>::quiet_NaN();
+    }
+    else if (x < one)
+    {
+        // Handle reflection, the [+/-] zero-pole, and non-pole, negative x.
+        if (x > zero)
+        {
+            result = -log10(one / x);
+        }
+        else if ((x == zero) || (-x == zero))
+        {
+            // Actually, this should be equivalent to -HUGE_VAL.
+
+            result = -std::numeric_limits<T>::infinity();
+        }
+        else
+        {
+            result = std::numeric_limits<T>::quiet_NaN();
+        }
+    }
+    else if(x > one)
+    {
+        // The algorithm for base-10 logarithm is based on Chapter 5, pages 35-36
+        // of Cody and Waite, Software Manual for the Elementary Functions,
+        // Prentice Hall, 1980.
+
+        int exp10val { };
+
+        T g { frexp10(x, &exp10val) };
+
+        while(g > one)
+        {
+            // TODO(ckormanyos) Do we really have to multiply individual powers
+            // of 10 here? Or is there a reliable way to simply use pow10()?
+            // For instance count intagral orders of exp10val and divide by pow10().
+
+            g /= 10;
+
+            ++exp10val;
+        }
+
+        if(g == one)
+        {
+            result = T { exp10val };
+        }
+        else
+        {
+            constexpr T inv_sqrt10 { UINT64_C(3162277660168379332), -19 };
+
+            const bool reduce_sqrt10 { g < inv_sqrt10 };
+
+            if (reduce_sqrt10)
+            {
+                constexpr T sqrt10 { UINT64_C(3162277660168379332), -18 };
+
+                g *= sqrt10;
+            }
+
+            const T s   { (g - one) / (g + one) };
+            const T z   { s + s };
+            const T zsq { z * z };
+
+            result = z * fma(detail::log_series_expansion(zsq), zsq, one);
+
+            result /= numbers::ln10_v<T>;
+
+            if(reduce_sqrt10)
+            {
+                constexpr T half { 5, -1 };
+
+                result -= half;
+            }
+
+            result += static_cast<T>(exp10val);
+        }
+    }
+    else
+    {
+        result = zero;
+    }
+
+    return result;
 }
 
-} // namespace detail
+} //namespace detail
 
 BOOST_DECIMAL_EXPORT template <typename T>
 constexpr auto log10(T x) noexcept

test/Jamfile

@@ -106,6 +106,7 @@ run test_legendre.cpp ;
 run test_literals.cpp ;
 run test_lgamma.cpp ;
 run test_log.cpp ;
+run test_log10.cpp ;
 run test_log1p.cpp ;
 run test_pow.cpp ;
 run test_promotion.cpp ;

test/test_log10.cpp

@@ -0,0 +1,200 @@
+// Copyright 2024 Matt Borland
+// Copyright 2024 Christopher Kormanyos
+// Distributed under the Boost Software License, Version 1.0.
+// https://www.boost.org/LICENSE_1_0.txt
+
+#include <chrono>
+#include <iomanip>
+#include <iostream>
+#include <limits>
+#include <random>
+#include <string>
+
+#include <boost/decimal.hpp>
+
+#if defined(__clang__)
+#  pragma clang diagnostic push
+#  pragma clang diagnostic ignored "-Wfloat-equal"
+#elif defined(__GNUC__)
+#  pragma GCC diagnostic push
+#  pragma GCC diagnostic ignored "-Wfloat-equal"
+#endif
+
+#include <boost/core/lightweight_test.hpp>
+
+namespace local
+{
+  template<typename IntegralTimePointType,
+           typename ClockType = std::chrono::high_resolution_clock>
+  auto time_point() noexcept -> IntegralTimePointType
+  {
+    using local_integral_time_point_type = IntegralTimePointType;
+    using local_clock_type               = ClockType;
+
+    const auto current_now =
+      static_cast<std::uintmax_t>
+      (
+        std::chrono::duration_cast<std::chrono::nanoseconds>
+        (
+          local_clock_type::now().time_since_epoch()
+        ).count()
+      );
+
+    return static_cast<local_integral_time_point_type>(current_now);
+  }
+
+  template<typename NumericType>
+  auto is_close_fraction(const NumericType& a,
+                         const NumericType& b,
+                         const NumericType& tol) noexcept -> bool
+  {
+    using std::fabs;
+
+    auto result_is_ok = bool { };
+
+    NumericType delta { };
+
+    if(b == static_cast<NumericType>(0))
+    {
+      delta = fabs(a - b); // LCOV_EXCL_LINE
+
+      result_is_ok = (delta < tol); // LCOV_EXCL_LINE
+    }
+    else
+    {
+      delta = fabs(1 - (a / b));
+
+      result_is_ok = (delta < tol);
+    }
+
+    // LCOV_EXCL_START
+    if (!result_is_ok)
+    {
+      std::cerr << std::setprecision(std::numeric_limits<NumericType>::digits10) << "a: " << a
+                << "\nb: " << b
+                << "\ndelta: " << delta
+                << "\ntol: " << tol << std::endl;
+    }
+    // LCOV_EXCL_STOP
+
+    return result_is_ok;
+  }
+
+  template<typename DecimalType, typename FloatType>
+  auto test_log10(const int tol_factor) -> bool
+  {
+    using decimal_type = DecimalType;
+    using float_type   = FloatType;
+
+    std::random_device rd;
+    std::mt19937_64 gen(rd());
+
+    gen.seed(time_point<typename std::mt19937_64::result_type>());
+
+    auto dis_r =
+      std::uniform_real_distribution<float_type>
+      {
+        static_cast<float_type>(1.2L),
+        static_cast<float_type>(8.9L)
+      };
+
+    bool result_is_ok { true };
+
+    auto trials = static_cast<std::uint32_t>(UINT8_C(0));
+
+    #if !defined(BOOST_DECIMAL_REDUCE_TEST_DEPTH)
+    constexpr auto count = (sizeof(decimal_type) == static_cast<std::size_t>(UINT8_C(4))) ? static_cast<std::uint32_t>(UINT32_C(0x200)) : static_cast<std::uint32_t>(UINT32_C(0x40));
+    #else
+    constexpr auto count = (sizeof(decimal_type) == static_cast<std::size_t>(UINT8_C(4))) ? static_cast<std::uint32_t>(UINT32_C(0x40)) : static_cast<std::uint32_t>(UINT32_C(0x4));
+    #endif
+
+    for( ; trials < count; ++trials)
+    {
+      auto x_flt = dis_r(gen);
+
+      auto dis_n =
+        std::uniform_int_distribution<int>
+        {
+          static_cast<int>(INT8_C(-17)),
+          static_cast<int>(INT8_C(17))
+        };
+
+      std::string str_e { "1.0E" + std::to_string(dis_n(gen)) };
+
+      char* p_end;
+
+      x_flt *= static_cast<float>(strtold(str_e.c_str(), &p_end));
+
+      static_cast<void>(p_end);
+
+      const auto x_dec = static_cast<decimal_type>(x_flt);
+
+      using std::log10;
+
+      const auto val_flt = log10(x_flt);
+      const auto val_dec = log10(x_dec);
+
+      const auto result_log_is_ok = is_close_fraction(val_flt, static_cast<float_type>(val_dec), std::numeric_limits<float_type>::epsilon() * static_cast<float_type>(tol_factor));
+
+      result_is_ok = (result_log_is_ok && result_is_ok);
+
+      if(!result_log_is_ok)
+      {
+        // LCOV_EXCL_START
+        std::cerr << "x_flt  : " << std::scientific << std::setprecision(std::numeric_limits<float_type>::digits10) << x_flt   << std::endl;
+        std::cerr << "val_flt: " << std::scientific << std::setprecision(std::numeric_limits<float_type>::digits10) << val_flt << std::endl;
+        std::cerr << "val_dec: " << std::scientific << std::setprecision(std::numeric_limits<float_type>::digits10) << val_dec << std::endl;
+
+        break;
+        // LCOV_EXCL_STOP
+      }
+    }
+
+    BOOST_TEST(result_is_ok);
+
+    return result_is_ok;
+  }
+
+  template<typename DecimalType, typename FloatType>
+  auto test_log10_pow10(const int tol_factor) -> bool
+  {
+    using decimal_type = DecimalType;
+    using float_type   = FloatType;
+
+    bool result_is_ok { true };
+
+    for(int i = static_cast<int>(INT8_C(-23)); i <= static_cast<int>(INT8_C(23)); ++i)
+    {
+      const decimal_type x_arg { 1, i };
+
+      const decimal_type log_val_p10 = log10(x_arg);
+
+      const auto result_log10_pow10_is_ok = (static_cast<float_type>(log_val_p10) == static_cast<float_type>(i));
+
+      result_is_ok = (result_log10_pow10_is_ok && result_is_ok);
+    }
+
+    BOOST_TEST(result_is_ok);
+
+    return result_is_ok;
+  }
+
+} // namespace local
+
+auto main() -> int
+{
+  auto result_is_ok = true;
+
+  {
+    using decimal_type = boost::decimal::decimal32;
+    using float_type   = float;
+
+    const auto test_log10_is_ok = local::test_log10<decimal_type, float_type>(128);
+
+    result_is_ok = (test_log10_is_ok && result_is_ok);
+  }
+
+  result_is_ok = ((boost::report_errors() == 0) && result_is_ok);
+
+  return (result_is_ok ? 0 : -1);
+}