From 6b48f06f0710da3b8a49123b2bfbdb60de4e1b50 Mon Sep 17 00:00:00 2001
From: Matt Borland <matt@mattborland.com>
Date: Tue, 3 Jan 2023 17:12:49 +0100
Subject: [PATCH] Improve handling of NANs in fdim

---
 include/boost/math/ccmath/fdim.hpp | 23 +++++++++++++++--------
 1 file changed, 15 insertions(+), 8 deletions(-)

diff --git a/include/boost/math/ccmath/fdim.hpp b/include/boost/math/ccmath/fdim.hpp
index 5b03d05c8..0ef4e6e06 100644
--- a/include/boost/math/ccmath/fdim.hpp
+++ b/include/boost/math/ccmath/fdim.hpp
@@ -17,7 +17,7 @@ namespace boost::math::ccmath {
 namespace detail {
 
 template <typename T>
-inline constexpr T fdim_impl(const T x, const T y) noexcept
+constexpr T fdim_impl(const T x, const T y) noexcept
 {
     if (x <= y)
     {
@@ -36,13 +36,20 @@ inline constexpr T fdim_impl(const T x, const T y) noexcept
 } // Namespace detail
 
 template <typename Real, std::enable_if_t<!std::is_integral_v<Real>, bool> = true>
-inline constexpr Real fdim(Real x, Real y) noexcept
+constexpr Real fdim(Real x, Real y) noexcept
 {
     if (BOOST_MATH_IS_CONSTANT_EVALUATED(x))
     {
-        return boost::math::ccmath::isnan(x) ? std::numeric_limits<Real>::quiet_NaN() :
-               boost::math::ccmath::isnan(y) ? std::numeric_limits<Real>::quiet_NaN() :
-               boost::math::ccmath::detail::fdim_impl(x, y);
+        if (boost::math::ccmath::isnan(x))
+        {
+            return x;
+        }
+        else if (boost::math::ccmath::isnan(y))
+        {
+            return y;
+        }
+
+        return boost::math::ccmath::detail::fdim_impl(x, y);
     }
     else
     {
@@ -52,7 +59,7 @@ inline constexpr Real fdim(Real x, Real y) noexcept
 }
 
 template <typename T1, typename T2>
-inline constexpr auto fdim(T1 x, T2 y) noexcept
+constexpr auto fdim(T1 x, T2 y) noexcept
 {
     if (BOOST_MATH_IS_CONSTANT_EVALUATED(x))
     {
@@ -83,13 +90,13 @@ inline constexpr auto fdim(T1 x, T2 y) noexcept
     }
 }
 
-inline constexpr float fdimf(float x, float y) noexcept
+constexpr float fdimf(float x, float y) noexcept
 {
     return boost::math::ccmath::fdim(x, y);
 }
 
 #ifndef BOOST_MATH_NO_LONG_DOUBLE_MATH_FUNCTIONS
-inline constexpr long double fdiml(long double x, long double y) noexcept
+constexpr long double fdiml(long double x, long double y) noexcept
 {
     return boost::math::ccmath::fdim(x, y);
 }