rsqrt.hpp 1.2 KB

123456789101112131415161718192021222324252627282930313233343536373839404142
  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. #ifndef BOOST_MATH_SPECIAL_FUNCTIONS_RSQRT_HPP
  6. #define BOOST_MATH_SPECIAL_FUNCTIONS_RSQRT_HPP
  7. #include <cmath>
  8. namespace boost::math {
  9. template<typename Real>
  10. inline Real rsqrt(Real const & x)
  11. {
  12. using std::sqrt;
  13. if constexpr (std::is_same_v<Real, float> || std::is_same_v<Real, double> || std::is_same_v<Real, long double>)
  14. {
  15. return 1/sqrt(x);
  16. }
  17. else
  18. {
  19. // if it's so tiny it rounds to 0 as long double,
  20. // no performance gains are possible:
  21. if (x < std::numeric_limits<long double>::denorm_min() || x > std::numeric_limits<long double>::max()) {
  22. return 1/sqrt(x);
  23. }
  24. Real x0 = 1/sqrt(static_cast<long double>(x));
  25. // Divide by 512 for leeway:
  26. Real s = sqrt(std::numeric_limits<Real>::epsilon())*x0/512;
  27. Real x1 = x0 + x0*(1-x*x0*x0)/2;
  28. while(abs(x1 - x0) > s) {
  29. x0 = x1;
  30. x1 = x0 + x0*(1-x*x0*x0)/2;
  31. }
  32. // Final iteration get ~2ULPs:
  33. return x1 + x1*(1-x*x1*x1)/2;;
  34. }
  35. }
  36. }
  37. #endif