Skip to content

Commit 84785bf

Browse files
ckormanyosckormanyos
authored
Merge pull request #532 from cppalliance/tanh128_2
Fix #531 and fix #533 via tanh 128-bit
2 parents 7f74bc0 + 7fdcf98 commit 84785bf

File tree

17 files changed

+612
-169
lines changed

17 files changed

+612
-169
lines changed

‎.github/workflows/codecov.yml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -27,7 +27,7 @@ jobs:
2727
standard: [ c++20 ]
2828
compiler: [ g++ ]
2929
steps:
30-
- uses: actions/checkout@v3
30+
- uses: actions/checkout@v4
3131
with:
3232
fetch-depth: '0'
3333
- name: update-tools

‎.github/workflows/fuzz.yml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -115,7 +115,7 @@ jobs:
115115
fi
116116
fi
117117
git config --global pack.threads 0
118-
- uses: actions/checkout@v3
118+
- uses: actions/checkout@v4
119119

120120
- name: Install packages
121121
if: matrix.install

‎.github/workflows/metal.yml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -22,7 +22,7 @@ jobs:
2222
run:
2323
shell: bash
2424
steps:
25-
- uses: actions/checkout@v3
25+
- uses: actions/checkout@v4
2626
with:
2727
fetch-depth: '0'
2828
- name: update-tools

‎.github/workflows/tidy.yml‎

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -27,7 +27,7 @@ jobs:
2727
standard: [ c++14 ]
2828
compiler: [ g++ ]
2929
steps:
30-
- uses: actions/checkout@v3
30+
- uses: actions/checkout@v4
3131
with:
3232
fetch-depth: '0'
3333
- name: update-tools

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

Lines changed: 1 addition & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -68,11 +68,10 @@ constexpr auto cosh_impl(T x) noexcept
6868
else
6969
{
7070
const auto exp_pos_val = exp(x);
71-
const auto exp_neg_val = one / exp_pos_val;
7271

7372
constexpr T two { 2, 0 };
7473

75-
result = (exp_pos_val + exp_neg_val) / two;
74+
result = (exp_pos_val + (one / exp_pos_val)) / two;
7675
}
7776
}
7877
}

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

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -76,7 +76,7 @@ constexpr auto exp_impl(T x) noexcept
7676
x -= numbers::ln2_v<T> * nf2;
7777
}
7878

79-
result = detail::exp_pade_appxroximant(x);
79+
result = detail::exp_pade_appxroximant_or_series(x);
8080

8181
if (nf2 > 0)
8282
{

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

Lines changed: 73 additions & 44 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,7 @@
77
#define BOOST_DECIMAL_DETAIL_CMATH_IMPL_EXP_IMPL_HPP
88

99
#include <boost/decimal/detail/concepts.hpp>
10+
#include <boost/decimal/detail/cmath/impl/taylor_series_result.hpp>
1011

1112
#ifndef BOOST_DECIMAL_BUILD_MODULE
1213
#include <array>
@@ -18,32 +19,76 @@ namespace boost {
1819
namespace decimal {
1920
namespace detail {
2021

22+
namespace exp_detail {
23+
24+
template <bool b>
25+
struct exp_table_imp
26+
{
27+
private:
28+
using d128_coeffs_t = std::array<decimal128, 17>;
29+
30+
public:
31+
static constexpr d128_coeffs_t d128_coeffs =
32+
{{
33+
// Series[Exp[x] - 1, {x, 0, 18}]
34+
// (1), // * x
35+
::boost::decimal::decimal128 { 5, -1 }, // * x^2
36+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(90350181040458), UINT64_C(12964998083131386532) }, -34 },
37+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(225875452601146), UINT64_C(13965751134118914724) }, -35 },
38+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(451750905202293), UINT64_C(9484758194528277842) }, -36 },
39+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(75291817533715), UINT64_C(10804165069276155440) }, -36 },
40+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(107559739333879), UINT64_C(7528774067376128516) }, -37 },
41+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(134449674167349), UINT64_C(4799281565792772746) }, -38 },
42+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(149388526852610), UINT64_C(5332535073103080820) }, -39 },
43+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(149388526852610), UINT64_C(5332535073103080820) }, -40 },
44+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(135807751684191), UINT64_C(3170782423392841514) }, -41 },
45+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(113173126403492), UINT64_C(11865690723015477068) }, -42 },
46+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(87056251079609), UINT64_C(13384395342406417346) }, -43 },
47+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(62183036485435), UINT64_C(9560282387433155251) }, -44 },
48+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(414553576569570), UINT64_C(2246069003855862950) }, -46 },
49+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(259095985355981), UINT64_C(6015479145837302244) }, -47 },
50+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(152409403150577), UINT64_C(4623619737181327888) }, -48 },
51+
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(84671890639209), UINT64_C(10767230553416093986) }, -49 },
52+
}};
53+
};
54+
55+
#if !(defined(__cpp_inline_variables) && __cpp_inline_variables >= 201606L) && (!defined(_MSC_VER) || _MSC_VER != 1900)
56+
57+
template <bool b>
58+
constexpr typename exp_table_imp<b>::d128_coeffs_t exp_table_imp<b>::d128_coeffs;
59+
60+
#endif
61+
62+
} //namespace exp_detail
63+
64+
using exp_table = exp_detail::exp_table_imp<true>;
65+
2166
template <BOOST_DECIMAL_DECIMAL_FLOATING_TYPE T>
22-
constexpr auto exp_pade_appxroximant(T x) noexcept;
67+
constexpr auto exp_pade_appxroximant_or_series(T x) noexcept;
2368

2469
template <>
25-
constexpr auto exp_pade_appxroximant<decimal32>(decimal32 x) noexcept
70+
constexpr auto exp_pade_appxroximant_or_series<decimal32>(decimal32 x) noexcept
2671
{
27-
// TODO: Chris: At 32-bit, reduce the number of coefficients in the Pade appxorimant of the exp() function.
28-
2972
using local_float_t = decimal32;
3073

31-
// PadeApproximant[Exp[x] - 1, {x, 0, {6, 6}}]
74+
// PadeApproximant[Exp[x] - 1, {x, 0, {3, 4}}]
3275
// FullSimplify[%]
33-
// (84 x (7920 + 240 x^2 + x^4))
34-
// / (665280 + x (-332640 + x (75600 + x (-10080 + x (840 + (-42 + x) x)))))
76+
// (40 x (42 + x^2))
77+
// / (1680 + x (-840 + x (180 + (-20 + x) x)))
3578

3679
const auto x2 = x * x;
3780

3881
// 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))));
82+
const local_float_t top { local_float_t { UINT8_C(40), 0 } * x * (local_float_t { UINT8_C(42), 0 } + x2) };
83+
const local_float_t bot { local_float_t { UINT16_C(1680), 0 } + x * (local_float_t { INT16_C(-840), 0 } + x * (local_float_t { UINT8_C(180), 0 } + (local_float_t { INT8_C(-20), 0 } + x) * x)) };
4184

42-
return local_float_t { 1 } + (top / bot);
85+
constexpr local_float_t one { 1 };
86+
87+
return one + (top / bot);
4388
}
4489

4590
template <>
46-
constexpr auto exp_pade_appxroximant<decimal64>(decimal64 x) noexcept
91+
constexpr auto exp_pade_appxroximant_or_series<decimal64>(decimal64 x) noexcept
4792
{
4893
using local_float_t = decimal64;
4994

@@ -58,54 +103,38 @@ constexpr auto exp_pade_appxroximant<decimal64>(decimal64 x) noexcept
58103
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);
59104
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))));
60105

61-
return local_float_t { 1 } + (top / bot);
106+
constexpr local_float_t one { 1 };
107+
108+
return one + (top / bot);
62109
}
63110

64111
template <>
65-
constexpr auto exp_pade_appxroximant<decimal128>(decimal128 x) noexcept
112+
constexpr auto exp_pade_appxroximant_or_series<decimal128>(decimal128 x) noexcept
66113
{
67114
// 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}].
115+
// Use argument scaling in combination with a Taylor series expansion to order-18.
71116

72117
using local_float_t = decimal128;
73118

74119
// Rescale the argument even further (and note the three squarings below).
75120
x /= 8;
76121

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) };
122+
constexpr local_float_t one { 1, 0 };
103123

104-
result *= result;
124+
// Note: The Taylor series expansion begins with coefficients of order-2.
125+
// So we need to multiply by x^2 and add the two skipped terms (1 + x).
126+
127+
local_float_t
128+
result
129+
{
130+
one + (x * (one + (x * taylor_series_result(x, exp_table::d128_coeffs))))
131+
};
132+
133+
// Scale up with three squarings in order to obtain the result.
105134
result *= result;
106135
result *= result;
107136

108-
return result;
137+
return result *= result;
109138
}
110139

111140
} //namespace detail

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

Lines changed: 1 addition & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -59,7 +59,7 @@ struct sinh_table_imp
5959

6060
static constexpr d128_coeffs_t d128_coeffs =
6161
{{
62-
// Series[Cosh[x], {x, 0, 34}]
62+
// Series[Sinh[x], {x, 0, 34}]
6363
// (1), // * x
6464
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(90350181040458), UINT64_C(12964998083131386532) }, -34 }, // * x^3
6565
::boost::decimal::decimal128 { boost::decimal::detail::uint128 { UINT64_C(451750905202293), UINT64_C(9484758194528277842) }, -36 }, // * x^5

0 commit comments

Comments
 (0)