math/special_functions/cbrt.hpp

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