216 lines · plain
1// (C) Copyright John Maddock 2006.2// (C) Copyright Matt Borland 2024.3// Use, modification and distribution are subject to the4// Boost Software License, Version 1.0. (See accompanying file5// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)6 7#ifndef BOOST_MATH_SF_CBRT_HPP8#define BOOST_MATH_SF_CBRT_HPP9 10#ifdef _MSC_VER11#pragma once12#endif13 14#include <boost/math/tools/config.hpp>15 16#ifndef BOOST_MATH_HAS_NVRTC17 18#include <boost/math/tools/rational.hpp>19#include <boost/math/tools/type_traits.hpp>20#include <boost/math/tools/cstdint.hpp>21#include <boost/math/policies/error_handling.hpp>22#include <boost/math/special_functions/math_fwd.hpp>23#include <boost/math/special_functions/fpclassify.hpp>24 25namespace boost{ namespace math{26 27namespace detail28{29 30struct big_int_type31{32 operator std::uintmax_t() const;33};34 35template <typename T>36struct largest_cbrt_int_type37{38 using type = typename std::conditional<39 std::is_convertible<big_int_type, T>::value,40 std::uintmax_t,41 unsigned int42 >::type;43};44 45template <typename T, typename Policy>46BOOST_MATH_GPU_ENABLED T cbrt_imp(T z, const Policy& pol)47{48 BOOST_MATH_STD_USING49 //50 // cbrt approximation for z in the range [0.5,1]51 // It's hard to say what number of terms gives the optimum52 // trade off between precision and performance, this seems53 // to be about the best for double precision.54 //55 // Maximum Deviation Found: 1.231e-00656 // Expected Error Term: -1.231e-00657 // Maximum Relative Change in Control Points: 5.982e-00458 //59 BOOST_MATH_STATIC const T P[] = { 60 static_cast<T>(0.37568269008611818),61 static_cast<T>(1.3304968705558024),62 static_cast<T>(-1.4897101632445036),63 static_cast<T>(1.2875573098219835),64 static_cast<T>(-0.6398703759826468),65 static_cast<T>(0.13584489959258635),66 };67 BOOST_MATH_STATIC const T correction[] = {68 static_cast<T>(0.62996052494743658238360530363911), // 2^-2/369 static_cast<T>(0.79370052598409973737585281963615), // 2^-1/370 static_cast<T>(1),71 static_cast<T>(1.2599210498948731647672106072782), // 2^1/372 static_cast<T>(1.5874010519681994747517056392723), // 2^2/373 };74 if((boost::math::isinf)(z) || (z == 0))75 return z;76 if(!(boost::math::isfinite)(z))77 {78 return policies::raise_domain_error("boost::math::cbrt<%1%>(%1%)", "Argument to function must be finite but got %1%.", z, pol);79 }80 81 int i_exp, sign(1);82 if(z < 0)83 {84 z = -z;85 sign = -sign;86 }87 88 T guess = frexp(z, &i_exp);89 int original_i_exp = i_exp; // save for later90 guess = tools::evaluate_polynomial(P, guess);91 int i_exp3 = i_exp / 3;92 93 using shift_type = typename largest_cbrt_int_type<T>::type;94 95 static_assert( ::std::numeric_limits<shift_type>::radix == 2, "The radix of the type to shift to must be 2.");96 97 if(abs(i_exp3) < std::numeric_limits<shift_type>::digits)98 {99 if(i_exp3 > 0)100 guess *= shift_type(1u) << i_exp3;101 else102 guess /= shift_type(1u) << -i_exp3;103 }104 else105 {106 guess = ldexp(guess, i_exp3);107 }108 i_exp %= 3;109 guess *= correction[i_exp + 2];110 //111 // Now inline Halley iteration.112 // We do this here rather than calling tools::halley_iterate since we can113 // simplify the expressions algebraically, and don't need most of the error114 // checking of the boilerplate version as we know in advance that the function115 // is well behaved...116 //117 using prec = typename policies::precision<T, Policy>::type;118 constexpr auto prec3 = prec::value / 3;119 constexpr auto new_prec = prec3 + 3;120 using new_policy = typename policies::normalise<Policy, policies::digits2<new_prec>>::type;121 //122 // Epsilon calculation uses compile time arithmetic when it's available for type T,123 // otherwise uses ldexp to calculate at runtime:124 //125 T eps = (new_prec > 3) ? policies::get_epsilon<T, new_policy>() : ldexp(T(1), -2 - tools::digits<T>() / 3);126 T diff;127 128 if(original_i_exp < std::numeric_limits<T>::max_exponent - 3)129 {130 //131 // Safe from overflow, use the fast method:132 //133 do134 {135 T g3 = guess * guess * guess;136 diff = (g3 + z + z) / (g3 + g3 + z);137 guess *= diff;138 }139 while(fabs(1 - diff) > eps);140 }141 else142 {143 //144 // Either we're ready to overflow, or we can't tell because numeric_limits isn't145 // available for type T:146 //147 do148 {149 T g2 = guess * guess;150 diff = (g2 - z / guess) / (2 * guess + z / g2);151 guess -= diff;152 }153 while((guess * eps) < fabs(diff));154 }155 156 return sign * guess;157}158 159} // namespace detail160 161template <typename T, typename Policy>162BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type cbrt(T z, const Policy& pol)163{164 using result_type = typename tools::promote_args<T>::type;165 using value_type = typename policies::evaluation<result_type, Policy>::type;166 return static_cast<result_type>(detail::cbrt_imp(value_type(z), pol));167}168 169template <typename T>170BOOST_MATH_GPU_ENABLED inline typename tools::promote_args<T>::type cbrt(T z)171{172 return cbrt(z, policies::policy<>());173}174 175} // namespace math176} // namespace boost177 178#else // Special NVRTC handling179 180namespace boost {181namespace math {182 183template <typename T>184BOOST_MATH_GPU_ENABLED double cbrt(T x)185{186 return ::cbrt(x);187}188 189BOOST_MATH_GPU_ENABLED inline float cbrt(float x)190{191 return ::cbrtf(x);192}193 194template <typename T, typename Policy>195BOOST_MATH_GPU_ENABLED double cbrt(T x, const Policy&)196{197 return ::cbrt(x);198}199 200template <typename Policy>201BOOST_MATH_GPU_ENABLED float cbrt(float x, const Policy&)202{203 return ::cbrtf(x);204}205 206} // namespace math207} // namespace boost208 209#endif // NVRTC210 211#endif // BOOST_MATH_SF_CBRT_HPP212 213 214 215 216