1// (C) Copyright John Maddock 2006.
2// Use, modification and distribution are subject to the
3// Boost Software License, Version 1.0. (See accompanying file
4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
5
6#ifndef BOOST_MATH_SF_CBRT_HPP
7#define BOOST_MATH_SF_CBRT_HPP
8
9#ifdef _MSC_VER
10#pragma once
11#endif
12
13#include <boost/math/tools/rational.hpp>
14#include <boost/math/policies/error_handling.hpp>
15#include <boost/math/special_functions/math_fwd.hpp>
16#include <boost/math/special_functions/fpclassify.hpp>
17#include <type_traits>
18#include <cstdint>
19
20namespace boost{ namespace math{
21
22namespace detail
23{
24
25struct big_int_type
26{
27 operator std::uintmax_t() const;
28};
29
30template <typename T>
31struct largest_cbrt_int_type
32{
33 using type = typename std::conditional<
34 std::is_convertible<big_int_type, T>::value,
35 std::uintmax_t,
36 unsigned int
37 >::type;
38};
39
40template <typename T, typename Policy>
41T cbrt_imp(T z, const Policy& pol)
42{
43 BOOST_MATH_STD_USING
44 //
45 // cbrt approximation for z in the range [0.5,1]
46 // It's hard to say what number of terms gives the optimum
47 // trade off between precision and performance, this seems
48 // to be about the best for double precision.
49 //
50 // Maximum Deviation Found: 1.231e-006
51 // Expected Error Term: -1.231e-006
52 // Maximum Relative Change in Control Points: 5.982e-004
53 //
54 static const T P[] = {
55 static_cast<T>(0.37568269008611818),
56 static_cast<T>(1.3304968705558024),
57 static_cast<T>(-1.4897101632445036),
58 static_cast<T>(1.2875573098219835),
59 static_cast<T>(-0.6398703759826468),
60 static_cast<T>(0.13584489959258635),
61 };
62 static const T correction[] = {
63 static_cast<T>(0.62996052494743658238360530363911), // 2^-2/3
64 static_cast<T>(0.79370052598409973737585281963615), // 2^-1/3
65 static_cast<T>(1),
66 static_cast<T>(1.2599210498948731647672106072782), // 2^1/3
67 static_cast<T>(1.5874010519681994747517056392723), // 2^2/3
68 };
69 if((boost::math::isinf)(z) || (z == 0))
70 return z;
71 if(!(boost::math::isfinite)(z))
72 {
73 return policies::raise_domain_error("boost::math::cbrt<%1%>(%1%)", "Argument to function must be finite but got %1%.", z, pol);
74 }
75
76 int i_exp, sign(1);
77 if(z < 0)
78 {
79 z = -z;
80 sign = -sign;
81 }
82
83 T guess = frexp(z, &i_exp);
84 int original_i_exp = i_exp; // save for later
85 guess = tools::evaluate_polynomial(P, guess);
86 int i_exp3 = i_exp / 3;
87
88 using shift_type = typename largest_cbrt_int_type<T>::type;
89
90 static_assert( ::std::numeric_limits<shift_type>::radix == 2, "The radix of the type to shift to must be 2.");
91
92 if(abs(i_exp3) < std::numeric_limits<shift_type>::digits)
93 {
94 if(i_exp3 > 0)
95 guess *= shift_type(1u) << i_exp3;
96 else
97 guess /= shift_type(1u) << -i_exp3;
98 }
99 else
100 {
101 guess = ldexp(guess, i_exp3);
102 }
103 i_exp %= 3;
104 guess *= correction[i_exp + 2];
105 //
106 // Now inline Halley iteration.
107 // We do this here rather than calling tools::halley_iterate since we can
108 // simplify the expressions algebraically, and don't need most of the error
109 // checking of the boilerplate version as we know in advance that the function
110 // is well behaved...
111 //
112 using prec = typename policies::precision<T, Policy>::type;
113 constexpr auto prec3 = prec::value / 3;
114 constexpr auto new_prec = prec3 + 3;
115 using new_policy = typename policies::normalise<Policy, policies::digits2<new_prec>>::type;
116 //
117 // Epsilon calculation uses compile time arithmetic when it's available for type T,
118 // otherwise uses ldexp to calculate at runtime:
119 //
120 T eps = (new_prec > 3) ? policies::get_epsilon<T, new_policy>() : ldexp(T(1), -2 - tools::digits<T>() / 3);
121 T diff;
122
123 if(original_i_exp < std::numeric_limits<T>::max_exponent - 3)
124 {
125 //
126 // Safe from overflow, use the fast method:
127 //
128 do
129 {
130 T g3 = guess * guess * guess;
131 diff = (g3 + z + z) / (g3 + g3 + z);
132 guess *= diff;
133 }
134 while(fabs(1 - diff) > eps);
135 }
136 else
137 {
138 //
139 // Either we're ready to overflow, or we can't tell because numeric_limits isn't
140 // available for type T:
141 //
142 do
143 {
144 T g2 = guess * guess;
145 diff = (g2 - z / guess) / (2 * guess + z / g2);
146 guess -= diff;
147 }
148 while((guess * eps) < fabs(diff));
149 }
150
151 return sign * guess;
152}
153
154} // namespace detail
155
156template <typename T, typename Policy>
157inline typename tools::promote_args<T>::type cbrt(T z, const Policy& pol)
158{
159 using result_type = typename tools::promote_args<T>::type;
160 using value_type = typename policies::evaluation<result_type, Policy>::type;
161 return static_cast<result_type>(detail::cbrt_imp(value_type(z), pol));
162}
163
164template <typename T>
165inline typename tools::promote_args<T>::type cbrt(T z)
166{
167 return cbrt(z, policies::policy<>());
168}
169
170} // namespace math
171} // namespace boost
172
173#endif // BOOST_MATH_SF_CBRT_HPP
174
175
176
177
178

source code of boost/libs/math/include/boost/math/special_functions/cbrt.hpp