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)
6
#ifndef BOOST_MATH_SF_CBRT_HPP
7
#define BOOST_MATH_SF_CBRT_HPP
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 <boost/mpl/divides.hpp>
18
#include <boost/mpl/plus.hpp>
19
#include <boost/mpl/if.hpp>
20
#include <boost/type_traits/is_convertible.hpp>
22
namespace boost{ namespace math{
29
operator boost::uintmax_t()const;
33
struct largest_cbrt_int_type
35
typedef typename mpl::if_<
36
boost::is_convertible<big_int_type, T>,
42
template <class T, class Policy>
43
T cbrt_imp(T z, const Policy& pol)
47
// cbrt approximation for z in the range [0.5,1]
48
// It's hard to say what number of terms gives the optimum
49
// trade off between precision and performance, this seems
50
// to be about the best for double precision.
52
// Maximum Deviation Found: 1.231e-006
53
// Expected Error Term: -1.231e-006
54
// Maximum Relative Change in Control Points: 5.982e-004
56
static const T P[] = {
57
static_cast<T>(0.37568269008611818),
58
static_cast<T>(1.3304968705558024),
59
static_cast<T>(-1.4897101632445036),
60
static_cast<T>(1.2875573098219835),
61
static_cast<T>(-0.6398703759826468),
62
static_cast<T>(0.13584489959258635),
64
static const T correction[] = {
65
static_cast<T>(0.62996052494743658238360530363911), // 2^-2/3
66
static_cast<T>(0.79370052598409973737585281963615), // 2^-1/3
68
static_cast<T>(1.2599210498948731647672106072782), // 2^1/3
69
static_cast<T>(1.5874010519681994747517056392723), // 2^2/3
72
if(!(boost::math::isfinite)(z))
74
return policies::raise_domain_error("boost::math::cbrt<%1%>(%1%)", "Argument to function must be finite but got %1%.", z, pol);
86
T guess = frexp(z, &i_exp);
87
int original_i_exp = i_exp; // save for later
88
guess = tools::evaluate_polynomial(P, guess);
89
int i_exp3 = i_exp / 3;
91
typedef typename largest_cbrt_int_type<T>::type shift_type;
93
BOOST_STATIC_ASSERT( ::std::numeric_limits<shift_type>::radix == 2);
95
if(abs(i_exp3) < std::numeric_limits<shift_type>::digits)
98
guess *= shift_type(1u) << i_exp3;
100
guess /= shift_type(1u) << -i_exp3;
104
guess = ldexp(guess, i_exp3);
107
guess *= correction[i_exp + 2];
109
// Now inline Halley iteration.
110
// We do this here rather than calling tools::halley_iterate since we can
111
// simplify the expressions algebraically, and don't need most of the error
112
// checking of the boilerplate version as we know in advance that the function
113
// is well behaved...
115
typedef typename policies::precision<T, Policy>::type prec;
116
typedef typename mpl::divides<prec, mpl::int_<3> >::type prec3;
117
typedef typename mpl::plus<prec3, mpl::int_<3> >::type new_prec;
118
typedef typename policies::normalise<Policy, policies::digits2<new_prec::value> >::type new_policy;
120
// Epsilon calculation uses compile time arithmetic when it's available for type T,
121
// otherwise uses ldexp to calculate at runtime:
123
T eps = (new_prec::value > 3) ? policies::get_epsilon<T, new_policy>() : ldexp(T(1), -2 - tools::digits<T>() / 3);
126
if(original_i_exp < std::numeric_limits<T>::max_exponent - 3)
129
// Safe from overflow, use the fast method:
133
T g3 = guess * guess * guess;
134
diff = (g3 + z + z) / (g3 + g3 + z);
137
while(fabs(1 - diff) > eps);
142
// Either we're ready to overflow, or we can't tell because numeric_limits isn't
143
// available for type T:
147
T g2 = guess * guess;
148
diff = (g2 - z / guess) / (2 * guess + z / g2);
151
while((guess * eps) < fabs(diff));
157
} // namespace detail
159
template <class T, class Policy>
160
inline typename tools::promote_args<T>::type cbrt(T z, const Policy& pol)
162
typedef typename tools::promote_args<T>::type result_type;
163
typedef typename policies::evaluation<result_type, Policy>::type value_type;
164
return static_cast<result_type>(detail::cbrt_imp(value_type(z), pol));
168
inline typename tools::promote_args<T>::type cbrt(T z)
170
return cbrt(z, policies::policy<>());
176
#endif // BOOST_MATH_SF_CBRT_HPP