// (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 #include #include #include #include #include #include #include namespace boost{ namespace math{ namespace detail { struct big_int_type { operator boost::uintmax_t()const; }; template struct largest_cbrt_int_type { typedef typename mpl::if_< boost::is_convertible, boost::uintmax_t, unsigned int >::type type; }; template 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(0.37568269008611818), static_cast(1.3304968705558024), static_cast(-1.4897101632445036), static_cast(1.2875573098219835), static_cast(-0.6398703759826468), static_cast(0.13584489959258635), }; static const T correction[] = { static_cast(0.62996052494743658238360530363911), // 2^-2/3 static_cast(0.79370052598409973737585281963615), // 2^-1/3 static_cast(1), static_cast(1.2599210498948731647672106072782), // 2^1/3 static_cast(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::type shift_type; BOOST_STATIC_ASSERT( ::std::numeric_limits::radix == 2); if(abs(i_exp3) < std::numeric_limits::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::type prec; typedef typename mpl::divides >::type prec3; typedef typename mpl::plus >::type new_prec; typedef typename policies::normalise >::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() : ldexp(T(1), -2 - tools::digits() / 3); T diff; if(original_i_exp < std::numeric_limits::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 inline typename tools::promote_args::type cbrt(T z, const Policy& pol) { typedef typename tools::promote_args::type result_type; typedef typename policies::evaluation::type value_type; return static_cast(detail::cbrt_imp(value_type(z), pol)); } template inline typename tools::promote_args::type cbrt(T z) { return cbrt(z, policies::policy<>()); } } // namespace math } // namespace boost #endif // BOOST_MATH_SF_CBRT_HPP