//  (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