123456789101112131415161718192021222324252627282930313233343536373839404142 |
- // (C) Copyright Nick Thompson 2020.
- // 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_SPECIAL_FUNCTIONS_RSQRT_HPP
- #define BOOST_MATH_SPECIAL_FUNCTIONS_RSQRT_HPP
- #include <cmath>
- namespace boost::math {
- template<typename Real>
- inline Real rsqrt(Real const & x)
- {
- using std::sqrt;
- if constexpr (std::is_same_v<Real, float> || std::is_same_v<Real, double> || std::is_same_v<Real, long double>)
- {
- return 1/sqrt(x);
- }
- else
- {
- // if it's so tiny it rounds to 0 as long double,
- // no performance gains are possible:
- if (x < std::numeric_limits<long double>::denorm_min() || x > std::numeric_limits<long double>::max()) {
- return 1/sqrt(x);
- }
- Real x0 = 1/sqrt(static_cast<long double>(x));
- // Divide by 512 for leeway:
- Real s = sqrt(std::numeric_limits<Real>::epsilon())*x0/512;
- Real x1 = x0 + x0*(1-x*x0*x0)/2;
- while(abs(x1 - x0) > s) {
- x0 = x1;
- x1 = x0 + x0*(1-x*x0*x0)/2;
- }
- // Final iteration get ~2ULPs:
- return x1 + x1*(1-x*x1*x1)/2;;
- }
- }
- }
- #endif
|