1 // (C) Copyright Nick Thompson 2020. 2 // Use, modification and distribution are subject to the 3 // Boost Software License, Version 1.0. (See accompanying file 4 // LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt) 5 6 #ifndef BOOST_MATH_SPECIAL_FUNCTIONS_RSQRT_HPP 7 #define BOOST_MATH_SPECIAL_FUNCTIONS_RSQRT_HPP 8 #include <cmath> 9 10 namespace boost::math { 11 12 template<typename Real> rsqrt(Real const & x)13inline Real rsqrt(Real const & x) 14 { 15 using std::sqrt; 16 if constexpr (std::is_same_v<Real, float> || std::is_same_v<Real, double> || std::is_same_v<Real, long double>) 17 { 18 return 1/sqrt(x); 19 } 20 else 21 { 22 // if it's so tiny it rounds to 0 as long double, 23 // no performance gains are possible: 24 if (x < std::numeric_limits<long double>::denorm_min() || x > std::numeric_limits<long double>::max()) { 25 return 1/sqrt(x); 26 } 27 Real x0 = 1/sqrt(static_cast<long double>(x)); 28 // Divide by 512 for leeway: 29 Real s = sqrt(std::numeric_limits<Real>::epsilon())*x0/512; 30 Real x1 = x0 + x0*(1-x*x0*x0)/2; 31 while(abs(x1 - x0) > s) { 32 x0 = x1; 33 x1 = x0 + x0*(1-x*x0*x0)/2; 34 } 35 // Final iteration get ~2ULPs: 36 return x1 + x1*(1-x*x1*x1)/2;; 37 } 38 } 39 40 41 } 42 #endif 43