77#define BOOST_DECIMAL_DETAIL_CMATH_LOG10_HPP
88
99#include < boost/decimal/fwd.hpp> // NOLINT(llvm-include-order)
10- #include < boost/decimal/detail/type_traits .hpp>
10+ #include < boost/decimal/detail/cmath/impl/log_impl .hpp>
1111#include < boost/decimal/detail/concepts.hpp>
1212#include < boost/decimal/detail/config.hpp>
13+ #include < boost/decimal/detail/type_traits.hpp>
1314#include < boost/decimal/numbers.hpp>
1415
1516#ifndef BOOST_DECIMAL_BUILD_MODULE
@@ -26,21 +27,102 @@ template <typename T>
2627constexpr auto log10_impl (T x) noexcept
2728 BOOST_DECIMAL_REQUIRES(detail::is_decimal_floating_point_v, T)
2829{
29- // TODO(ckormanyos) Actually this is eintirely preliminary. The implementation
30- // of log10 will/should-be entirely re-worked with specific base-10-relevant details.
31-
32- // TODO(ckormanyos) Handle non-normal varguments.
33-
34- // TODO(ckormanyos) Put in a basic check for pure powers of 10, resulting
35- // in an exact result.
36-
37- // Starting point: See implementation of log() and adapt to the following series.
38- // Series[Log[10] Log[10, (1 + (z/2))/(1 - (z/2))], {z, 0, 17}]
39-
40- return log (x) / numbers::ln10_v<T>;
30+ constexpr T zero { 0 , 0 };
31+ constexpr T one { 1 , 0 };
32+
33+ T result { };
34+
35+ const auto fpc = fpclassify (x);
36+
37+ if (fpc == FP_ZERO)
38+ {
39+ result = -std::numeric_limits<T>::infinity ();
40+ }
41+ else if (signbit (x) || (fpc == FP_NAN))
42+ {
43+ result = std::numeric_limits<T>::quiet_NaN ();
44+ }
45+ else if (fpc == FP_INFINITE)
46+ {
47+ result = std::numeric_limits<T>::infinity ();
48+ }
49+ else
50+ {
51+ int exp10val { };
52+
53+ const auto gn { frexp10 (x, &exp10val) };
54+
55+ const auto
56+ zeros_removal
57+ {
58+ remove_trailing_zeros (gn)
59+ };
60+
61+ const bool is_pure { static_cast <unsigned >(zeros_removal.trimmed_number ) == 1U };
62+
63+ if (is_pure)
64+ {
65+ // Here, a pure power-of-10 argument gets a pure integral result.
66+ const int p10 { exp10val + static_cast <int >(zeros_removal.number_of_removed_zeros ) };
67+
68+ result = T { p10 };
69+ }
70+ else
71+ {
72+ if (x < one)
73+ {
74+ // Handle reflection.
75+ result = -log10 (one / x);
76+ }
77+ else if (x > one)
78+ {
79+ // The algorithm for base-10 logarithm is based on Chapter 5,
80+ // pages 35-36 of Cody and Waite, "Software Manual for the
81+ // Elementary Functions", Prentice Hall, 1980.
82+
83+ // In this implementation, however, we use 2s (as for natural
84+ // logarithm in Cody and Waite) even though we are computing
85+ // the base-10 logarithm.
86+
87+ T g { gn, -std::numeric_limits<T>::digits10 };
88+
89+ exp10val += std::numeric_limits<T>::digits10;
90+
91+ int reduce_sqrt2 { };
92+
93+ while (g < numbers::inv_sqrt2_v<T>)
94+ {
95+ g += g;
96+
97+ ++reduce_sqrt2;
98+ }
99+
100+ const T s { (g - one) / (g + one) };
101+ const T z { s + s };
102+ const T zsq { z * z };
103+
104+ result = z * fma (detail::log_series_expansion (zsq), zsq, one);
105+
106+ result /= numbers::ln10_v<T>;
107+
108+ for (int i = 0 ; i < reduce_sqrt2; ++i)
109+ {
110+ result -= numbers::log10_2_v<T>;
111+ }
112+
113+ result += static_cast <T>(exp10val);
114+ }
115+ else
116+ {
117+ result = zero;
118+ }
119+ }
120+ }
121+
122+ return result;
41123}
42124
43- } // namespace detail
125+ } // namespace detail
44126
45127BOOST_DECIMAL_EXPORT template <typename T>
46128constexpr auto log10 (T x) noexcept
0 commit comments