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)13 inline 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