Optimize FP64 comparison functions on GPUs with slow FP64 performance

libcudacxx/benchmarks/bench/comparison/basic.cu

@@ -0,0 +1,101 @@
+//===----------------------------------------------------------------------===//
+//
+// Part of libcu++, the C++ Standard Library for your entire system,
+// under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+// SPDX-FileCopyrightText: Copyright (c) 2026 NVIDIA CORPORATION & AFFILIATES.
+//
+//===----------------------------------------------------------------------===//
+
+#include <cuda/std/__floating_point/storage.h>
+#include <cuda/std/cmath>
+
+#include <cuda_runtime_api.h>
+#include <device_side_benchmark.cuh>
+
+#include <nvbench/nvbench.cuh>
+
+struct double_pair
+{
+  double x{};
+  double y{};
+};
+
+struct isless_op
+{
+  __device__ __forceinline__ bool operator()(double x, double y) const
+  {
+    return cuda::std::isless(x, y);
+  }
+};
+
+struct islessequal_op
+{
+  __device__ __forceinline__ bool operator()(double x, double y) const
+  {
+    return cuda::std::islessequal(x, y);
+  }
+};
+
+struct isgreater_op
+{
+  __device__ __forceinline__ bool operator()(double x, double y) const
+  {
+    return cuda::std::isgreater(x, y);
+  }
+};
+
+struct isgreaterequal_op
+{
+  __device__ __forceinline__ bool operator()(double x, double y) const
+  {
+    return cuda::std::isgreaterequal(x, y);
+  }
+};
+
+struct islessgreater_op
+{
+  __device__ __forceinline__ bool operator()(double x, double y) const
+  {
+    return cuda::std::islessgreater(x, y);
+  }
+};
+
+template <typename CompareOp>
+struct benchmark_op_t
+{
+  __device__ __forceinline__ double_pair operator()(double_pair value) const
+  {
+    // both operands change every iteration, so the compiler cannot hoist work on a constant operand
+    // storage is used to prevent code elimination by the compiler (same approach as in isnan)
+    // __fp_get_storage and __fp_from_storage are no op
+    const auto storage = cuda::std::__fp_get_storage(value.x);
+    return {value.y, cuda::std::__fp_from_storage<double>(storage + CompareOp{}(value.x, value.y))};
+  }
+};
+
+template <typename CompareOp>
+void comparison_bench(nvbench::state& state, nvbench::type_list<CompareOp>)
+{
+  constexpr int block_size    = 256;
+  constexpr int unroll_factor = 1024;
+  const auto& kernel    = benchmark_kernel<block_size, unroll_factor, benchmark_op_t<CompareOp>, double_pair, false>;
+  const int num_SMs     = state.get_device().value().get_number_of_sms(); // NOLINT(bugprone-unchecked-optional-access)
+  int max_blocks_per_SM = 0;
+  NVBENCH_CUDA_CALL_NOEXCEPT(cudaOccupancyMaxActiveBlocksPerMultiprocessor(&max_blocks_per_SM, kernel, block_size, 0));
+  const int grid_size = max_blocks_per_SM * num_SMs;
+  state.exec(nvbench::exec_tag::gpu | nvbench::exec_tag::no_batch, [&](nvbench::launch&) {
+    kernel<<<grid_size, block_size>>>(benchmark_op_t<CompareOp>{});
+  });
+}
+
+NVBENCH_DECLARE_TYPE_STRINGS(isless_op, "isless", "");
+NVBENCH_DECLARE_TYPE_STRINGS(islessequal_op, "islessequal", "");
+NVBENCH_DECLARE_TYPE_STRINGS(isgreater_op, "isgreater", "");
+NVBENCH_DECLARE_TYPE_STRINGS(isgreaterequal_op, "isgreaterequal", "");
+NVBENCH_DECLARE_TYPE_STRINGS(islessgreater_op, "islessgreater", "");
+
+using ops = nvbench::type_list<isless_op, islessequal_op, isgreater_op, isgreaterequal_op, islessgreater_op>;
+
+NVBENCH_BENCH_TYPES(comparison_bench, NVBENCH_TYPE_AXES(ops)).set_name("base").set_type_axes_names({"op"});

libcudacxx/benchmarks/bench/fmin_fmax/basic.cu

@@ -0,0 +1,76 @@
+//===----------------------------------------------------------------------===//
+//
+// Part of libcu++, the C++ Standard Library for your entire system,
+// under the Apache License v2.0 with LLVM Exceptions.
+// See https://llvm.org/LICENSE.txt for license information.
+// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
+// SPDX-FileCopyrightText: Copyright (c) 2026 NVIDIA CORPORATION & AFFILIATES.
+//
+//===----------------------------------------------------------------------===//
+
+#include <cuda/std/__floating_point/storage.h>
+#include <cuda/std/cmath>
+
+#include <cuda_runtime_api.h>
+#include <device_side_benchmark.cuh>
+
+#include <nvbench/nvbench.cuh>
+
+struct double_pair
+{
+  double x{};
+  double y{};
+};
+
+struct fmax_op
+{
+  __device__ __forceinline__ double operator()(double x, double y) const
+  {
+    return cuda::std::fmax(x, y);
+  }
+};
+
+struct fmin_op
+{
+  __device__ __forceinline__ double operator()(double x, double y) const
+  {
+    return cuda::std::fmin(x, y);
+  }
+};
+
+template <typename MinMaxOp>
+struct benchmark_op_t
+{
+  __device__ __forceinline__ double_pair operator()(double_pair value) const
+  {
+    // both operands change every iteration, so the compiler cannot hoist work on a constant operand
+    // storage is used to prevent code elimination by the compiler (adapted from isnan benchmark)
+    // the result is mixed with the bits of x, otherwise the chain settles on a single value (x == y)
+    // __fp_get_storage and __fp_from_storage are no op
+    const auto storage = cuda::std::__fp_get_storage(value.x);
+    const auto result  = cuda::std::__fp_get_storage(MinMaxOp{}(value.x, value.y));
+    return {value.y, cuda::std::__fp_from_storage<double>(result ^ (storage >> 7))};
+  }
+};
+
+template <typename MinMaxOp>
+void min_max_bench(nvbench::state& state, nvbench::type_list<MinMaxOp>)
+{
+  constexpr int block_size    = 256;
+  constexpr int unroll_factor = 1024;
+  const auto& kernel    = benchmark_kernel<block_size, unroll_factor, benchmark_op_t<MinMaxOp>, double_pair, false>;
+  const int num_SMs     = state.get_device().value().get_number_of_sms(); // NOLINT(bugprone-unchecked-optional-access)
+  int max_blocks_per_SM = 0;
+  NVBENCH_CUDA_CALL_NOEXCEPT(cudaOccupancyMaxActiveBlocksPerMultiprocessor(&max_blocks_per_SM, kernel, block_size, 0));
+  const int grid_size = max_blocks_per_SM * num_SMs;
+  state.exec(nvbench::exec_tag::gpu | nvbench::exec_tag::no_batch, [&](nvbench::launch&) {
+    kernel<<<grid_size, block_size>>>(benchmark_op_t<MinMaxOp>{});
+  });
+}
+
+NVBENCH_DECLARE_TYPE_STRINGS(fmax_op, "fmax", "");
+NVBENCH_DECLARE_TYPE_STRINGS(fmin_op, "fmin", "");
+
+using ops = nvbench::type_list<fmax_op, fmin_op>;
+
+NVBENCH_BENCH_TYPES(min_max_bench, NVBENCH_TYPE_AXES(ops)).set_name("base").set_type_axes_names({"op"});

libcudacxx/include/cuda/std/__cmath/min_max.h

@@ -24,18 +24,45 @@
 #include <cuda/__type_traits/is_floating_point.h>
 #include <cuda/std/__cmath/isnan.h>
 #include <cuda/std/__concepts/concept_macros.h>
+#include <cuda/std/__floating_point/storage.h>
+#include <cuda/std/__limits/numeric_limits.h>
 #include <cuda/std/__type_traits/conditional.h>
 #include <cuda/std/__type_traits/is_extended_arithmetic.h>
 #include <cuda/std/__type_traits/is_integral.h>
 #include <cuda/std/__type_traits/is_same.h>
 #include <cuda/std/__type_traits/promote.h>
+#include <cuda/std/cstdint>
 
 #include <nv/target>
 
 #include <cuda/std/__cccl/prologue.h>
 
 _CCCL_BEGIN_NAMESPACE_CUDA_STD
 
+// maps the bits of a double to an int64_t in IEEE totalOrder: -0 < +0, NaNs are placed beyond +-inf
+[[nodiscard]] _CCCL_API constexpr int64_t __fp64_to_total_order_int(double __x) noexcept
+{
+  const auto __x_as_int = static_cast<int64_t>(::cuda::std::__fp_get_storage(__x));
+  return __x_as_int < 0 ? __x_as_int ^ ::cuda::std::numeric_limits<int64_t>::max() : __x_as_int;
+}
+
+// fmax (_IsMax) / fmin for fp64 with integer instructions only
+// -0 < +0 is preserved
+// single NaN is ignored; two NaNs return a fixed quiet NaN
+template <bool _IsMax>
+[[nodiscard]] _CCCL_API constexpr double __fp64_fmax_fmin_int(double __x, double __y) noexcept
+{
+  const auto __x_bits     = ::cuda::std::__fp_get_storage(__x);
+  const auto __y_bits     = ::cuda::std::__fp_get_storage(__y);
+  const bool __x_nan      = ::cuda::std::isnan(__x);
+  const bool __y_nan      = ::cuda::std::isnan(__y);
+  const bool __x_lt_y     = ::cuda::std::__fp64_to_total_order_int(__x) < ::cuda::std::__fp64_to_total_order_int(__y);
+  const bool __pick_y     = __x_nan ? !__y_nan : (!__y_nan && __x_lt_y == _IsMax);
+  const auto __result     = __pick_y ? __y_bits : __x_bits;
+  const auto __quiet_bits = ::cuda::std::__fp_get_storage(::cuda::std::numeric_limits<double>::quiet_NaN());
+  return ::cuda::std::__fp_from_storage<double>((__x_nan && __y_nan) ? __quiet_bits : __result);
+}
+
 /***********************************************************************************************************************
  * fmax
  **********************************************************************************************************************/
@@ -47,7 +74,7 @@ _CCCL_REQUIRES(is_floating_point_v<_Tp>)
 [[nodiscard]] _CCCL_API _Tp __with_builtin_fmax(_Tp __x, _Tp __y) noexcept
 {
 #  if _CCCL_COMPILER(CLANG) && _CCCL_HOST_ARCH(ARM64) && _CCCL_HOST_COMPILATION()
-  // Clang lowers fmax to AArch64 fmaxnm, which returns quite NaN for signaling NaNs,
+  // Clang lowers fmax to AArch64 fmaxnm, which returns quiet NaN for signaling NaNs,
   // See https://github.com/llvm/llvm-project/issues/176624
   if (::cuda::std::isnan(__x))
   {
@@ -126,6 +153,10 @@ _CCCL_REQUIRES(__is_extended_arithmetic_v<_Tp>)
           NV_IF_TARGET(NV_PROVIDES_SM_100, (return ::__nv_fp128_fmax(__x, __y);))
         }
 #endif // _CCCL_HAS_FLOAT128()
+        if constexpr (is_same_v<_Tp, double>)
+        {
+          _CCCL_DISPATCH_SLOW_FP64_TARGET(return ::cuda::std::__fp64_fmax_fmin_int<true>(__x, __y);)
+        }
 #if _CCCL_USE_BUILTIN_FMAX()
         if constexpr (is_floating_point_v<_Tp>)
         {
@@ -185,7 +216,7 @@ _CCCL_REQUIRES(is_floating_point_v<_Tp>)
 [[nodiscard]] _CCCL_API _Tp __with_builtin_fmin(_Tp __x, _Tp __y) noexcept
 {
 #  if _CCCL_COMPILER(CLANG) && _CCCL_HOST_ARCH(ARM64) && _CCCL_HOST_COMPILATION()
-  // Clang lowers fmin to AArch64 fminnm, which returns quite NaN for signaling NaNs,
+  // Clang lowers fmin to AArch64 fminnm, which returns quiet NaN for signaling NaNs,
   // See https://github.com/llvm/llvm-project/issues/176624
   if (::cuda::std::isnan(__x))
   {
@@ -264,6 +295,10 @@ _CCCL_REQUIRES(__is_extended_arithmetic_v<_Tp>)
           NV_IF_TARGET(NV_PROVIDES_SM_100, (return ::__nv_fp128_fmin(__x, __y);))
         }
 #endif // _CCCL_HAS_FLOAT128()
+        if constexpr (is_same_v<_Tp, double>)
+        {
+          _CCCL_DISPATCH_SLOW_FP64_TARGET(return ::cuda::std::__fp64_fmax_fmin_int<false>(__x, __y);)
+        }
 #if _CCCL_USE_BUILTIN_FMIN()
         if constexpr (is_floating_point_v<_Tp>)
         {