Skip to content

Commit a145fdc

Browse files
ckormanyosckormanyos
committed
Impl expm1 64/128 and other small repairs
1 parent 43a9bf1 commit a145fdc

File tree

6 files changed

+323
-86
lines changed

6 files changed

+323
-86
lines changed

‎include/boost/decimal/detail/cmath/expm1.hpp‎

Lines changed: 4 additions & 62 deletions
Original file line numberDiff line numberDiff line change
@@ -7,13 +7,14 @@
77
#define BOOST_DECIMAL_DETAIL_CMATH_EXPM1_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/expm1_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
16-
#include <array>
17+
#include <cmath>
1718
#include <type_traits>
1819
#endif
1920

@@ -70,66 +71,7 @@ constexpr auto expm1_impl(T x) noexcept
7071
}
7172
else
7273
{
73-
// Specifically derive a polynomial expansion for Exp[x] - 1 for this work.
74-
// Table[{x, Exp[x] - 1}, {x, -Log[2], Log[2], 1/60}]
75-
// N[%, 48]
76-
// Fit[%, {x, x^2, x^3, x^4, x^5, x^6, x^7, x^8, x^9, x^10, x^11, x^12, x^13, x^14}, x]
77-
78-
// 0.1000000000000000003213692169066381945529176657E+01 x
79-
// + 0.4999999999999999998389405741198241227273662223E+00 x^2
80-
// + 0.1666666666666664035765593562709186076539985328E+00 x^3
81-
// + 0.4166666666666666934614928838666442575683452206E-01 x^4
82-
// + 0.8333333333339521841328202617206868796855583809E-02 x^5
83-
// + 0.1388888888888953513176946682731620625302469979E-02 x^6
84-
// + 0.1984126983488689186859793276462049624439889135E-03 x^7
85-
// + 0.2480158730001499149369647648735612017495156006E-04 x^8
86-
// + 0.2755732258782898252481007286813761544775538366E-05 x^9
87-
// + 0.2755732043147979013276287368071846972098889744E-06 x^10
88-
// + 0.2505116286861719378770371641094067075234027345E-07 x^11
89-
// + 0.2087632598463662328337672597832718168295672334E-08 x^12
90-
// + 0.1619385892296180390338553597911165126625722485E-09 x^13
91-
// + 0.1154399218598221557485183048765404442959841646E-10 x^14
92-
93-
using coefficient_array_type = std::array<T, static_cast<std::size_t>(UINT8_C(14))>;
94-
95-
#if (defined(__clang__) && (__clang__ < 6))
96-
# pragma clang diagnostic push
97-
# pragma clang diagnostic ignored "-Wmissing-braces"
98-
#endif
99-
100-
constexpr auto coefficient_table =
101-
coefficient_array_type
102-
{
103-
T { UINT64_C(100000000000000000), -17 - 0 }, // * x
104-
T { UINT64_C(500000000000000000), -18 - 0 }, // * x^2
105-
T { UINT64_C(166666666666666404), -18 - 0 }, // * x^3
106-
T { UINT64_C(416666666666666693), -18 - 1 }, // * x^4
107-
T { UINT64_C(833333333333952184), -18 - 2 }, // * x^5
108-
T { UINT64_C(138888888888895351), -18 - 2 }, // * x^6
109-
T { UINT64_C(198412698348868919), -18 - 3 }, // * x^7
110-
T { UINT64_C(248015873000149915), -18 - 4 }, // * x^8
111-
T { UINT64_C(275573225878289825), -18 - 5 }, // * x^9
112-
T { UINT64_C(275573204314797901), -18 - 6 }, // * x^10
113-
T { UINT64_C(250511628686171938), -18 - 7 }, // * x^11
114-
T { UINT64_C(208763259846366233), -18 - 8 }, // * x^12
115-
T { UINT64_C(161938589229618039), -18 - 9 }, // * x^13
116-
T { UINT64_C(115439921859822156), -18 - 10 } // * x^14
117-
};
118-
119-
#if (defined(__clang__) && (__clang__ < 6))
120-
# pragma clang diagnostic pop
121-
#endif
122-
123-
auto rit = coefficient_table.crbegin() + static_cast<std::size_t>((sizeof(T) == 4U) ? 5U : 0U);
124-
125-
result = *rit;
126-
127-
while(rit != coefficient_table.crend())
128-
{
129-
result = fma(result, x, *rit++);
130-
}
131-
132-
result *= x;
74+
result = x * detail::expm1_series_expansion(x);
13375
}
13476
}
13577

‎include/boost/decimal/detail/cmath/ilogb.hpp‎

Lines changed: 15 additions & 11 deletions
Original file line numberDiff line numberDiff line change
@@ -22,26 +22,30 @@ BOOST_DECIMAL_EXPORT template <typename T>
2222
constexpr auto ilogb(T d) noexcept
2323
BOOST_DECIMAL_REQUIRES_RETURN(detail::is_decimal_floating_point_v, T, int)
2424
{
25-
const auto fpc_d = fpclassify(d);
25+
const auto fpc = fpclassify(d);
2626

27-
if (fpc_d == FP_ZERO)
27+
int result { };
28+
29+
if (fpc == FP_ZERO)
2830
{
29-
return FP_ILOGB0;
31+
result = static_cast<int>(FP_ILOGB0);
3032
}
31-
else if (fpc_d == FP_INFINITE)
33+
else if (fpc == FP_INFINITE)
3234
{
33-
return INT_MAX;
35+
result = static_cast<int>(INT_MAX);
3436
}
35-
else if (fpc_d == FP_NAN)
37+
else if (fpc == FP_NAN)
3638
{
37-
return FP_ILOGBNAN;
39+
result = static_cast<int>(FP_ILOGBNAN);
3840
}
41+
else
42+
{
43+
const auto offset = detail::num_digits(d.full_significand()) - 1;
3944

40-
const auto offset = detail::num_digits(d.full_significand()) - 1;
41-
42-
const auto expval = static_cast<int>(static_cast<int>(d.unbiased_exponent()) + offset);
45+
result = static_cast<int>(static_cast<int>(d.unbiased_exponent()) + offset);
46+
}
4347

44-
return expval;
48+
return result;
4549
}
4650

4751
} // namespace decimal

‎include/boost/decimal/detail/cmath/impl/expm1_impl.hpp‎

Lines changed: 158 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,158 @@
1+
// Copyright 2023 - 2024 Matt Borland
2+
// Copyright 2023 - 2024 Christopher Kormanyos
3+
// Distributed under the Boost Software License, Version 1.0.
4+
// https://www.boost.org/LICENSE_1_0.txt
5+
6+
#ifndef BOOST_DECIMAL_DETAIL_CMATH_IMPL_EXPM1_IMPL_HPP
7+
#define BOOST_DECIMAL_DETAIL_CMATH_IMPL_EXPM1_IMPL_HPP
8+
9+
#include <boost/decimal/detail/concepts.hpp>
10+
#include <boost/decimal/detail/cmath/impl/taylor_series_result.hpp>
11+
12+
#ifndef BOOST_DECIMAL_BUILD_MODULE
13+
#include <array>
14+
#include <cstddef>
15+
#include <cstdint>
16+
#endif
17+
18+
namespace boost {
19+
namespace decimal {
20+
namespace detail {
21+
22+
namespace expm1_detail {
23+
24+
template <bool b>
25+
struct expm1_table_imp
26+
{
27+
private:
28+
using d32_coeffs_t = std::array<decimal32, 10>;
29+
using d64_coeffs_t = std::array<decimal64, 14>;
30+
using d128_coeffs_t = std::array<decimal128, 32>;
31+
32+
public:
33+
static constexpr d32_coeffs_t d32_coeffs =
34+
{{
35+
// Specifically derive a polynomial expansion for Exp[x] - 1 for this work.
36+
// Table[{x, Exp[x] - 1}, {x, -Log[2], Log[2], 1/60}]
37+
// N[%, 48]
38+
// Fit[%, {x, x^2, x^3, x^4, x^5, x^6, x^7, x^8, x^9, x^10}, x]
39+
40+
::boost::decimal::decimal32 { UINT64_C(1000000000005449334), - 19 + 1 }, // * x
41+
::boost::decimal::decimal32 { UINT64_C(5000000000003881336), - 19 - 0 }, // * x^2
42+
::boost::decimal::decimal32 { UINT64_C(1666666664242981149), - 19 - 0 }, // * x^3
43+
::boost::decimal::decimal32 { UINT64_C(4166666665026072773), - 19 - 1 }, // * x^4
44+
::boost::decimal::decimal32 { UINT64_C(8333336317448167991), - 19 - 2 }, // * x^5
45+
::boost::decimal::decimal32 { UINT64_C(1388889096793935619), - 19 - 2 }, // * x^6
46+
::boost::decimal::decimal32 { UINT64_C(1983978347911205530), - 19 - 3 }, // * x^7
47+
::boost::decimal::decimal32 { UINT64_C(2480049494648544583), - 19 - 4 }, // * x^8
48+
::boost::decimal::decimal32 { UINT64_C(2787876201220259352), - 19 - 5 }, // * x^9
49+
::boost::decimal::decimal32 { UINT64_C(2780855729673643225), - 19 - 6 }, // * x^10
50+
}};
51+
52+
static constexpr d64_coeffs_t d64_coeffs =
53+
{{
54+
// Specifically derive a polynomial expansion for Exp[x] - 1 for this work.
55+
// Table[{x, Exp[x] - 1}, {x, -Log[2], Log[2], 1/60}]
56+
// N[%, 48]
57+
// Fit[%, {x, x^2, x^3, x^4, x^5, x^6, x^7, x^8, x^9, x^10, x^11, x^12, x^13, x^14}, x]
58+
59+
::boost::decimal::decimal64 { UINT64_C(1000000000000000003), - 19 + 1 }, // * x
60+
::boost::decimal::decimal64 { UINT64_C(4999999999999999998), - 19 - 0 }, // * x^2
61+
::boost::decimal::decimal64 { UINT64_C(1666666666666664035), - 19 - 0 }, // * x^3
62+
::boost::decimal::decimal64 { UINT64_C(4166666666666666934), - 19 - 1 }, // * x^4
63+
::boost::decimal::decimal64 { UINT64_C(8333333333339521841), - 19 - 2 }, // * x^5
64+
::boost::decimal::decimal64 { UINT64_C(1388888888888953513), - 19 - 2 }, // * x^6
65+
::boost::decimal::decimal64 { UINT64_C(1984126983488689186), - 19 - 3 }, // * x^7
66+
::boost::decimal::decimal64 { UINT64_C(2480158730001499149), - 19 - 4 }, // * x^8
67+
::boost::decimal::decimal64 { UINT64_C(2755732258782898252), - 19 - 5 }, // * x^9
68+
::boost::decimal::decimal64 { UINT64_C(2755732043147979013), - 19 - 6 }, // * x^10
69+
::boost::decimal::decimal64 { UINT64_C(2505116286861719378), - 19 - 7 }, // * x^11
70+
::boost::decimal::decimal64 { UINT64_C(2087632598463662328), - 19 - 8 }, // * x^12
71+
::boost::decimal::decimal64 { UINT64_C(1619385892296180390), - 19 - 9 }, // * x^13
72+
::boost::decimal::decimal64 { UINT64_C(1154399218598221557), - 19 - 10 } // * x^14
73+
}};
74+
75+
static constexpr d128_coeffs_t d128_coeffs =
76+
{{
77+
// Specifically derive a polynomial expansion for Exp[x] - 1 for this work.
78+
// Table[{x, Exp[x] - 1}, {x, -Log[2], Log[2], 1/60}]
79+
// N[%, 48]
80+
// Fit[%, {x, x^2, x^3, x^4, x^5, x^6, x^7, x^8, x^9, x^10, x^11, x^12, x^13, x^14, x^15, x^16, x^17, x^18, x^19, x^20, x^21, x^22, x^23, x^24, x^25, x^26, x^27, x^28, x^29, x^30, x^31, x^32 }, x]
81+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(54210108624275), UINT64_C(4089650035136921600) }, -33 }, // * x
82+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(271050543121376), UINT64_C(2001506101975056384) }, -34 }, // * x^2
83+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(90350181040458), UINT64_C(12964998083131386532) }, -34 }, // * x^3
84+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(225875452601146), UINT64_C(13965751134118914724) }, -35 }, // * x^4
85+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(451750905202293), UINT64_C(9484758194528277842) }, -36 }, // * x^5
86+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(75291817533715), UINT64_C(10804165069276155440) }, -36 }, // * x^6
87+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(107559739333879), UINT64_C(7528774067376128516) }, -37 }, // * x^7
88+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(134449674167349), UINT64_C(4799281565792772746) }, -38 }, // * x^8
89+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(149388526852610), UINT64_C(5332535073103080820) }, -39 }, // * x^9
90+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(149388526852610), UINT64_C(5332535073103080820) }, -40 }, // * x^10
91+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(135807751684191), UINT64_C(3170782423392841514) }, -41 }, // * x^11
92+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(113173126403492), UINT64_C(11865690723015477068) }, -42 }, // * x^12
93+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(87056251079609), UINT64_C(13384395342406416636) }, -43 }, // * x^13
94+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(62183036485435), UINT64_C(9560282387433156335) }, -44 }, // * x^14
95+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(414553576569570), UINT64_C(2246069003862680020) }, -46 }, // * x^15
96+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(259095985355981), UINT64_C(6015479145828949264) }, -47 }, // * x^16
97+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(152409403150577), UINT64_C(4623619732418095578) }, -48 }, // * x^17
98+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(84671890639209), UINT64_C(10767230558026320466) }, -49 }, // * x^18
99+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(445641529680050), UINT64_C(8125595620937745600) }, -51 }, // * x^19
100+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(222820764840025), UINT64_C(4062767274683195140) }, -52 }, // * x^20
101+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(106105126114297), UINT64_C(13344759429965740488) }, -53 }, // * x^21
102+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(482296027792262), UINT64_C(7088674266265745598) }, -55 }, // * x^22
103+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(209693925127072), UINT64_C(336105452763225878) }, -56 }, // * x^23
104+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(87372468802945), UINT64_C(10013088901203012320) }, -57 }, // * x^24
105+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(349489875208886), UINT64_C(9445768661182748344) }, -59 }, // * x^25
106+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(134419182774415), UINT64_C(9680981560342232810) }, -60 }, // * x^26
107+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(497848829278818), UINT64_C(16288994997110182382) }, -62 }, // * x^27
108+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(177803151475355), UINT64_C(16680206430774781810) }, -63 }, // * x^28
109+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(61311025561137), UINT64_C(7837795588749518446) }, -64 }, // * x^29
110+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(204371229207757), UINT64_C(18366861741830034248) }, -66 }, // * x^30
111+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(66162682638108), UINT64_C(6755035083974089930) }, -67 }, // * x^31
112+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(206436477688751), UINT64_C(15666750779045089894) }, -69 }, // * x^32
113+
}};
114+
};
115+
116+
#if !(defined(__cpp_inline_variables) && __cpp_inline_variables >= 201606L) && (!defined(_MSC_VER) || _MSC_VER != 1900)
117+
118+
template <bool b>
119+
constexpr typename expm1_table_imp<b>::d32_coeffs_t expm1_table_imp<b>::d32_coeffs;
120+
121+
template <bool b>
122+
constexpr typename expm1_table_imp<b>::d64_coeffs_t expm1_table_imp<b>::d64_coeffs;
123+
124+
template <bool b>
125+
constexpr typename expm1_table_imp<b>::d128_coeffs_t expm1_table_imp<b>::d128_coeffs;
126+
127+
#endif
128+
129+
} //namespace expm1_detail
130+
131+
using expm1_table = expm1_detail::expm1_table_imp<true>;
132+
133+
template <BOOST_DECIMAL_DECIMAL_FLOATING_TYPE T>
134+
constexpr auto expm1_series_expansion(T x) noexcept;
135+
136+
template <>
137+
constexpr auto expm1_series_expansion<decimal32>(decimal32 x) noexcept
138+
{
139+
return taylor_series_result(x, expm1_table::d32_coeffs);
140+
}
141+
142+
template <>
143+
constexpr auto expm1_series_expansion<decimal64>(decimal64 x) noexcept
144+
{
145+
return taylor_series_result(x, expm1_table::d64_coeffs);
146+
}
147+
148+
template <>
149+
constexpr auto expm1_series_expansion<decimal128>(decimal128 x) noexcept
150+
{
151+
return taylor_series_result(x, expm1_table::d128_coeffs);
152+
}
153+
154+
} //namespace detail
155+
} //namespace decimal
156+
} //namespace boost
157+
158+
#endif //BOOST_DECIMAL_DETAIL_CMATH_IMPL_EXPM1_IMPL_HPP

‎include/boost/decimal/detail/cmath/log.hpp‎

Lines changed: 5 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -52,6 +52,11 @@ constexpr auto log_impl(T x) noexcept
5252
}
5353
else if(x > one)
5454
{
55+
// Use the implementation of log10 in order to compute the natural
56+
// logarithm. The base of the boost::decimal library is, in fact,
57+
// base-10. And so, somewhat uncommonly, the fastest and most accurate
58+
// logarithm in this system is log10 in base-10.
59+
5560
result = log10(x) * numbers::ln10_v<T>;
5661
}
5762
else

0 commit comments

Comments
 (0)