Optimize FP64 fabs on GPUs with low FP64 performance

libcudacxx/benchmarks/bench/fabs/basic.cu

@@ -0,0 +1,47 @@
+//===----------------------------------------------------------------------===//
+//
+// 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/bit>
+#include <cuda/std/cmath>
+
+#include <cuda_runtime_api.h>
+#include <device_side_benchmark.cuh>
+
+#include <nvbench/nvbench.cuh>
+
+struct benchmark_op_t
+{
+  __device__ __forceinline__ double operator()(double value) const
+  {
+    // storage is used to prevent code elimination by the compiler (adapted from isnan benchmark)
+    // the rotation moves a mantissa bit into the sign bit, so the next input can be negative again
+    // __fp_get_storage and __fp_from_storage are no op
+    const auto storage = cuda::std::__fp_get_storage(value);
+    const auto result  = cuda::std::__fp_get_storage(cuda::std::fabs(value));
+    return cuda::std::__fp_from_storage<double>(result ^ cuda::std::rotr(storage, 7));
+  }
+};
+
+void fabs_bench(nvbench::state& state)
+{
+  constexpr int block_size    = 256;
+  constexpr int unroll_factor = 1024;
+  const auto& kernel          = benchmark_kernel<block_size, unroll_factor, benchmark_op_t, double, 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{});
+  });
+}
+
+NVBENCH_BENCH(fabs_bench).set_name("base");

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

@@ -25,8 +25,11 @@
 #include <cuda/std/__floating_point/fp.h>
 #include <cuda/std/__type_traits/is_extended_arithmetic.h>
 #include <cuda/std/__type_traits/is_integral.h>
+#include <cuda/std/cstdint>
 #include <cuda/std/limits>
 
+#include <nv/target>
+
 #include <cuda/std/__cccl/prologue.h>
 
 #if _CCCL_CHECK_BUILTIN(builtin_fabs) || _CCCL_COMPILER(GCC)
@@ -56,6 +59,20 @@ _CCCL_BEGIN_NAMESPACE_CUDA_STD
 
 // fabs
 
+// fabs for fp64 with integer instructions only: clears the sign bit as IEEE 754 abs, NaNs keep their payload and are
+// not quieted
+[[nodiscard]] _CCCL_DEVICE_API inline double __fp64_fabs_int(double __x) noexcept
+{
+  // On SM100+, the compiler recognizes clearing the sign bit of a 64-bit integer as fabs, which reintroduces fp64
+  // instructions. Clearing it in the high 32-bit half avoids that
+  NV_IF_ELSE_TARGET(
+    NV_PROVIDES_SM_100,
+    (return ::__hiloint2double(::__double2hiint(__x) & ::cuda::std::numeric_limits<int32_t>::max(),
+                               ::__double2loint(__x));),
+    (return ::cuda::std::__fp_from_storage<double>(
+              ::cuda::std::__fp_get_storage(__x) & ::cuda::std::__fp_exp_mant_mask_of_v<double>);))
+}
+
 _CCCL_TEMPLATE(class _Tp)
 _CCCL_REQUIRES(__is_extended_arithmetic_v<_Tp>)
 [[nodiscard]] _CCCL_API constexpr auto fabs(_Tp __x) noexcept
@@ -77,6 +94,14 @@ _CCCL_REQUIRES(__is_extended_arithmetic_v<_Tp>)
   }
   else
   {
+    if constexpr (is_same_v<_Tp, double>)
+    {
+      // fabs for fp64 with integer instructions on slow FP64 targets
+      _CCCL_IF_NOT_CONSTEVAL_DEFAULT
+      {
+        _CCCL_DISPATCH_SLOW_FP64_TARGET(return ::cuda::std::__fp64_fabs_int(__x);)
+      }
+    }
 #ifdef _CCCL_BUILTIN_FABSF
     if constexpr (is_same_v<_Tp, float>)
     {

libcudacxx/test/libcudacxx/std/numerics/c.math/fp_abs.pass.cpp

@@ -14,8 +14,10 @@
 #include <disable_nvfp_conversions_and_operators.h>
 // clang-format on
 
+#include <cuda/std/bit>
 #include <cuda/std/cassert>
 #include <cuda/std/cmath>
+#include <cuda/std/cstdint>
 #include <cuda/std/limits>
 #include <cuda/std/type_traits>
 #include <cuda/type_traits>
@@ -210,10 +212,39 @@ TEST_FUNC constexpr bool test(float val)
   return true;
 }
 
+// double has an integer implementation on slow fp64 targets
+// the NaN sign and payload after fabs are implementation-defined on device, only check that the result is a NaN
+TEST_FUNC void test_fp64_nan(float val)
+{
+  // the inputs are built at runtime: fabs of a constant sNaN is evaluated at compile time, which quiets it
+  const auto runtime               = static_cast<cuda::std::uint64_t>(val);
+  const cuda::std::uint64_t nans[] = {
+    0x7ff8'0000'0000'0000 * runtime, // +qNaN
+    0xfff8'0000'0000'0000 * runtime, // -qNaN
+    0x7ff4'0000'0000'0000 * runtime, // +sNaN
+    0xfff4'0000'0000'0000 * runtime, // -sNaN
+    0xfff0'0000'0000'0001 * runtime, // -sNaN, smallest payload
+    0xfff8'1234'5678'0001 * runtime, // -qNaN with payload
+  };
+
+  for (const auto nan_bits : nans)
+  {
+    assert(cuda::std::isnan(cuda::std::fabs(cuda::std::bit_cast<double>(nan_bits))));
+  }
+}
+
+// -0.0 == 0.0, so the sign bit of fabs(-0.0) has to be checked explicitly
+TEST_FUNC void test_fp64_neg_zero(float val)
+{
+  assert(!cuda::std::signbit(cuda::std::fabs(-0.0 * val)));
+}
+
 int main(int, char**)
 {
   volatile float val = 1.0f;
   test(val);
+  test_fp64_nan(val);
+  test_fp64_neg_zero(val);
 #if _CCCL_HAS_CONSTEXPR_BIT_CAST()
   static_assert(test(1.0f));
 #endif // _CCCL_HAS_CONSTEXPR_BIT_CAST()