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 <boost/mpl/divides.hpp>
18#include <boost/mpl/plus.hpp>
19#include <boost/mpl/if.hpp>
20#include <boost/type_traits/is_convertible.hpp>
21
22namespace boost{ namespace math{
23
24namespace detail
25{
26
27struct big_int_type
28{
29 operator boost::uintmax_t()const;
30};
31
32template <class T>
33struct largest_cbrt_int_type
34{
35 typedef typename mpl::if_c<
36 boost::is_convertible<big_int_type, T>::value,
37 boost::uintmax_t,
38 unsigned int
39 >::type type;
40};
41
42template <class T, class Policy>
43T cbrt_imp(T z, const Policy& pol)
44{
45 BOOST_MATH_STD_USING
46 //
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.
51 //
52 // Maximum Deviation Found: 1.231e-006
53 // Expected Error Term: -1.231e-006
54 // Maximum Relative Change in Control Points: 5.982e-004
55 //
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),
63 };
64 static const T correction[] = {
65 static_cast<T>(0.62996052494743658238360530363911), // 2^-2/3
66 static_cast<T>(0.79370052598409973737585281963615), // 2^-1/3
67 static_cast<T>(1),
68 static_cast<T>(1.2599210498948731647672106072782), // 2^1/3
69 static_cast<T>(1.5874010519681994747517056392723), // 2^2/3
70 };
71 if((boost::math::isinf)(z) || (z == 0))
72 return z;
73 if(!(boost::math::isfinite)(z))
74 {
75 return policies::raise_domain_error("boost::math::cbrt<%1%>(%1%)", "Argument to function must be finite but got %1%.", z, pol);
76 }
77
78 int i_exp, sign(1);
79 if(z < 0)
80 {
81 z = -z;
82 sign = -sign;
83 }
84
85 T guess = frexp(z, &i_exp);
86 int original_i_exp = i_exp; // save for later
87 guess = tools::evaluate_polynomial(P, guess);
88 int i_exp3 = i_exp / 3;
89
90 typedef typename largest_cbrt_int_type<T>::type shift_type;
91
92 BOOST_STATIC_ASSERT( ::std::numeric_limits<shift_type>::radix == 2);
93
94 if(abs(i_exp3) < std::numeric_limits<shift_type>::digits)
95 {
96 if(i_exp3 > 0)
97 guess *= shift_type(1u) << i_exp3;
98 else
99 guess /= shift_type(1u) << -i_exp3;
100 }
101 else
102 {
103 guess = ldexp(guess, i_exp3);
104 }
105 i_exp %= 3;
106 guess *= correction[i_exp + 2];
107 //
108 // Now inline Halley iteration.
109 // We do this here rather than calling tools::halley_iterate since we can
110 // simplify the expressions algebraically, and don't need most of the error
111 // checking of the boilerplate version as we know in advance that the function
112 // is well behaved...
113 //
114 typedef typename policies::precision<T, Policy>::type prec;
115 typedef typename mpl::divides<prec, boost::integral_constant<int, 3> >::type prec3;
116 typedef typename mpl::plus<prec3, boost::integral_constant<int, 3> >::type new_prec;
117 typedef typename policies::normalise<Policy, policies::digits2<new_prec::value> >::type new_policy;
118 //
119 // Epsilon calculation uses compile time arithmetic when it's available for type T,
120 // otherwise uses ldexp to calculate at runtime:
121 //
122 T eps = (new_prec::value > 3) ? policies::get_epsilon<T, new_policy>() : ldexp(T(1), -2 - tools::digits<T>() / 3);
123 T diff;
124
125 if(original_i_exp < std::numeric_limits<T>::max_exponent - 3)
126 {
127 //
128 // Safe from overflow, use the fast method:
129 //
130 do
131 {
132 T g3 = guess * guess * guess;
133 diff = (g3 + z + z) / (g3 + g3 + z);
134 guess *= diff;
135 }
136 while(fabs(1 - diff) > eps);
137 }
138 else
139 {
140 //
141 // Either we're ready to overflow, or we can't tell because numeric_limits isn't
142 // available for type T:
143 //
144 do
145 {
146 T g2 = guess * guess;
147 diff = (g2 - z / guess) / (2 * guess + z / g2);
148 guess -= diff;
149 }
150 while((guess * eps) < fabs(diff));
151 }
152
153 return sign * guess;
154}
155
156} // namespace detail
157
158template <class T, class Policy>
159inline typename tools::promote_args<T>::type cbrt(T z, const Policy& pol)
160{
161 typedef typename tools::promote_args<T>::type result_type;
162 typedef typename policies::evaluation<result_type, Policy>::type value_type;
163 return static_cast<result_type>(detail::cbrt_imp(value_type(z), pol));
164}
165
166template <class T>
167inline typename tools::promote_args<T>::type cbrt(T z)
168{
169 return cbrt(z, policies::policy<>());
170}
171
172} // namespace math
173} // namespace boost
174
175#endif // BOOST_MATH_SF_CBRT_HPP
176
177
178
179
180

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