| // (C) Copyright John Maddock 2006. |
| // Use, modification and distribution are subject to the |
| // Boost Software License, Version 1.0. (See accompanying file |
| // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) |
| |
| #ifndef BOOST_MATH_SF_CBRT_HPP |
| #define BOOST_MATH_SF_CBRT_HPP |
| |
| #ifdef _MSC_VER |
| #pragma once |
| #endif |
| |
| #include <boost/math/tools/rational.hpp> |
| #include <boost/math/policies/error_handling.hpp> |
| #include <boost/math/special_functions/math_fwd.hpp> |
| #include <boost/math/special_functions/fpclassify.hpp> |
| #include <boost/mpl/divides.hpp> |
| #include <boost/mpl/plus.hpp> |
| #include <boost/mpl/if.hpp> |
| #include <boost/type_traits/is_convertible.hpp> |
| |
| namespace boost{ namespace math{ |
| |
| namespace detail |
| { |
| |
| struct big_int_type |
| { |
| operator boost::uintmax_t()const; |
| }; |
| |
| template <class T> |
| struct largest_cbrt_int_type |
| { |
| typedef typename mpl::if_< |
| boost::is_convertible<big_int_type, T>, |
| boost::uintmax_t, |
| unsigned int |
| >::type type; |
| }; |
| |
| template <class T, class Policy> |
| T cbrt_imp(T z, const Policy& pol) |
| { |
| BOOST_MATH_STD_USING |
| // |
| // cbrt approximation for z in the range [0.5,1] |
| // It's hard to say what number of terms gives the optimum |
| // trade off between precision and performance, this seems |
| // to be about the best for double precision. |
| // |
| // Maximum Deviation Found: 1.231e-006 |
| // Expected Error Term: -1.231e-006 |
| // Maximum Relative Change in Control Points: 5.982e-004 |
| // |
| static const T P[] = { |
| static_cast<T>(0.37568269008611818), |
| static_cast<T>(1.3304968705558024), |
| static_cast<T>(-1.4897101632445036), |
| static_cast<T>(1.2875573098219835), |
| static_cast<T>(-0.6398703759826468), |
| static_cast<T>(0.13584489959258635), |
| }; |
| static const T correction[] = { |
| static_cast<T>(0.62996052494743658238360530363911), // 2^-2/3 |
| static_cast<T>(0.79370052598409973737585281963615), // 2^-1/3 |
| static_cast<T>(1), |
| static_cast<T>(1.2599210498948731647672106072782), // 2^1/3 |
| static_cast<T>(1.5874010519681994747517056392723), // 2^2/3 |
| }; |
| |
| if(!(boost::math::isfinite)(z)) |
| { |
| return policies::raise_domain_error("boost::math::cbrt<%1%>(%1%)", "Argument to function must be finite but got %1%.", z, pol); |
| } |
| |
| int i_exp, sign(1); |
| if(z < 0) |
| { |
| z = -z; |
| sign = -sign; |
| } |
| if(z == 0) |
| return 0; |
| |
| T guess = frexp(z, &i_exp); |
| int original_i_exp = i_exp; // save for later |
| guess = tools::evaluate_polynomial(P, guess); |
| int i_exp3 = i_exp / 3; |
| |
| typedef typename largest_cbrt_int_type<T>::type shift_type; |
| |
| if(abs(i_exp3) < std::numeric_limits<shift_type>::digits) |
| { |
| if(i_exp3 > 0) |
| guess *= shift_type(1u) << i_exp3; |
| else |
| guess /= shift_type(1u) << -i_exp3; |
| } |
| else |
| { |
| guess = ldexp(guess, i_exp3); |
| } |
| i_exp %= 3; |
| guess *= correction[i_exp + 2]; |
| // |
| // Now inline Halley iteration. |
| // We do this here rather than calling tools::halley_iterate since we can |
| // simplify the expressions algebraically, and don't need most of the error |
| // checking of the boilerplate version as we know in advance that the function |
| // is well behaved... |
| // |
| typedef typename policies::precision<T, Policy>::type prec; |
| typedef typename mpl::divides<prec, mpl::int_<3> >::type prec3; |
| typedef typename mpl::plus<prec3, mpl::int_<3> >::type new_prec; |
| typedef typename policies::normalise<Policy, policies::digits2<new_prec::value> >::type new_policy; |
| // |
| // Epsilon calculation uses compile time arithmetic when it's available for type T, |
| // otherwise uses ldexp to calculate at runtime: |
| // |
| T eps = (new_prec::value > 3) ? policies::get_epsilon<T, new_policy>() : ldexp(T(1), -2 - tools::digits<T>() / 3); |
| T diff; |
| |
| if(original_i_exp < std::numeric_limits<T>::max_exponent - 3) |
| { |
| // |
| // Safe from overflow, use the fast method: |
| // |
| do |
| { |
| T g3 = guess * guess * guess; |
| diff = (g3 + z + z) / (g3 + g3 + z); |
| guess *= diff; |
| } |
| while(fabs(1 - diff) > eps); |
| } |
| else |
| { |
| // |
| // Either we're ready to overflow, or we can't tell because numeric_limits isn't |
| // available for type T: |
| // |
| do |
| { |
| T g2 = guess * guess; |
| diff = (g2 - z / guess) / (2 * guess + z / g2); |
| guess -= diff; |
| } |
| while((guess * eps) < fabs(diff)); |
| } |
| |
| return sign * guess; |
| } |
| |
| } // namespace detail |
| |
| template <class T, class Policy> |
| inline typename tools::promote_args<T>::type cbrt(T z, const Policy& pol) |
| { |
| typedef typename tools::promote_args<T>::type result_type; |
| typedef typename policies::evaluation<result_type, Policy>::type value_type; |
| return static_cast<result_type>(detail::cbrt_imp(value_type(z), pol)); |
| } |
| |
| template <class T> |
| inline typename tools::promote_args<T>::type cbrt(T z) |
| { |
| return cbrt(z, policies::policy<>()); |
| } |
| |
| } // namespace math |
| } // namespace boost |
| |
| #endif // BOOST_MATH_SF_CBRT_HPP |
| |
| |
| |
| |