Skip to content

Commit de795d4

Browse files
ckormanyosckormanyos
authored
Merge pull request #527 from cppalliance/cosh128
Fix #522 via impl of cosh/exp at 128 bits and a few tests
2 parents a199bd0 + 49486d5 commit de795d4

File tree

7 files changed

+390
-59
lines changed

7 files changed

+390
-59
lines changed

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

Lines changed: 3 additions & 41 deletions
Original file line numberDiff line numberDiff line change
@@ -7,9 +7,10 @@
77
#define BOOST_DECIMAL_DETAIL_CMATH_COSH_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/cosh_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
@@ -58,48 +59,9 @@ constexpr auto cosh_impl(T x) noexcept
5859
{
5960
if (x < one)
6061
{
61-
using coefficient_array_type = std::array<T, static_cast<std::size_t>(UINT8_C(9))>;
62-
63-
#if (defined(__clang__) && (__clang__ < 6))
64-
# pragma clang diagnostic push
65-
# pragma clang diagnostic ignored "-Wmissing-braces"
66-
#endif
67-
68-
constexpr auto coefficient_table =
69-
coefficient_array_type
70-
{
71-
// Series[Cosh[x], {x, 0, 18}]
72-
// (1), // * 1
73-
T { 5, -1 }, // * x^2
74-
T { UINT64_C(416666666666666667), - 18 - 1 }, // * x^4
75-
T { UINT64_C(138888888888888889), - 18 - 2 }, // * x^6
76-
T { UINT64_C(248015873015873016), - 18 - 4 }, // * x^8
77-
T { UINT64_C(275573192239858907), - 18 - 6 }, // * x^10
78-
T { UINT64_C(208767569878680990), - 18 - 8 }, // * x^12
79-
T { UINT64_C(114707455977297247), - 18 - 10 }, // * x^14
80-
T { UINT64_C(477947733238738530), - 18 - 13 }, // * x^16
81-
T { UINT64_C(156192069685862265), - 18 - 15 } // * x^18
82-
};
83-
84-
#if (defined(__clang__) && (__clang__ < 6))
85-
# pragma clang diagnostic pop
86-
#endif
87-
88-
auto rit =
89-
coefficient_table.crbegin()
90-
+ static_cast<std::size_t>
91-
(
92-
(sizeof(T) == static_cast<std::size_t>(UINT8_C(4))) ? 4U : 0U
93-
);
94-
95-
result = *rit;
96-
9762
const auto xsq = x * x;
9863

99-
while(rit != coefficient_table.crend())
100-
{
101-
result = fma(result, xsq, *rit++);
102-
}
64+
result = detail::cosh_series_expansion(xsq);
10365

10466
result = fma(result, xsq, one);
10567
}

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

Lines changed: 2 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@
77
#define BOOST_DECIMAL_DETAIL_CMATH_EXP_HPP
88

99
#include <boost/decimal/fwd.hpp> // NOLINT(llvm-include-order)
10+
#include <boost/decimal/detail/cmath/impl/exp_impl.hpp>
1011
#include <boost/decimal/detail/cmath/impl/pow_impl.hpp>
1112
#include <boost/decimal/detail/type_traits.hpp>
1213
#include <boost/decimal/detail/concepts.hpp>
@@ -75,18 +76,7 @@ constexpr auto exp_impl(T x) noexcept
7576
x -= numbers::ln2_v<T> * nf2;
7677
}
7778

78-
// PadeApproximant[Exp[x] - 1, {x, 0, {6, 6}}]
79-
// FullSimplify[%]
80-
// (84 x (7920 + 240 x^2 + x^4))
81-
// / (665280 + x (-332640 + x (75600 + x (-10080 + x (840 + (-42 + x) x)))))
82-
83-
const auto x2 = x * x;
84-
85-
// Use the small-argument Pade approximation having coefficients shown above.
86-
const T top = T { UINT8_C(84), 0 } * x * ( T { UINT16_C(7920), 0 } + ( T { UINT8_C(240), 0 } + x2) * x2);
87-
const T bot = T { UINT32_C(665280), 0 } + x * (T { INT32_C(-332640), 0 } + x * (T { UINT32_C(75600), 0 } + x * (T { INT16_C(-10080), 0 } + x * (T { UINT16_C(840), 0 } + (T { INT8_C(-42), 0 } + x) * x))));
88-
89-
result = one + (top / bot);
79+
result = detail::exp_pade_appxroximant(x);
9080

9181
if (nf2 > 0)
9282
{

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

Lines changed: 126 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,126 @@
1+
// Copyright 2024 Matt Borland
2+
// Copyright 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_COSH_IMPL_HPP
7+
#define BOOST_DECIMAL_DETAIL_CMATH_IMPL_COSH_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 cosh_detail {
23+
24+
template <bool b>
25+
struct cosh_table_imp
26+
{
27+
private:
28+
using d32_coeffs_t = std::array<decimal32, 6>;
29+
using d64_coeffs_t = std::array<decimal64, 9>;
30+
using d128_coeffs_t = std::array<decimal128, 17>;
31+
32+
public:
33+
static constexpr d32_coeffs_t d32_coeffs =
34+
{{
35+
// Series[Cosh[x], {x, 0, 10}]
36+
// (1), // * 1
37+
::boost::decimal::decimal32 { 5, -1 }, // * x^2
38+
::boost::decimal::decimal32 { UINT64_C(4166666666666666667), - 19 - 1 }, // * x^6
39+
::boost::decimal::decimal32 { UINT64_C(1388888888888888889), - 19 - 2 }, // * x^8
40+
::boost::decimal::decimal32 { UINT64_C(2480158730158730159), - 19 - 4 }, // * x^10
41+
::boost::decimal::decimal32 { UINT64_C(2755731922398589065), - 19 - 6 }, // * x^12
42+
::boost::decimal::decimal32 { UINT64_C(2087675698786809898), - 19 - 8 }, // * x^12
43+
}};
44+
45+
static constexpr d64_coeffs_t d64_coeffs =
46+
{{
47+
// Series[Cosh[x], {x, 0, 18}]
48+
// (1), // * 1
49+
::boost::decimal::decimal64 { 5, -1 }, // * x^2
50+
::boost::decimal::decimal64 { UINT64_C(4166666666666666667), - 19 - 1 }, // * x^4
51+
::boost::decimal::decimal64 { UINT64_C(1388888888888888889), - 19 - 2 }, // * x^6
52+
::boost::decimal::decimal64 { UINT64_C(2480158730158730159), - 19 - 4 }, // * x^8
53+
::boost::decimal::decimal64 { UINT64_C(2755731922398589065), - 19 - 6 }, // * x^10
54+
::boost::decimal::decimal64 { UINT64_C(2087675698786809898), - 19 - 8 }, // * x^12
55+
::boost::decimal::decimal64 { UINT64_C(1147074559772972471), - 19 - 10 }, // * x^14
56+
::boost::decimal::decimal64 { UINT64_C(4779477332387385297), - 19 - 13 }, // * x^16
57+
::boost::decimal::decimal64 { UINT64_C(1561920696858622646), - 19 - 15 } // * x^18
58+
}};
59+
60+
static constexpr d128_coeffs_t d128_coeffs =
61+
{{
62+
// Series[Cosh[x], {x, 0, 34}]
63+
// (1), // * 1
64+
::boost::decimal::decimal128 { 5, -1 }, // * x^2
65+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(225875452601146), UINT64_C(13965751134118914724) }, -35 }, // * x^4
66+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(75291817533715), UINT64_C(10804165069276155440) }, -36 }, // * x^6
67+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(134449674167349), UINT64_C(4799281565792772746) }, -38 }, // * x^8
68+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(149388526852610), UINT64_C(5332535073103080820) }, -40 }, // * x^10
69+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(113173126403492), UINT64_C(11865690723015477068) }, -42 }, // * x^12
70+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(62183036485435), UINT64_C(9560282387433155251) }, -44 }, // * x^14
71+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(259095985355981), UINT64_C(6015479145837302244) }, -47 }, // * x^16
72+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(84671890639209), UINT64_C(10767230553416093986) }, -49 }, // * x^18
73+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(222820764840025), UINT64_C(4062785569898205740) }, -52 }, // * x^20
74+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(482296027792262), UINT64_C(7037075391028107068) }, -55 }, // * x^22
75+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(87372468802946), UINT64_C(1542176615384940434) }, -57 }, // * x^24
76+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(134419182773763), UINT64_C(3791559721646796942) }, -60 }, // * x^26
77+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(177803151817147), UINT64_C(1794430560736952558) }, -63 }, // * x^28
78+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(204371438870284), UINT64_C(366311534299067156) }, -66 }, // * x^30
79+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(206019595635366), UINT64_C(17625897212400736954) }, -69 }, // * x^32
80+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(183618177928134), UINT64_C(9987905770721758456) }, -72 }, // * x^34
81+
}};
82+
};
83+
84+
#if !(defined(__cpp_inline_variables) && __cpp_inline_variables >= 201606L) && (!defined(_MSC_VER) || _MSC_VER != 1900)
85+
86+
template <bool b>
87+
constexpr typename cosh_table_imp<b>::d32_coeffs_t cosh_table_imp<b>::d32_coeffs;
88+
89+
template <bool b>
90+
constexpr typename cosh_table_imp<b>::d64_coeffs_t cosh_table_imp<b>::d64_coeffs;
91+
92+
template <bool b>
93+
constexpr typename cosh_table_imp<b>::d128_coeffs_t cosh_table_imp<b>::d128_coeffs;
94+
95+
#endif
96+
97+
} //namespace cosh_detail
98+
99+
using cosh_table = cosh_detail::cosh_table_imp<true>;
100+
101+
template <BOOST_DECIMAL_DECIMAL_FLOATING_TYPE T>
102+
constexpr auto cosh_series_expansion(T z) noexcept;
103+
104+
template <>
105+
constexpr auto cosh_series_expansion<decimal32>(decimal32 z2) noexcept
106+
{
107+
return taylor_series_result(z2, cosh_table::d32_coeffs);
108+
}
109+
110+
template <>
111+
constexpr auto cosh_series_expansion<decimal64>(decimal64 z2) noexcept
112+
{
113+
return taylor_series_result(z2, cosh_table::d64_coeffs);
114+
}
115+
116+
template <>
117+
constexpr auto cosh_series_expansion<decimal128>(decimal128 z2) noexcept
118+
{
119+
return taylor_series_result(z2, cosh_table::d128_coeffs);
120+
}
121+
122+
} //namespace detail
123+
} //namespace decimal
124+
} //namespace boost
125+
126+
#endif //BOOST_DECIMAL_DETAIL_CMATH_IMPL_COSH_IMPL_HPP

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

Lines changed: 115 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,115 @@
1+
// Copyright 2024 Matt Borland
2+
// Copyright 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_EXP_IMPL_HPP
7+
#define BOOST_DECIMAL_DETAIL_CMATH_IMPL_EXP_IMPL_HPP
8+
9+
#include <boost/decimal/detail/concepts.hpp>
10+
11+
#ifndef BOOST_DECIMAL_BUILD_MODULE
12+
#include <array>
13+
#include <cstddef>
14+
#include <cstdint>
15+
#endif
16+
17+
namespace boost {
18+
namespace decimal {
19+
namespace detail {
20+
21+
template <BOOST_DECIMAL_DECIMAL_FLOATING_TYPE T>
22+
constexpr auto exp_pade_appxroximant(T x) noexcept;
23+
24+
template <>
25+
constexpr auto exp_pade_appxroximant<decimal32>(decimal32 x) noexcept
26+
{
27+
// TODO: Chris: At 32-bit, reduce the number of coefficients in the Pade appxorimant of the exp() function.
28+
29+
using local_float_t = decimal32;
30+
31+
// PadeApproximant[Exp[x] - 1, {x, 0, {6, 6}}]
32+
// FullSimplify[%]
33+
// (84 x (7920 + 240 x^2 + x^4))
34+
// / (665280 + x (-332640 + x (75600 + x (-10080 + x (840 + (-42 + x) x)))))
35+
36+
const auto x2 = x * x;
37+
38+
// Use the small-argument Pade approximation having coefficients shown above.
39+
const local_float_t top = local_float_t { UINT8_C(84), 0 } * x * ( local_float_t { UINT16_C(7920), 0 } + ( local_float_t { UINT8_C(240), 0 } + x2) * x2);
40+
const local_float_t bot = local_float_t { UINT32_C(665280), 0 } + x * (local_float_t { INT32_C(-332640), 0 } + x * (local_float_t { UINT32_C(75600), 0 } + x * (local_float_t { INT16_C(-10080), 0 } + x * (local_float_t { UINT16_C(840), 0 } + (local_float_t { INT8_C(-42), 0 } + x) * x))));
41+
42+
return local_float_t { 1 } + (top / bot);
43+
}
44+
45+
template <>
46+
constexpr auto exp_pade_appxroximant<decimal64>(decimal64 x) noexcept
47+
{
48+
using local_float_t = decimal64;
49+
50+
// PadeApproximant[Exp[x] - 1, {x, 0, {6, 6}}]
51+
// FullSimplify[%]
52+
// (84 x (7920 + 240 x^2 + x^4))
53+
// / (665280 + x (-332640 + x (75600 + x (-10080 + x (840 + (-42 + x) x)))))
54+
55+
const auto x2 = x * x;
56+
57+
// Use the small-argument Pade approximation having coefficients shown above.
58+
const local_float_t top = local_float_t { UINT8_C(84), 0 } * x * ( local_float_t { UINT16_C(7920), 0 } + ( local_float_t { UINT8_C(240), 0 } + x2) * x2);
59+
const local_float_t bot = local_float_t { UINT32_C(665280), 0 } + x * (local_float_t { INT32_C(-332640), 0 } + x * (local_float_t { UINT32_C(75600), 0 } + x * (local_float_t { INT16_C(-10080), 0 } + x * (local_float_t { UINT16_C(840), 0 } + (local_float_t { INT8_C(-42), 0 } + x) * x))));
60+
61+
return local_float_t { 1 } + (top / bot);
62+
}
63+
64+
template <>
65+
constexpr auto exp_pade_appxroximant<decimal128>(decimal128 x) noexcept
66+
{
67+
// Compute exp(x) - 1 for x small.
68+
69+
// Use an order-12 Pade approximation of the exponential function.
70+
// PadeApproximant[Exp[x] - 1, {x, 0, 12, 12}].
71+
72+
using local_float_t = decimal128;
73+
74+
// Rescale the argument even further (and note the three squarings below).
75+
x /= 8;
76+
77+
const local_float_t x2 = (x * x);
78+
79+
const local_float_t top = ((((( local_float_t { boost::decimal::detail::uint128 { UINT64_C(130576843339991), UINT64_C(2348781707059460614) }, -46 } * x2
80+
+ local_float_t { boost::decimal::detail::uint128 { UINT64_C(502720846858965), UINT64_C(15499169997977266440) }, -43 } ) * x2
81+
+ local_float_t { boost::decimal::detail::uint128 { UINT64_C(492264253244299), UINT64_C(6469924059228430936) }, -40 } ) * x2
82+
+ local_float_t { boost::decimal::detail::uint128 { UINT64_C(168354374609550), UINT64_C(6971973999273187690) }, -37 } ) * x2
83+
+ local_float_t { boost::decimal::detail::uint128 { UINT64_C(196413437044475), UINT64_C(8133969665818718980) }, -35 } ) * x2
84+
+ local_float_t { boost::decimal::detail::uint128 { UINT64_C(54210108624275), UINT64_C(4089650035136921600) }, -33 } )
85+
;
86+
87+
const local_float_t bot = (((((((((((( local_float_t( +boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(418515523525612), UINT64_C(10839100561497421498) }, -49 } ) * x
88+
+ local_float_t( -boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(65288421669995), UINT64_C(10397762890384506116) }, -46 } )) * x
89+
+ local_float_t( +boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(502720846858965), UINT64_C(15499169997977266440) }, -45 } )) * x
90+
+ local_float_t( -boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(251360423429482), UINT64_C(16972957035843409028) }, -43 } )) * x
91+
+ local_float_t( +boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(90489752434613), UINT64_C(15702571451232594082) }, -41 } )) * x
92+
+ local_float_t( -boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(246132126622149), UINT64_C(12458334066468991276) }, -40 } )) * x
93+
+ local_float_t( +boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(516877465906514), UINT64_C(5871083058504374896) }, -39 } )) * x
94+
+ local_float_t( -boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(84177187304775), UINT64_C(3485986999636593840) }, -37 } )) * x
95+
+ local_float_t( +boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(105221484130968), UINT64_C(18192541804827906022) }, -36 } )) * x
96+
+ local_float_t( -boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(98206718522237), UINT64_C(13290356869764135298) }, -35 } )) * x
97+
+ local_float_t( +boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(64816434224676), UINT64_C(16519268045002340975) }, -34 } )) * x
98+
+ local_float_t( -boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(271050543121376), UINT64_C(2001506101975056384) }, -34 } )) * x
99+
+ local_float_t( +boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(54210108624275), UINT64_C(4089650035136921600) }, -33 } ))
100+
;
101+
102+
local_float_t result { local_float_t { 1 } + ((x * top) / bot) };
103+
104+
result *= result;
105+
result *= result;
106+
result *= result;
107+
108+
return result;
109+
}
110+
111+
} //namespace detail
112+
} //namespace decimal
113+
} //namespace boost
114+
115+
#endif //BOOST_DECIMAL_DETAIL_CMATH_IMPL_EXP_IMPL_HPP

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

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -9,9 +9,9 @@
99
#include <boost/decimal/fwd.hpp> // NOLINT(llvm-include-order)
1010
#include <boost/decimal/detail/cmath/impl/tgamma_impl.hpp>
1111
#include <boost/decimal/detail/cmath/sin.hpp>
12+
#include <boost/decimal/detail/config.hpp>
1213
#include <boost/decimal/detail/type_traits.hpp>
1314
#include <boost/decimal/numbers.hpp>
14-
#include <boost/decimal/detail/config.hpp>
1515

1616
#ifndef BOOST_DECIMAL_BUILD_MODULE
1717
#include <iterator>

0 commit comments

Comments
 (0)