math/ccmath/sqrt.hpp

0001 //  (C) Copyright Matt Borland 2021.
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 //  Constexpr implementation of sqrt function
0007
0008 #ifndef BOOST_MATH_CCMATH_SQRT
0009 #define BOOST_MATH_CCMATH_SQRT
0010
0011 #include <boost/math/ccmath/detail/config.hpp>
0012
0013 #ifdef BOOST_MATH_NO_CCMATH
0014 #error "The header <boost/math/sqrt.hpp> can only be used in C++17 and later."
0015 #endif
0016
0017 #include <boost/math/ccmath/abs.hpp>
0018 #include <boost/math/ccmath/isnan.hpp>
0019 #include <boost/math/ccmath/isinf.hpp>
0020 #include <boost/math/tools/is_constant_evaluated.hpp>
0021
0022 namespace boost::math::ccmath {
0023
0024 namespace detail {
0025
0026 template <typename Real>
0027 constexpr Real sqrt_impl_2(Real x, Real s, Real s2)
0028 {
0029     return !(s < s2) ? s2 : sqrt_impl_2(x, (x / s + s) / 2, s);
0030 }
0031
0032 template <typename Real>
0033 constexpr Real sqrt_impl_1(Real x, Real s)
0034 {
0035     return sqrt_impl_2(x, (x / s + s) / 2, s);
0036 }
0037
0038 template <typename Real>
0039 constexpr Real sqrt_impl(Real x)
0040 {
0041     return sqrt_impl_1(x, x > 1 ? x : Real(1));
0042 }
0043
0044 } // namespace detail
0045
0046 template <typename Real, std::enable_if_t<!std::is_integral_v<Real>, bool> = true>
0047 constexpr Real sqrt(Real x)
0048 {
0049     if(BOOST_MATH_IS_CONSTANT_EVALUATED(x))
0050     {
0051         if (boost::math::ccmath::isnan(x) ||
0052            (boost::math::ccmath::isinf(x) && x > 0) ||
0053             boost::math::ccmath::abs(x) == Real(0))
0054         {
0055             return x;
0056         }
0057         // Domain error is implementation defined so return NAN
0058         else if (boost::math::ccmath::isinf(x) && x < 0)
0059         {
0060             return std::numeric_limits<Real>::quiet_NaN();
0061         }
0062
0063         return detail::sqrt_impl<Real>(x);
0064     }
0065     else
0066     {
0067         using std::sqrt;
0068         return sqrt(x);
0069     }
0070 }
0071
0072 template <typename Z, std::enable_if_t<std::is_integral_v<Z>, bool> = true>
0073 constexpr double sqrt(Z x)
0074 {
0075     return detail::sqrt_impl<double>(static_cast<double>(x));
0076 }
0077
0078 } // Namespaces
0079
0080 #endif // BOOST_MATH_CCMATH_SQRT