math/special_functions/rsqrt.hpp

0001 //  (C) Copyright Nick Thompson 2020.
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_SPECIAL_FUNCTIONS_RSQRT_HPP
0007 #define BOOST_MATH_SPECIAL_FUNCTIONS_RSQRT_HPP
0008 #include <cmath>
0009 #include <type_traits>
0010 #include <limits>
0011
0012 #include <boost/math/tools/is_standalone.hpp>
0013 #ifndef BOOST_MATH_STANDALONE
0014 #  include <boost/config.hpp>
0015 #  ifdef BOOST_NO_CXX17_IF_CONSTEXPR
0016 #    error "The header <boost/math/rqrt.hpp> can only be used in C++17 and later."
0017 #  endif
0018 #endif
0019
0020 namespace boost::math {
0021
0022 template<typename Real>
0023 inline Real rsqrt(Real const & x)
0024 {
0025     using std::sqrt;
0026     if constexpr (std::is_arithmetic_v<Real> && !std::is_integral_v<Real>)
0027     {
0028         return 1/sqrt(x);
0029     }
0030     else
0031     {
0032         // if it's so tiny it rounds to 0 as long double,
0033         // no performance gains are possible:
0034         if (x < std::numeric_limits<long double>::denorm_min() || x > (std::numeric_limits<long double>::max)()) {
0035             return 1/sqrt(x);
0036         }
0037         Real x0 = 1/sqrt(static_cast<long double>(x));
0038         // Divide by 512 for leeway:
0039         Real s = sqrt(std::numeric_limits<Real>::epsilon())*x0/512;
0040         Real x1 = x0 + x0*(1-x*x0*x0)/2;
0041         while(abs(x1 - x0) > s) {
0042             x0 = x1;
0043             x1 = x0 + x0*(1-x*x0*x0)/2;
0044         }
0045         // Final iteration get ~2ULPs:
0046         return  x1 + x1*(1-x*x1*x1)/2;;
0047     }
0048 }
0049
0050
0051 }
0052 #endif