Log10

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,119 @@ 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 { };
+
+    const auto fpc = fpclassify(x);
+
+    if (fpc == FP_ZERO)
+    {
+        result = -std::numeric_limits<T>::infinity();
+    }
+    else if (signbit(x) || (fpc == FP_NAN))
+    {
+        result = std::numeric_limits<T>::quiet_NaN();
+    }
+    else if (fpc == FP_INFINITE)
+    {
+        result = std::numeric_limits<T>::infinity();
+    }
+    else
+    {
+        int exp10val { };
+
+        using representation_type = std::conditional_t<std::is_same<T, decimal32>::value, std::uint32_t,
+                                    std::conditional_t<std::is_same<T, decimal64>::value, std::uint64_t, detail::uint128>>;
+
+        const representation_type gn { frexp10(x, &exp10val) };
+
+        const remove_trailing_zeros_return<representation_type>
+            zeros_removal
+            {
+                remove_trailing_zeros(gn)
+            };
+
+        const bool is_pure
+        {
+            (zeros_removal.trimmed_number == static_cast<representation_type>(UINT8_C(1)))
+        };
+
+        if(is_pure)
+        {
+            // Here, a pure power-of-10 argument gets a pure integral result.
+            const int p10 { exp10val + static_cast<int>(zeros_removal.number_of_removed_zeros) };
+
+            result = T { p10 };
+        }
+        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))
+                {
+                    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.
+
+                // In this implementation, however, we use 2s (as for natural
+                // logarithm in Cody and Waite) even though we are computing
+                // the base-10 logarithm.
+
+                T g { gn, -std::numeric_limits<T>::digits10 };
+
+                exp10val += std::numeric_limits<T>::digits10;
+
+                int reduce_sqrt2 { };
+
+                while (g < numbers::inv_sqrt2_v<T>)
+                {
+                    g += g;
+
+                    ++reduce_sqrt2;
+                }
+
+                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>;
+
+                for(int i = 0; i < reduce_sqrt2; ++i)
+                {
+                    result -= numbers::log10_2_v<T>;
+                }
+
+                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

include/boost/decimal/numbers.hpp

@@ -38,6 +38,13 @@ BOOST_DECIMAL_EXPORT template <>
 BOOST_DECIMAL_CONSTEXPR_VARIABLE_SPECIALIZATION decimal128 log10e_v<decimal128> = decimal128{detail::uint128{UINT64_C(235431510388986),
                                                                                                              UINT64_C(2047877485384264674)}, -34};
 
+BOOST_DECIMAL_EXPORT template <BOOST_DECIMAL_DECIMAL_FLOATING_TYPE Dec, std::enable_if_t<detail::is_decimal_floating_point_v<Dec>, bool> = true>
+BOOST_DECIMAL_CONSTEXPR_VARIABLE Dec log10_2_v = Dec{UINT64_C(3010299956639811952), -19};
+
+BOOST_DECIMAL_EXPORT template <>
+BOOST_DECIMAL_CONSTEXPR_VARIABLE_SPECIALIZATION decimal128 log10_2_v<decimal128> = decimal128{detail::uint128{UINT64_C(163188687641095),
+                                                                                                              UINT64_C(3612628795761985410)}, -34};
+
 BOOST_DECIMAL_EXPORT template <BOOST_DECIMAL_DECIMAL_FLOATING_TYPE Dec, std::enable_if_t<detail::is_decimal_floating_point_v<Dec>, bool> = true>
 BOOST_DECIMAL_CONSTEXPR_VARIABLE Dec pi_v = Dec{UINT64_C(3141592653589793238), -18};
 
@@ -127,8 +134,9 @@ BOOST_DECIMAL_CONSTEXPR_VARIABLE_SPECIALIZATION decimal128 phi_v<decimal128> = d
                                                                                                           UINT64_C(2061523135646567614)}, -33};
 
 BOOST_DECIMAL_EXPORT BOOST_DECIMAL_CONSTEXPR_VARIABLE auto e {e_v<decimal64>};
-BOOST_DECIMAL_EXPORT BOOST_DECIMAL_CONSTEXPR_VARIABLE auto log2e {log2e_v<decimal64>};
+BOOST_DECIMAL_EXPORT BOOST_DECIMAL_CONSTEXPR_VARIABLE auto log10_2 {log10_2_v<decimal64>};
 BOOST_DECIMAL_EXPORT BOOST_DECIMAL_CONSTEXPR_VARIABLE auto log10e {log10e_v<decimal64>};
+BOOST_DECIMAL_EXPORT BOOST_DECIMAL_CONSTEXPR_VARIABLE auto log2e {log2e_v<decimal64>};
 BOOST_DECIMAL_EXPORT BOOST_DECIMAL_CONSTEXPR_VARIABLE auto pi {pi_v<decimal64>};
 BOOST_DECIMAL_EXPORT BOOST_DECIMAL_CONSTEXPR_VARIABLE auto inv_pi {inv_pi_v<decimal64>};
 BOOST_DECIMAL_EXPORT BOOST_DECIMAL_CONSTEXPR_VARIABLE auto inv_sqrtpi {inv_sqrtpi_v<decimal64>};

test/Jamfile

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