math/special_functions/cbrt.hpp

0001 //  (C) Copyright John Maddock 2006.
0002 //  Use, modification and distribution are subject to the
0003 //  Boost Software License, Version 1.0. (See accompanying file
0004 //  LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)
0005
0006 #ifndef BOOST_MATH_SF_CBRT_HPP
0007 #define BOOST_MATH_SF_CBRT_HPP
0008
0009 #ifdef _MSC_VER
0010 #pragma once
0011 #endif
0012
0013 #include <boost/math/tools/rational.hpp>
0014 #include <boost/math/policies/error_handling.hpp>
0015 #include <boost/math/special_functions/math_fwd.hpp>
0016 #include <boost/math/special_functions/fpclassify.hpp>
0017 #include <type_traits>
0018 #include <cstdint>
0019
0020 namespace boost{ namespace math{
0021
0022 namespace detail
0023 {
0024
0025 struct big_int_type
0026 {
0027    operator std::uintmax_t() const;
0028 };
0029
0030 template <typename T>
0031 struct largest_cbrt_int_type
0032 {
0033    using type = typename std::conditional<
0034       std::is_convertible<big_int_type, T>::value,
0035       std::uintmax_t,
0036       unsigned int
0037    >::type;
0038 };
0039
0040 template <typename T, typename Policy>
0041 T cbrt_imp(T z, const Policy& pol)
0042 {
0043    BOOST_MATH_STD_USING
0044    //
0045    // cbrt approximation for z in the range [0.5,1]
0046    // It's hard to say what number of terms gives the optimum
0047    // trade off between precision and performance, this seems
0048    // to be about the best for double precision.
0049    //
0050    // Maximum Deviation Found:                     1.231e-006
0051    // Expected Error Term:                         -1.231e-006
0052    // Maximum Relative Change in Control Points:   5.982e-004
0053    //
0054    static const T P[] = {
0055       static_cast<T>(0.37568269008611818),
0056       static_cast<T>(1.3304968705558024),
0057       static_cast<T>(-1.4897101632445036),
0058       static_cast<T>(1.2875573098219835),
0059       static_cast<T>(-0.6398703759826468),
0060       static_cast<T>(0.13584489959258635),
0061    };
0062    static const T correction[] = {
0063       static_cast<T>(0.62996052494743658238360530363911),  // 2^-2/3
0064       static_cast<T>(0.79370052598409973737585281963615),  // 2^-1/3
0065       static_cast<T>(1),
0066       static_cast<T>(1.2599210498948731647672106072782),   // 2^1/3
0067       static_cast<T>(1.5874010519681994747517056392723),   // 2^2/3
0068    };
0069    if((boost::math::isinf)(z) || (z == 0))
0070       return z;
0071    if(!(boost::math::isfinite)(z))
0072    {
0073       return policies::raise_domain_error("boost::math::cbrt<%1%>(%1%)", "Argument to function must be finite but got %1%.", z, pol);
0074    }
0075
0076    int i_exp, sign(1);
0077    if(z < 0)
0078    {
0079       z = -z;
0080       sign = -sign;
0081    }
0082
0083    T guess = frexp(z, &i_exp);
0084    int original_i_exp = i_exp; // save for later
0085    guess = tools::evaluate_polynomial(P, guess);
0086    int i_exp3 = i_exp / 3;
0087
0088    using shift_type = typename largest_cbrt_int_type<T>::type;
0089
0090    static_assert( ::std::numeric_limits<shift_type>::radix == 2, "The radix of the type to shift to must be 2.");
0091
0092    if(abs(i_exp3) < std::numeric_limits<shift_type>::digits)
0093    {
0094       if(i_exp3 > 0)
0095          guess *= shift_type(1u) << i_exp3;
0096       else
0097          guess /= shift_type(1u) << -i_exp3;
0098    }
0099    else
0100    {
0101       guess = ldexp(guess, i_exp3);
0102    }
0103    i_exp %= 3;
0104    guess *= correction[i_exp + 2];
0105    //
0106    // Now inline Halley iteration.
0107    // We do this here rather than calling tools::halley_iterate since we can
0108    // simplify the expressions algebraically, and don't need most of the error
0109    // checking of the boilerplate version as we know in advance that the function
0110    // is well behaved...
0111    //
0112    using prec = typename policies::precision<T, Policy>::type;
0113    constexpr auto prec3 = prec::value / 3;
0114    constexpr auto new_prec = prec3 + 3;
0115    using new_policy = typename policies::normalise<Policy, policies::digits2<new_prec>>::type;
0116    //
0117    // Epsilon calculation uses compile time arithmetic when it's available for type T,
0118    // otherwise uses ldexp to calculate at runtime:
0119    //
0120    T eps = (new_prec > 3) ? policies::get_epsilon<T, new_policy>() : ldexp(T(1), -2 - tools::digits<T>() / 3);
0121    T diff;
0122
0123    if(original_i_exp < std::numeric_limits<T>::max_exponent - 3)
0124    {
0125       //
0126       // Safe from overflow, use the fast method:
0127       //
0128       do
0129       {
0130          T g3 = guess * guess * guess;
0131          diff = (g3 + z + z) / (g3 + g3 + z);
0132          guess *= diff;
0133       }
0134       while(fabs(1 - diff) > eps);
0135    }
0136    else
0137    {
0138       //
0139       // Either we're ready to overflow, or we can't tell because numeric_limits isn't
0140       // available for type T:
0141       //
0142       do
0143       {
0144          T g2 = guess * guess;
0145          diff = (g2 - z / guess) / (2 * guess + z / g2);
0146          guess -= diff;
0147       }
0148       while((guess * eps) < fabs(diff));
0149    }
0150
0151    return sign * guess;
0152 }
0153
0154 } // namespace detail
0155
0156 template <typename T, typename Policy>
0157 inline typename tools::promote_args<T>::type cbrt(T z, const Policy& pol)
0158 {
0159    using result_type = typename tools::promote_args<T>::type;
0160    using value_type = typename policies::evaluation<result_type, Policy>::type;
0161    return static_cast<result_type>(detail::cbrt_imp(value_type(z), pol));
0162 }
0163
0164 template <typename T>
0165 inline typename tools::promote_args<T>::type cbrt(T z)
0166 {
0167    return cbrt(z, policies::policy<>());
0168 }
0169
0170 } // namespace math
0171 } // namespace boost
0172
0173 #endif // BOOST_MATH_SF_CBRT_HPP
0174
0175
0176
0177