Skip to content

Commit e544227

Browse files
ckormanyosckormanyos
committed
Simplify and finish atan()
1 parent 4426c57 commit e544227

File tree

3 files changed

+186
-91
lines changed

3 files changed

+186
-91
lines changed

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

Lines changed: 6 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -51,7 +51,7 @@ constexpr auto atan_impl(T x) noexcept
5151

5252
if (x <= T { 48 })
5353
{
54-
const bool is_smallish { x <= T { 6 } };
54+
const bool is_smallish { x <= T { 3 } };
5555

5656
// The portion of the algorithm for arc-tangent regarding scaling large-valued
5757
// argument is based on Chapter 11, page 194 of Cody and Waite, "Software Manual
@@ -69,11 +69,10 @@ constexpr auto atan_impl(T x) noexcept
6969
constexpr T three_halves { 15, -1 };
7070

7171
result = (fx_arg <= std::numeric_limits<T>::epsilon()) ? fx_arg
72-
: (fx_arg <= T { 4375, -4 }) ? detail::atan_series_small (fx_arg)
73-
: (fx_arg <= T { 6875, -4 }) ? detail::atan_values<T>(0U) + detail::atan_series_small((fx_arg - half) / (one + fx_arg / 2))
74-
: (fx_arg <= T { 11875, -4 }) ? detail::atan_values<T>(1U) + detail::atan_series_small((fx_arg - one) / (fx_arg + one))
75-
: (fx_arg <= T { 24375, -4 }) ? detail::atan_values<T>(2U) + detail::atan_series_small((fx_arg - three_halves) / (one + three_halves * fx_arg))
76-
: detail::atan_series_med ( fx_arg)
72+
: (fx_arg <= T { 4375, -4 }) ? detail::atan_series (fx_arg)
73+
: (fx_arg <= T { 6875, -4 }) ? detail::atan_values<T>(0U) + detail::atan_series((fx_arg - half) / (one + fx_arg / 2))
74+
: (fx_arg <= T { 11875, -4 }) ? detail::atan_values<T>(1U) + detail::atan_series((fx_arg - one) / (fx_arg + one))
75+
: detail::atan_values<T>(2U) + detail::atan_series((fx_arg - three_halves) / (one + three_halves * fx_arg))
7776
;
7877

7978
if(!is_smallish)
@@ -85,7 +84,7 @@ constexpr auto atan_impl(T x) noexcept
8584
}
8685
else
8786
{
88-
result = my_pi_half - detail::atan_series_small(one / x);
87+
result = my_pi_half - detail::atan_series(one / x);
8988
}
9089
}
9190

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

Lines changed: 50 additions & 83 deletions
Original file line numberDiff line numberDiff line change
@@ -48,7 +48,7 @@ struct atan_table_imp
4848

4949
// 10th degree remez polynomial calculated from 0, 0.4375
5050
// Estimated max error: 2.3032664387910605e-12
51-
static constexpr std::array<::boost::decimal::decimal32, 11> d32_coeffs_small =
51+
static constexpr std::array<::boost::decimal::decimal32, 11> d32_coeffs =
5252
{{
5353
::boost::decimal::decimal32 { UINT64_C(61037779951304161), -18, true },
5454
::boost::decimal::decimal32 { UINT64_C(10723099589331457), -17 },
@@ -63,7 +63,7 @@ struct atan_table_imp
6363
::boost::decimal::decimal32 { UINT64_C(23032664387910606), -29 },
6464
}};
6565

66-
static constexpr std::array<::boost::decimal::decimal64, 11> d64_coeffs_small =
66+
static constexpr std::array<::boost::decimal::decimal64, 11> d64_coeffs =
6767
{{
6868
::boost::decimal::decimal64 { UINT64_C(61037779951304161), -18, true },
6969
::boost::decimal::decimal64 { UINT64_C(10723099589331457), -17 },
@@ -77,83 +77,15 @@ struct atan_table_imp
7777
::boost::decimal::decimal64 { UINT64_C(99999999877886492), -17 },
7878
::boost::decimal::decimal64 { UINT64_C(23032664387910606), -29 },
7979
}};
80-
81-
static constexpr std::array<::boost::decimal::decimal128, 11> d128_coeffs_small =
82-
{{
83-
::boost::decimal::decimal128 { UINT64_C(61037779951304161), -18, true },
84-
::boost::decimal::decimal128 { UINT64_C(10723099589331457), -17 },
85-
::boost::decimal::decimal128 { UINT64_C(22515613909953665), -18 },
86-
::boost::decimal::decimal128 { UINT64_C(15540713402718176), -17, true },
87-
::boost::decimal::decimal128 { UINT64_C(35999727706986597), -19 },
88-
::boost::decimal::decimal128 { UINT64_C(19938867353282852), -17 },
89-
::boost::decimal::decimal128 { UINT64_C(62252075283915644), -22 },
90-
::boost::decimal::decimal128 { UINT64_C(33333695504913247), -17, true },
91-
::boost::decimal::decimal128 { UINT64_C(10680927642397763), -24 },
92-
::boost::decimal::decimal128 { UINT64_C(99999999877886492), -17 },
93-
::boost::decimal::decimal128 { UINT64_C(23032664387910606), -29 },
94-
}};
95-
96-
// 10th degree remez polynomial from 2.4375, 6
97-
// Estimated max error: 2.6239664084435361e-9
98-
static constexpr std::array<::boost::decimal::decimal32, 11> d32_coeffs_med =
99-
{{
100-
::boost::decimal::decimal32 { UINT64_C(35895331641408534), -24, true },
101-
::boost::decimal::decimal32 { UINT64_C(1734850544519432), -21 },
102-
::boost::decimal::decimal32 { UINT64_C(38064221484608425), -21, true },
103-
::boost::decimal::decimal32 { UINT64_C(50135011697517902), -20 },
104-
::boost::decimal::decimal32 { UINT64_C(44154446962804779), -19, true },
105-
::boost::decimal::decimal32 { UINT64_C(27400572833763747), -18 },
106-
::boost::decimal::decimal32 { UINT64_C(12289830364128736), -17, true },
107-
::boost::decimal::decimal32 { UINT64_C(40157034119189716), -17 },
108-
::boost::decimal::decimal32 { UINT64_C(94842703533437844), -17, true },
109-
::boost::decimal::decimal32 { UINT64_C(15738722141839421), -16 },
110-
::boost::decimal::decimal32 { UINT64_C(1425850153011925), -16, true },
111-
}};
112-
113-
static constexpr std::array<::boost::decimal::decimal64, 11> d64_coeffs_med =
114-
{{
115-
::boost::decimal::decimal64 { UINT64_C(35895331641408534), -24, true },
116-
::boost::decimal::decimal64 { UINT64_C(1734850544519432), -21 },
117-
::boost::decimal::decimal64 { UINT64_C(38064221484608425), -21, true },
118-
::boost::decimal::decimal64 { UINT64_C(50135011697517902), -20 },
119-
::boost::decimal::decimal64 { UINT64_C(44154446962804779), -19, true },
120-
::boost::decimal::decimal64 { UINT64_C(27400572833763747), -18 },
121-
::boost::decimal::decimal64 { UINT64_C(12289830364128736), -17, true },
122-
::boost::decimal::decimal64 { UINT64_C(40157034119189716), -17 },
123-
::boost::decimal::decimal64 { UINT64_C(94842703533437844), -17, true },
124-
::boost::decimal::decimal64 { UINT64_C(15738722141839421), -16 },
125-
::boost::decimal::decimal64 { UINT64_C(1425850153011925), -16, true },
126-
}};
127-
128-
static constexpr std::array<::boost::decimal::decimal128, 11> d128_coeffs_med =
129-
{{
130-
::boost::decimal::decimal128 { UINT64_C(35895331641408534), -24, true },
131-
::boost::decimal::decimal128 { UINT64_C(1734850544519432), -21 },
132-
::boost::decimal::decimal128 { UINT64_C(38064221484608425), -21, true },
133-
::boost::decimal::decimal128 { UINT64_C(50135011697517902), -20 },
134-
::boost::decimal::decimal128 { UINT64_C(44154446962804779), -19, true },
135-
::boost::decimal::decimal128 { UINT64_C(27400572833763747), -18 },
136-
::boost::decimal::decimal128 { UINT64_C(12289830364128736), -17, true },
137-
::boost::decimal::decimal128 { UINT64_C(40157034119189716), -17 },
138-
::boost::decimal::decimal128 { UINT64_C(94842703533437844), -17, true },
139-
::boost::decimal::decimal128 { UINT64_C(15738722141839421), -16 },
140-
::boost::decimal::decimal128 { UINT64_C(1425850153011925), -16, true },
141-
}};
142-
14380
};
14481

14582
#if !(defined(__cpp_inline_variables) && __cpp_inline_variables >= 201606L) && (!defined(_MSC_VER) || _MSC_VER != 1900)
14683

147-
template <bool b> constexpr std::array<decimal32, 11> atan_table_imp<b>::d32_coeffs_small;
148-
template <bool b> constexpr std::array<decimal32, 11> atan_table_imp<b>::d32_coeffs_med;
149-
template <bool b> constexpr std::array<decimal32, 3> atan_table_imp<b>::d32_atan_values;
84+
template <bool b> constexpr std::array<decimal32, 11> atan_table_imp<b>::d32_coeffs;
85+
template <bool b> constexpr std::array<decimal64, 11> atan_table_imp<b>::d64_coeffs;
15086

151-
template <bool b> constexpr std::array<decimal64, 11> atan_table_imp<b>::d64_coeffs_small;
152-
template <bool b> constexpr std::array<decimal64, 11> atan_table_imp<b>::d64_coeffs_med;
87+
template <bool b> constexpr std::array<decimal32, 3> atan_table_imp<b>::d32_atan_values;
15388
template <bool b> constexpr std::array<decimal64, 3> atan_table_imp<b>::d64_atan_values;
154-
155-
template <bool b> constexpr std::array<decimal128, 11> atan_table_imp<b>::d128_coeffs_small;
156-
template <bool b> constexpr std::array<decimal128, 11> atan_table_imp<b>::d128_coeffs_med;
15789
template <bool b> constexpr std::array<decimal128, 3> atan_table_imp<b>::d128_atan_values;
15890

15991
#endif
@@ -163,21 +95,56 @@ using atan_table = atan_table_imp<true>;
16395
} //namespace atan_detail
16496

16597
template <BOOST_DECIMAL_DECIMAL_FLOATING_TYPE T>
166-
constexpr auto atan_series_small(T x) noexcept;
167-
168-
template <BOOST_DECIMAL_DECIMAL_FLOATING_TYPE T>
169-
constexpr auto atan_series_med(T x) noexcept;
98+
constexpr auto atan_series(T x) noexcept;
17099

171100
template <BOOST_DECIMAL_DECIMAL_FLOATING_TYPE T>
172101
constexpr auto atan_values(std::size_t idx) noexcept -> T;
173102

174-
template <> constexpr auto atan_series_small<decimal32> (decimal32 x) noexcept { return remez_series_result(x, atan_detail::atan_table::d32_coeffs_small); }
175-
template <> constexpr auto atan_series_small<decimal64> (decimal64 x) noexcept { return remez_series_result(x, atan_detail::atan_table::d64_coeffs_small); }
176-
template <> constexpr auto atan_series_small<decimal128>(decimal128 x) noexcept { return remez_series_result(x, atan_detail::atan_table::d128_coeffs_small); }
103+
template <> constexpr auto atan_series<decimal32> (decimal32 x) noexcept { return remez_series_result(x, atan_detail::atan_table::d32_coeffs); }
104+
template <> constexpr auto atan_series<decimal64> (decimal64 x) noexcept { return remez_series_result(x, atan_detail::atan_table::d64_coeffs); }
177105

178-
template <> constexpr auto atan_series_med <decimal32> (decimal32 x) noexcept { return remez_series_result(x, atan_detail::atan_table::d32_coeffs_med ); }
179-
template <> constexpr auto atan_series_med <decimal64> (decimal64 x) noexcept { return remez_series_result(x, atan_detail::atan_table::d64_coeffs_med ); }
180-
template <> constexpr auto atan_series_med <decimal128>(decimal128 x) noexcept { return remez_series_result(x, atan_detail::atan_table::d128_coeffs_med ); }
106+
template <>
107+
constexpr auto atan_series<decimal128>(decimal128 x) noexcept
108+
{
109+
// PadeApproximant[ArcTan[x]/x, {x, 0, {18, 18}}]
110+
// FullSimplify[%]
111+
// HornerForm[Numerator[Out[2]]]
112+
// HornerForm[Denominator[Out[2]]]
113+
114+
const decimal128 x2 { x * x };
115+
116+
const decimal128
117+
top
118+
{
119+
decimal128 { UINT64_C(21427381364263875) }
120+
+ x2 * (decimal128 { UINT64_C(91886788553059500) }
121+
+ x2 * (decimal128 { UINT64_C(163675410390191700) }
122+
+ x2 * (decimal128 { UINT64_C(156671838074852100) }
123+
+ x2 * (decimal128 { UINT64_C(87054123957610810) }
124+
+ x2 * (decimal128 { UINT64_C(28283323008669300) }
125+
+ x2 * (decimal128 { UINT64_C(5134145876036100) }
126+
+ x2 * (decimal128 { UINT64_C(463911017673180) }
127+
+ x2 * (decimal128 { UINT64_C(16016872057515) }
128+
+ decimal128 { UINT64_C(90194313216) } * x2))))))))
129+
};
130+
131+
const decimal128
132+
bot
133+
{
134+
decimal128 { UINT64_C(21427381364263875) }
135+
+ x2 * (decimal128 { UINT64_C(99029249007814125) }
136+
+ x2 * (decimal128 { UINT64_C(192399683786610300) }
137+
+ x2 * (decimal128 { UINT64_C(204060270682768500) }
138+
+ x2 * (decimal128 { UINT64_C(128360492848838250) }
139+
+ x2 * (decimal128 { UINT64_C(48688462804731750) }
140+
+ x2 * (decimal128 { UINT64_C(10819658401051500) }
141+
+ x2 * (decimal128 { UINT64_C(1298359008126180) }
142+
+ x2 * (decimal128 { UINT64_C(70562989572075) }
143+
+ decimal128 { UINT64_C(1120047453525) } * x2))))))))
144+
};
145+
146+
return (x * top) / bot;
147+
}
181148

182149
template <> constexpr auto atan_values<decimal32> (std::size_t idx) noexcept -> decimal32 { return atan_detail::atan_table::d32_atan_values [idx]; }
183150
template <> constexpr auto atan_values<decimal64> (std::size_t idx) noexcept -> decimal64 { return atan_detail::atan_table::d64_atan_values [idx]; }

‎test/test_atan.cpp‎

Lines changed: 130 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -1,4 +1,5 @@
1-
// Copyright 2023 Matt Borland
1+
// Copyright 2023 - 2024 Matt Borland
2+
// Copyright 2023 - 2024 Christopher Kormanyos
23
// Distributed under the Boost Software License, Version 1.0.
34
// https://www.boost.org/LICENSE_1_0.txt
45

@@ -35,6 +36,47 @@ static constexpr auto N = static_cast<std::size_t>(128U); // Number of trials
3536
static constexpr auto N = static_cast<std::size_t>(128U >> 4U); // Number of trials
3637
#endif
3738

39+
namespace local {
40+
41+
template<typename NumericType>
42+
auto is_close_fraction(const NumericType& a,
43+
const NumericType& b,
44+
const NumericType& tol) noexcept -> bool
45+
{
46+
using std::fabs;
47+
48+
auto result_is_ok = bool { };
49+
50+
NumericType delta { };
51+
52+
if (b == static_cast<NumericType>(0))
53+
{
54+
delta = fabs(a - b); // LCOV_EXCL_LINE
55+
56+
result_is_ok = (delta < tol); // LCOV_EXCL_LINE
57+
}
58+
else
59+
{
60+
delta = fabs(1 - (a / b));
61+
62+
result_is_ok = (delta < tol);
63+
}
64+
65+
// LCOV_EXCL_START
66+
if (!result_is_ok)
67+
{
68+
std::cerr << std::setprecision(std::numeric_limits<NumericType>::digits10) << "a: " << a
69+
<< "\nb: " << b
70+
<< "\ndelta: " << delta
71+
<< "\ntol: " << tol << std::endl;
72+
}
73+
// LCOV_EXCL_STOP
74+
75+
return result_is_ok;
76+
}
77+
78+
} // namespace local
79+
3880
static std::mt19937_64 rng(42);
3981

4082
using namespace boost::decimal;
@@ -248,6 +290,87 @@ void test_atan()
248290
BOOST_TEST_EQ(atan(std::numeric_limits<Dec>::epsilon() * Dec(one(rng))), std::numeric_limits<Dec>::epsilon() * Dec(one(rng)));
249291
}
250292

293+
namespace local {
294+
295+
auto test_atan_128(const int tol_factor) -> bool
296+
{
297+
using decimal_type = boost::decimal::decimal128;
298+
299+
using str_ctrl_array_type = std::array<const char*, 31U>;
300+
301+
const str_ctrl_array_type ctrl_strings =
302+
{{
303+
// Table[N[ArcTan[n + (n + 1)/10], 36], {n, 0, 30, 1}]
304+
"0.0996686524911620273784461198780205902",
305+
"0.876058050598193423114047521128341339",
306+
"1.16066898625340562678011092078453218",
307+
"1.28474488507757839521660045035576124",
308+
"1.35212738092095465718914794138981285",
309+
"1.39408747072486000451142034998493574",
310+
"1.42263630606306524074609878711808965",
311+
"1.44328676857965836015625241442945425",
312+
"1.45890606062322050438578419322289952",
313+
"1.47112767430373459185287557176173085",
314+
"1.48094878710026889004729148794586818",
315+
"1.48901194617690061950062694387225355",
316+
"1.49574956319845608444166302563671922",
317+
"1.50146319310668803425604443731368245",
318+
"1.50636948736934306863178215633740183",
319+
"1.51062807563988690252010732385580363",
320+
"1.51435914848319329323768224134174846",
321+
"1.51765491794996116222569383300002088",
322+
"1.52058730451178540948599035411720525",
323+
"1.52321322351791322342928897562326592",
324+
"1.52557830188603652983022768173087618",
325+
"1.52771954287113490153496773234777134",
326+
"1.52966727041734455455350773831183303",
327+
"1.53144657040629129910455280928278240",
328+
"1.53307837432803240434373777942318448",
329+
"1.53458028472248503068348160776186418",
330+
"1.53596721141036644725997548254677890",
331+
"1.53725186723321821298607774053472251",
332+
"1.53844515819879435219259100906944202",
333+
"1.53955649336462834297760994674726047",
334+
"1.54059403307910435064686494555939664",
335+
}};
336+
337+
std::array<decimal_type, std::tuple_size<str_ctrl_array_type>::value> atan_values { };
338+
std::array<decimal_type, std::tuple_size<str_ctrl_array_type>::value> ctrl_values { };
339+
340+
int nx { 0 };
341+
342+
bool result_is_ok { true };
343+
344+
const decimal_type my_tol { std::numeric_limits<decimal_type>::epsilon() * static_cast<decimal_type>(tol_factor) };
345+
346+
for(auto i = static_cast<std::size_t>(UINT8_C(0)); i < std::tuple_size<str_ctrl_array_type>::value; ++i)
347+
{
348+
const decimal_type
349+
x_arg
350+
{
351+
decimal_type { nx }
352+
+ decimal_type { nx + 1, -1 }
353+
};
354+
355+
++nx;
356+
357+
atan_values[i] = atan(x_arg);
358+
359+
static_cast<void>
360+
(
361+
from_chars(ctrl_strings[i], ctrl_strings[i] + std::strlen(ctrl_strings[i]), ctrl_values[i])
362+
);
363+
364+
const auto result_atan_is_ok = is_close_fraction(atan_values[i], ctrl_values[i], my_tol);
365+
366+
result_is_ok = (result_atan_is_ok && result_is_ok);
367+
}
368+
369+
return result_is_ok;
370+
}
371+
372+
} // namespace local
373+
251374
int main()
252375
{
253376
test_atan<decimal32>();
@@ -256,5 +379,11 @@ int main()
256379
spot_test(0.344559F);
257380
spot_test(0.181179F);
258381

382+
{
383+
const auto result_pos128_is_ok = local::test_atan_128(800'000);
384+
385+
BOOST_TEST(result_pos128_is_ok);
386+
}
387+
259388
return boost::report_errors();
260389
}

0 commit comments

Comments
 (0)