Skip to content

Commit e72fd34

Browse files
ckormanyosckormanyos
authored
Merge pull request #1328 from boostorg/rework_sqrt_cbrt
Rework sqrt cbrt
2 parents be078b6 + 4f156cd commit e72fd34

4 files changed

Lines changed: 123 additions & 82 deletions

File tree

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

Lines changed: 28 additions & 18 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
1-
// Copyright 2023 - 2024 Matt Borland
2-
// Copyright 2023 - 2024 Christopher Kormanyos
1+
// Copyright 2023 - 2026 Matt Borland
2+
// Copyright 2023 - 2026 Christopher Kormanyos
33
// Distributed under the Boost Software License, Version 1.0.
44
// https://www.boost.org/LICENSE_1_0.txt
55

@@ -103,42 +103,47 @@ constexpr auto cbrt_impl(const T x) noexcept
103103
}
104104
else
105105
{
106-
// Scale the argument to the interval 1/10 <= x < 1.
106+
// Scale the argument to the interval 1/64 < x <= 1/8.
107107
T gx { gn, -std::numeric_limits<T>::digits10 };
108108

109+
const bool is_scaled_by_eight { (gx > T { 125, -3 }) };
110+
111+
if(is_scaled_by_eight)
112+
{
113+
gx /= 8;
114+
}
115+
109116
exp10val += std::numeric_limits<T>::digits10;
110117

111118
// For this work we perform an order-2 Pade approximation of the cube-root
112-
// at argument x = 1/2. This results in slightly more than 2 decimal digits
113-
// of accuracy over the interval 1/10 <= x < 1.
119+
// at argument x = 1/16. This results in slightly more than 4 decimal digits
120+
// of accuracy over the interval 1/64 < x <= 1/8.
114121

115-
// PadeApproximant[x^(1/3), {x, 1/2, {2, 2}}]
122+
// PadeApproximant[x^(1/3), {x, 1/16, {2, 2}}]
116123
// FullSimplify[%]
117124

118125
// HornerForm[Numerator[Out[2]]]
119126
// Results in:
120-
// 5 + x (70 + 56 x)
127+
// 5 + x (560 + 3584 x)
121128

122129
// HornerForm[Denominator[Out[2]]]
123130
// Results in:
124-
// 2^(1/3) (14 + x (70 + 20 x))
131+
// 2^(1/3) (28 + x (1120 + 2560 x))
125132

126-
constexpr T five { 5 };
127-
constexpr T fourteen { 14 };
128-
constexpr T seventy { 7, 1 };
133+
constexpr T five { 5 };
129134

130135
result =
131-
(five + gx * (seventy + gx * 56))
132-
/ (numbers::cbrt2_v<T> * (fourteen + gx * (seventy + gx * 20)));
136+
(five + gx * (560 + 3584 * gx))
137+
/ (numbers::cbrt2_v<T> * (28 + gx * (1120 + 2560 * gx)));
133138

134-
// Perform 3, 4 or 5 Newton-Raphson iterations depending on precision.
135-
// Note from above, we start with slightly more than 2 decimal digits
139+
// Perform 2, 3 or 4 Newton-Raphson iterations depending on precision.
140+
// Note from above, we start with slightly more than 4 decimal digits
136141
// of accuracy.
137142

138143
constexpr int iter_loops
139144
{
140-
std::numeric_limits<T>::digits10 < 10 ? 3
141-
: std::numeric_limits<T>::digits10 < 20 ? 4 : 5
145+
std::numeric_limits<T>::digits10 < 10 ? 2
146+
: std::numeric_limits<T>::digits10 < 20 ? 3 : 4
142147
};
143148

144149
for (int idx = 0; idx < iter_loops; ++idx)
@@ -172,6 +177,11 @@ constexpr auto cbrt_impl(const T x) noexcept
172177
break;
173178
}
174179
}
180+
181+
if(is_scaled_by_eight)
182+
{
183+
result *= 2;
184+
}
175185
}
176186
}
177187

@@ -181,7 +191,7 @@ constexpr auto cbrt_impl(const T x) noexcept
181191
} //namespace detail
182192

183193
BOOST_DECIMAL_EXPORT template <typename T>
184-
constexpr auto cbrt(const T val) noexcept
194+
constexpr auto cbrt(const T val) noexcept // LCOV_EXCL_LINE
185195
BOOST_DECIMAL_REQUIRES(detail::is_decimal_floating_point_v, T)
186196
{
187197
using evaluation_type = detail::evaluation_type_t<T>;

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

Lines changed: 29 additions & 15 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
1-
// Copyright 2023 - 2024 Matt Borland
2-
// Copyright 2023 - 2024 Christopher Kormanyos
1+
// Copyright 2023 - 2026 Matt Borland
2+
// Copyright 2023 - 2026 Christopher Kormanyos
33
// Distributed under the Boost Software License, Version 1.0.
44
// https://www.boost.org/LICENSE_1_0.txt
55

@@ -94,43 +94,52 @@ constexpr auto sqrt_impl(const T x) noexcept
9494
}
9595
else
9696
{
97-
// Scale the argument to the interval 1/10 <= x < 1.
97+
// Scale the argument to the interval 1/16 < x <= 1/4.
9898
T gx { gn, -std::numeric_limits<T>::digits10 };
9999

100+
const bool is_scaled_by_four { (gx > T { 25, -2 }) };
101+
102+
if(is_scaled_by_four)
103+
{
104+
gx /= 4;
105+
}
106+
100107
exp10val += std::numeric_limits<T>::digits10;
101108

102109
// For this work we perform an order-2 Pade approximation of the square root
103-
// at argument x = 1/2. This results in slightly more than 2 decimal digits
104-
// of accuracy over the interval 1/10 <= x < 1.
110+
// at argument x = 1/8. This results in slightly more than 4 decimal digits
111+
// of accuracy over the interval 1/16 < x <= 1/4.
105112

106113
// See also WolframAlpha at:
107114
// https://www.wolframalpha.com/input?i=FullSimplify%5BPadeApproximant%5BSqrt%5Bx%5D%2C+%7Bx%2C+1%2F2%2C+%7B2%2C+2%7D%7D%5D%5D
108115

109-
// PadeApproximant[Sqrt[x], {x, 1/2, {2, 2}}]
116+
// PadeApproximant[Sqrt[x], {x, 1/8, {2, 2}}]
110117
// FullSimplify[%]
111118

112119
// HornerForm[Numerator[Out[2]]]
113120
// Results in:
114-
// 1 + x (20 + 20 x)
121+
// 1 + x (80 + 320 x)
115122

116123
// HornerForm[Denominator[Out[2]]]
117124
// Results in:
118-
// 5 Sqrt[2] + x (20 Sqrt[2] + 4 Sqrt[2] x)
125+
// 10 Sqrt[2] + x (160 Sqrt[2] + 128 Sqrt[2] x)
126+
// which simplifies to
127+
// Sqrt[2] (10 + x (160 + 128 x))
119128

120-
constexpr T five { 5 };
129+
constexpr T ten { 1, 1 };
121130

122131
result =
123-
(one + gx * ((one + gx) * 20))
124-
/ (numbers::sqrt2_v<T> * ((gx * 4) * (five + gx) + five));
132+
(one + gx * (80 + 320 * gx))
133+
/ (numbers::sqrt2_v<T> * (ten + gx * (160 + 128 * gx)));
125134

126-
// Perform 3, 4 or 5 Newton-Raphson iterations depending on precision.
127-
// Note from above, we start with slightly more than 2 decimal digits
135+
// Perform 2, 3 or 4 Newton-Raphson iterations depending on precision.
136+
// Note from above, we start with slightly more than 4 decimal digits
128137
// of accuracy.
129138

130139
constexpr int iter_loops
131140
{
132-
std::numeric_limits<T>::digits10 < 10 ? 3
133-
: std::numeric_limits<T>::digits10 < 20 ? 4 : 5
141+
std::numeric_limits<T>::digits10 < 10 ? 2
142+
: std::numeric_limits<T>::digits10 < 20 ? 3 : 4
134143
};
135144

136145
for (int idx = 0; idx < iter_loops; ++idx)
@@ -155,6 +164,11 @@ constexpr auto sqrt_impl(const T x) noexcept
155164
result /= numbers::sqrt10_v<T>;
156165
}
157166
}
167+
168+
if(is_scaled_by_four)
169+
{
170+
result *= 2;
171+
}
158172
}
159173
}
160174

‎test/test_cbrt.cpp‎

Lines changed: 20 additions & 30 deletions
Original file line numberDiff line numberDiff line change
@@ -1,5 +1,5 @@
1-
// Copyright 2024 Matt Borland
2-
// Copyright 2024 Christopher Kormanyos
1+
// Copyright 2024 - 2026 Matt Borland
2+
// Copyright 2024 - 2026 Christopher Kormanyos
33
// Distributed under the Boost Software License, Version 1.0.
44
// https://www.boost.org/LICENSE_1_0.txt
55

@@ -50,22 +50,12 @@ namespace local
5050
{
5151
using std::fabs;
5252

53-
auto result_is_ok = bool { };
54-
55-
NumericType delta { };
56-
57-
if(b == static_cast<NumericType>(0))
58-
{
59-
delta = fabs(a - b); // LCOV_EXCL_LINE
60-
61-
result_is_ok = (delta < tol); // LCOV_EXCL_LINE
62-
}
63-
else
53+
const NumericType delta
6454
{
65-
delta = fabs(1 - (a / b));
55+
(b == static_cast<NumericType>(0)) ? fabs(a - b) : fabs(1 - (a / b))
56+
};
6657

67-
result_is_ok = (delta < tol);
68-
}
58+
const bool result_is_ok { (delta < tol) };
6959

7060
// LCOV_EXCL_START
7161
if (!result_is_ok)
@@ -110,9 +100,9 @@ namespace local
110100
auto trials = static_cast<std::uint32_t>(UINT8_C(0));
111101

112102
#if !defined(BOOST_DECIMAL_REDUCE_TEST_DEPTH)
113-
constexpr auto count = (sizeof(decimal_type) == static_cast<std::size_t>(UINT8_C(4))) ? UINT32_C(0x400) : UINT32_C(0x40);
103+
constexpr std::uint32_t count { ((std::numeric_limits<DecimalType>::digits10 < 10) ? UINT16_C(3200) : UINT16_C(1600)) };
114104
#else
115-
constexpr auto count = (sizeof(decimal_type) == static_cast<std::size_t>(UINT8_C(4))) ? UINT32_C(0x40) : UINT32_C(0x4);
105+
constexpr std::uint32_t count { ((std::numeric_limits<DecimalType>::digits10 < 10) ? UINT16_C(320) : UINT16_C(160)) };
116106
#endif
117107

118108
for( ; trials < count; ++trials)
@@ -341,9 +331,9 @@ int main()
341331
using decimal_type = boost::decimal::decimal32_t;
342332
using float_type = float;
343333

344-
const auto result_small_is_ok = local::test_cbrt<decimal_type, float_type>(INT32_C(16), 1.0E-26L, 1.0E-01L);
345-
const auto result_medium_is_ok = local::test_cbrt<decimal_type, float_type>(INT32_C(16), 0.9E-01L, 1.1E+01L);
346-
const auto result_large_is_ok = local::test_cbrt<decimal_type, float_type>(INT32_C(16), 1.0E+01L, 1.0E+26L);
334+
const auto result_small_is_ok = local::test_cbrt<decimal_type, float_type>(16, 1.0E-26L, 1.0E-01L);
335+
const auto result_medium_is_ok = local::test_cbrt<decimal_type, float_type>(16, 0.9E-02L, 1.1E+01L);
336+
const auto result_large_is_ok = local::test_cbrt<decimal_type, float_type>(16, 1.0E+01L, 1.0E+26L);
347337

348338
BOOST_TEST(result_small_is_ok);
349339
BOOST_TEST(result_medium_is_ok);
@@ -361,12 +351,12 @@ int main()
361351
}
362352

363353
{
364-
using decimal_type = boost::decimal::decimal_fast32_t;
365-
using float_type = float;
354+
using decimal_type = boost::decimal::decimal64_t;
355+
using float_type = double;
366356

367-
const auto result_small_is_ok = local::test_cbrt<decimal_type, float_type>(INT32_C(16), 1.0E-26L, 1.0E-01L);
368-
const auto result_medium_is_ok = local::test_cbrt<decimal_type, float_type>(INT32_C(16), 0.9E-01L, 1.1E+01L);
369-
const auto result_large_is_ok = local::test_cbrt<decimal_type, float_type>(INT32_C(16), 1.0E+01L, 1.0E+26L);
357+
const auto result_small_is_ok = local::test_cbrt<decimal_type, float_type>(16, 1.0E-26L, 1.0E-01L);
358+
const auto result_medium_is_ok = local::test_cbrt<decimal_type, float_type>(16, 0.9E-01L, 1.1E+01L);
359+
const auto result_large_is_ok = local::test_cbrt<decimal_type, float_type>(16, 1.0E+01L, 1.0E+26L);
370360

371361
BOOST_TEST(result_small_is_ok);
372362
BOOST_TEST(result_medium_is_ok);
@@ -384,7 +374,7 @@ int main()
384374
}
385375

386376
{
387-
using decimal_type = boost::decimal::decimal64_t;
377+
using decimal_type = boost::decimal::decimal_fast64_t;
388378
using float_type = double;
389379

390380
const auto result_small_is_ok = local::test_cbrt<decimal_type, float_type>(INT32_C(16), 1.0E-76L, 1.0E-01L);
@@ -407,14 +397,14 @@ int main()
407397
}
408398

409399
{
410-
const auto result_cbrt128_is_ok = local::test_cbrt_128(96);
400+
const auto result_cbrt128_is_ok = local::test_cbrt_128(16);
411401

412402
BOOST_TEST(result_cbrt128_is_ok);
413403

414404
result_is_ok = (result_cbrt128_is_ok && result_is_ok);
415405
}
416406

417-
result_is_ok = ((boost::report_errors() == 0) && result_is_ok);
407+
BOOST_TEST(result_is_ok);
418408

419-
return (result_is_ok ? 0 : -1);
409+
return boost::report_errors();
420410
}

0 commit comments

Comments
 (0)