Update remquo to better handle specific edge cases

include/ccmath/detail/basic/impl/remquo_double_impl.hpp

@@ -28,26 +28,32 @@ namespace ccm::internal
 		{
 			inline constexpr double remquo_double_impl(double x, double y, int * quo) noexcept
 			{
-				// Assume all edge case checking is done before calling this function.
 				std::int64_t x_i64{};
 				std::int64_t y_i64{};
 				std::uint64_t x_sign{};
 				std::uint64_t quotient_sign{};
 				int computed_quotient{};
 
-				x_i64		  = ccm::helpers::double_to_int64(x);
-				y_i64		  = ccm::helpers::double_to_int64(y);
+				x_i64 = ccm::helpers::double_to_int64(x);
+				y_i64 = ccm::helpers::double_to_int64(y);
 
 				// Determine the signs of x and the quotient.
 				x_sign		  = static_cast<std::uint64_t>(x_i64) & 0x8000000000000000ULL;
 				quotient_sign = x_sign ^ (static_cast<std::uint64_t>(y_i64) & 0x8000000000000000ULL);
 
 				// Clear the sign bits from the int64_t representations of x and y.
-				y_i64 &= 0x7fffffffffffffffULL;
 				x_i64 &= 0x7fffffffffffffffULL;
+				y_i64 &= 0x7fffffffffffffffULL;
+
+				// If y is zero.
+				if (CCM_UNLIKELY(y_i64 == 0)) { return (x * y) / (x * y); }
+
+				// If x is not finite or y is NaN.
+				if (CCM_UNLIKELY(x_i64 >= 0x7ff0000000000000ULL || y_i64 > 0x7ff0000000000000ULL)) { return (x * y) / (x * y); }
 
 				// b (or bit 54) represents the highest bit we can compare against for both signed and unsigned integers without causing an overflow.
-				// Here we are checking that if y_i64 is within the range of signed 64-bit integers that can be represented without setting the MSB (i.e., positive or zero).
+				// Here we are checking that if y_i64 is within the range of signed 64-bit integers that can be represented without setting the MSB (i.e.,
+				// positive or zero).
 				if (y_i64 <= 0x7fbfffffffffffffULL)
 				{
 					x = ccm::fmod(x, 8 * y); // now x < (8 * y)
@@ -63,7 +69,6 @@ namespace ccm::internal
 				y				  = ccm::helpers::int64_to_double(y_i64);
 				computed_quotient = 0;
 
-
 				if (y_i64 <= 0x7fcfffffffffffffULL && x >= 4 * y)
 				{
 					x -= 4 * y;
@@ -116,7 +121,7 @@ namespace ccm::internal
 	}	  // namespace
 
 	template <typename T, typename = std::enable_if_t<std::is_floating_point_v<T>>>
-	inline constexpr T remquo_double(T x, T y, int* quo) noexcept
+	inline constexpr T remquo_double(T x, T y, int * quo) noexcept
 	{
 		return static_cast<T>(impl::remquo_double_impl(x, y, quo));
 	}

include/ccmath/detail/basic/impl/remquo_float_impl.hpp

@@ -28,7 +28,6 @@ namespace ccm::internal
 		{
 			inline constexpr float remquo_float_impl(float x, float y, int * quo) noexcept
 			{
-				// Assume all edge case checking is done before calling this function.
 				std::int32_t x_i32{};
 				std::int32_t y_i32{};
 				std::uint32_t x_sign{};
@@ -43,8 +42,14 @@ namespace ccm::internal
 				quotient_sign = x_sign ^ (static_cast<std::uint32_t>(y_i32) & 0x80000000);
 
 				// Clear the sign bits from the int32_t representations of x and y.
-				y_i32 &= 0x7fffffff;
 				x_i32 &= 0x7fffffff;
+				y_i32 &= 0x7fffffff;
+
+				// If y is zero.
+				if (CCM_UNLIKELY(y_i32 == 0)) { return (x * y) / (x * y); }
+
+				// If x is not finite or y is NaN.
+				if (CCM_UNLIKELY(x_i32 >= 0x7f800000 || y_i32 > 0x7f800000)) { return (x * y) / (x * y); }
 
 				if (y_i32 <= 0x7dffffff)
 				{

include/ccmath/detail/basic/remquo.hpp

@@ -44,28 +44,6 @@ namespace ccm
 	template <typename T, std::enable_if_t<std::is_floating_point_v<T>, bool> = true>
 	inline constexpr T remquo(T x, T y, int * quo)
 	{
-		// We are not using __builtin_remquo with GCC due to a failure for it to work with static_assert
-		// Our implementation does not have this issue.
-
-		// If x is ±∞ and y is not NaN, NaN is returned.
-		if (CCM_UNLIKELY(ccm::isinf(x) && !ccm::isnan(y))) { return (x * y) / (x * y); }
-
-		// If y is ±0 and x is not NaN, NaN is returned.
-		if (CCM_UNLIKELY(y == 0.0 && !ccm::isnan(x))) { return (x * y) / (x * y); }
-
-		// If either x or y is NaN, NaN is returned.
-		if (CCM_UNLIKELY(ccm::isnan(x)))
-		{
-			if (ccm::signbit(x)) { return -std::numeric_limits<T>::quiet_NaN(); }
-			else { return std::numeric_limits<T>::quiet_NaN(); }
-		}
-
-		if (CCM_UNLIKELY(ccm::isnan(y)))
-		{
-			if (ccm::signbit(y)) { return -std::numeric_limits<T>::quiet_NaN(); }
-			else { return std::numeric_limits<T>::quiet_NaN(); }
-		}
-
 		if constexpr (std::is_same_v<T, float>) { return ccm::internal::remquo_float(x, y, quo); }
 		else { return ccm::internal::remquo_double(x, y, quo); }
 	}