53 lines · plain
1// (C) Copyright Nick Thompson 2020.2// Use, modification and distribution are subject to the3// Boost Software License, Version 1.0. (See accompanying file4// LICENSE_1_0.txt or copy at http://www.boost.org/LICENSE_1_0.txt)5 6#ifndef BOOST_MATH_SPECIAL_FUNCTIONS_RSQRT_HPP7#define BOOST_MATH_SPECIAL_FUNCTIONS_RSQRT_HPP8#include <cmath>9#include <type_traits>10#include <limits>11 12#include <boost/math/tools/is_standalone.hpp>13#ifndef BOOST_MATH_STANDALONE14# include <boost/config.hpp>15# ifdef BOOST_MATH_NO_CXX17_IF_CONSTEXPR16# error "The header <boost/math/rqrt.hpp> can only be used in C++17 and later."17# endif18#endif19 20namespace boost::math {21 22template<typename Real>23inline Real rsqrt(Real const & x)24{25 using std::sqrt;26 if constexpr (std::is_arithmetic_v<Real> && !std::is_integral_v<Real>)27 {28 return 1/sqrt(x);29 }30 else31 {32 // if it's so tiny it rounds to 0 as long double,33 // no performance gains are possible:34 if (x < std::numeric_limits<long double>::denorm_min() || x > (std::numeric_limits<long double>::max)()) {35 return 1/sqrt(x);36 }37 Real x0 = 1/sqrt(static_cast<long double>(x));38 // Divide by 512 for leeway:39 Real s = sqrt(std::numeric_limits<Real>::epsilon())*x0/512;40 Real x1 = x0 + x0*(1-x*x0*x0)/2;41 while(abs(x1 - x0) > s) {42 x0 = x1;43 x1 = x0 + x0*(1-x*x0*x0)/2;44 }45 // Final iteration get ~2ULPs:46 return x1 + x1*(1-x*x1*x1)/2;;47 }48}49 50 51}52#endif53