]>
Commit | Line | Data |
---|---|---|
7c673cae FG |
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 | ||
22 | namespace boost{ namespace math{ | |
23 | ||
24 | namespace detail | |
25 | { | |
26 | ||
27 | struct big_int_type | |
28 | { | |
29 | operator boost::uintmax_t()const; | |
30 | }; | |
31 | ||
32 | template <class T> | |
33 | struct largest_cbrt_int_type | |
34 | { | |
f67539c2 TL |
35 | typedef typename mpl::if_c< |
36 | boost::is_convertible<big_int_type, T>::value, | |
7c673cae FG |
37 | boost::uintmax_t, |
38 | unsigned int | |
39 | >::type type; | |
40 | }; | |
41 | ||
42 | template <class T, class Policy> | |
43 | T 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 | }; | |
b32b8144 FG |
71 | if((boost::math::isinf)(z) || (z == 0)) |
72 | return z; | |
7c673cae FG |
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 | } | |
7c673cae FG |
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; | |
f67539c2 TL |
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; | |
7c673cae FG |
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 | ||
158 | template <class T, class Policy> | |
159 | inline 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 | ||
166 | template <class T> | |
167 | inline 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 |