Implement xtd::rsqrt()

include/xtd/math.h

@@ -32,6 +32,7 @@
 #include "math/cbrt.h"
 #include "math/hypot.h"
 //hypot(x,y,z)
+#include "math/rsqrt.h"
 
 // Trigonometric functions
 #include "math/sin.h"

include/xtd/math/rsqrt.h

@@ -0,0 +1,68 @@
+/*
+ * Copyright 2026 European Organization for Nuclear Research (CERN)
+ * Authors: Andrea Bocci <andrea.bocci@cern.ch>, Aurora Perego <aurora.perego@cern.ch>, Simone Balducci <simone.balducci@cern.ch>
+ * SPDX-License-Identifier: MPL-2.0
+ */
+
+#pragma once
+
+#include <concepts>
+#include <cmath>
+
+#include "xtd/internal/defines.h"
+
+namespace xtd {
+
+  /* Computes the reciprocal of the square root of arg, in single precision.
+   */
+  XTD_DEVICE_FUNCTION inline constexpr float rsqrt(float arg) {
+#if defined(XTD_TARGET_CUDA)
+    // CUDA device code
+    // Note: __frsqrt_rn() is correctly rounded, while rsqrtf() is rounded to 2 ULPs
+    return ::__frsqrt_rn(arg);
+#elif defined(XTD_TARGET_HIP)
+    // HIP/ROCm device code
+    return ::rsqrtf(arg);
+#elif defined(XTD_TARGET_SYCL)
+    // SYCL device code
+    return sycl::rsqrt(arg);
+#else
+    // standard C/C++ code
+    return 1.f / ::sqrtf(arg);
+#endif
+  }
+
+  /* Computes the reciprocal of the square root of arg, in double precision.
+   */
+  XTD_DEVICE_FUNCTION inline constexpr double rsqrt(double arg) {
+#if defined(XTD_TARGET_CUDA)
+    // CUDA device code
+    return ::rsqrt(arg);
+#elif defined(XTD_TARGET_HIP)
+    // HIP/ROCm device code
+    return ::rsqrt(arg);
+#elif defined(XTD_TARGET_SYCL)
+    // SYCL device code
+    return sycl::rsqrt(arg);
+#else
+    // standard C/C++ code
+    return 1. / ::sqrt(arg);
+#endif
+  }
+
+  /* Computes the reciprocal of the square root of arg, in double precision.
+   */
+  XTD_DEVICE_FUNCTION inline constexpr double rsqrt(std::integral auto arg) {
+    return xtd::rsqrt(static_cast<double>(arg));
+  }
+
+  /* Computes the reciprocal of the square root of arg, in single precision.
+   */
+  XTD_DEVICE_FUNCTION inline constexpr float rsqrtf(std::floating_point auto arg) {
+    return xtd::rsqrt(static_cast<float>(arg));
+  }
+  XTD_DEVICE_FUNCTION inline constexpr float rsqrtf(std::integral auto arg) {
+    return xtd::rsqrt(static_cast<float>(arg));
+  }
+
+}  // namespace xtd

test/src/math/rsqrt/mpfr_rsqrt.h

@@ -0,0 +1,26 @@
+/*
+ * Copyright 2026 European Organization for Nuclear Research (CERN)
+ * Authors: Andrea Bocci <andrea.bocci@cern.ch>
+ * SPDX-License-Identifier: MPL-2.0
+ */
+
+// mpfr::real headers
+#include <real.hpp>
+
+// xtd headers
+#include <xtd/concepts/arithmetic.h>
+
+// test headers
+#include "common/mpfr.h"
+
+inline float mpfr_rsqrtf(xtd::arithmetic auto arg) {
+  float result;
+  (1. / mpfr::sqrt(static_cast<mpfr_single>(static_cast<float>(arg)))).conv(result);
+  return result;
+}
+
+inline double mpfr_rsqrt(xtd::arithmetic auto arg) {
+  double result;
+  (1. / mpfr::sqrt(static_cast<mpfr_double>(static_cast<double>(arg)))).conv(result);
+  return result;
+}

test/src/math/rsqrt/rsqrt_t.cc

@@ -0,0 +1,52 @@
+/*
+ * Copyright 2026 European Organization for Nuclear Research (CERN)
+ * Authors: Andrea Bocci <andrea.bocci@cern.ch>, Aurora Perego <aurora.perego@cern.ch>
+ * SPDX-License-Identifier: MPL-2.0
+ */
+
+// C++ standard headers
+#include <cmath>
+#include <vector>
+
+// Catch2 headers
+#include <catch.hpp>
+
+// xtd headers
+#include "xtd/math/rsqrt.h"
+
+// test headers
+#include "common/cpu/device.h"
+#include "common/cpu/validate.h"
+#include "mpfr_rsqrt.h"
+
+constexpr int ulps_single = 0;
+constexpr int ulps_double = 0;
+
+TEST_CASE("xtd::rsqrt", "[rsqrt][cpu]") {
+  const auto& device = test::cpu::device();
+  DYNAMIC_SECTION("CPU: " << device.name()) {
+    SECTION("float xtd::rsqrt(float)") {
+      validate<float, float, xtd::rsqrt, mpfr_rsqrtf>(device, ulps_single);
+    }
+
+    SECTION("double xtd::rsqrt(double)") {
+      validate<double, double, xtd::rsqrt, mpfr_rsqrt>(device, ulps_double);
+    }
+
+    SECTION("double xtd::rsqrt(int)") {
+      validate<double, int, xtd::rsqrt, mpfr_rsqrt>(device, ulps_double);
+    }
+
+    SECTION("float xtd::rsqrtf(float)") {
+      validate<float, float, xtd::rsqrtf, mpfr_rsqrtf>(device, ulps_single);
+    }
+
+    SECTION("float xtd::rsqrtf(double)") {
+      validate<float, double, xtd::rsqrtf, mpfr_rsqrtf>(device, ulps_single);
+    }
+
+    SECTION("float xtd::rsqrtf(int)") {
+      validate<float, int, xtd::rsqrtf, mpfr_rsqrtf>(device, ulps_single);
+    }
+  }
+}

test/src/math/rsqrt/rsqrt_t.cu

@@ -0,0 +1,53 @@
+/*
+ * Copyright 2026 European Organization for Nuclear Research (CERN)
+ * Authors: Andrea Bocci <andrea.bocci@cern.ch>, Aurora Perego <aurora.perego@cern.ch>
+ * SPDX-License-Identifier: MPL-2.0
+ */
+
+// Catch2 headers
+#define CATCH_CONFIG_NO_POSIX_SIGNALS
+#include <catch.hpp>
+
+// xtd headers
+#include "xtd/math/rsqrt.h"
+
+// test headers
+#include "common/cuda/platform.h"
+#include "common/cuda/validate.h"
+#include "mpfr_rsqrt.h"
+
+constexpr int ulps_single = 0;
+constexpr int ulps_double = 1;
+
+TEST_CASE("xtd::rsqrt", "[rsqrt][cuda]") {
+  const auto& platform = test::cuda::platform();
+  DYNAMIC_SECTION("CUDA platform: " << platform.name()) {
+    for (const auto& device : platform.devices()) {
+      DYNAMIC_SECTION("CUDA device " << device.index() << ": " << device.name()) {
+        SECTION("float xtd::rsqrt(float)") {
+          validate<float, float, xtd::rsqrt, mpfr_rsqrtf>(device, ulps_single);
+        }
+
+        SECTION("double xtd::rsqrt(double)") {
+          validate<double, double, xtd::rsqrt, mpfr_rsqrt>(device, ulps_double);
+        }
+
+        SECTION("double xtd::rsqrt(int)") {
+          validate<double, int, xtd::rsqrt, mpfr_rsqrt>(device, ulps_double);
+        }
+
+        SECTION("float xtd::rsqrtf(float)") {
+          validate<float, float, xtd::rsqrtf, mpfr_rsqrtf>(device, ulps_single);
+        }
+
+        SECTION("float xtd::rsqrtf(double)") {
+          validate<float, double, xtd::rsqrtf, mpfr_rsqrtf>(device, ulps_single);
+        }
+
+        SECTION("float xtd::rsqrtf(int)") {
+          validate<float, int, xtd::rsqrtf, mpfr_rsqrtf>(device, ulps_single);
+        }
+      }
+    }
+  }
+}

test/src/math/rsqrt/rsqrt_t.hip.cc

@@ -0,0 +1,53 @@
+/*
+ * Copyright 2026 European Organization for Nuclear Research (CERN)
+ * Authors: Andrea Bocci <andrea.bocci@cern.ch>, Aurora Perego <aurora.perego@cern.ch>
+ * SPDX-License-Identifier: MPL-2.0
+ */
+
+// Catch2 headers
+#define CATCH_CONFIG_NO_POSIX_SIGNALS
+#include <catch.hpp>
+
+// xtd headers
+#include "xtd/math/rsqrt.h"
+
+// test headers
+#include "common/hip/platform.h"
+#include "common/hip/validate.h"
+#include "mpfr_rsqrt.h"
+
+constexpr int ulps_single = 1;
+constexpr int ulps_double = 1;
+
+TEST_CASE("xtd::rsqrt", "[rsqrt][hip]") {
+  const auto& platform = test::hip::platform();
+  DYNAMIC_SECTION("HIP platform: " << platform.name()) {
+    for (const auto& device : platform.devices()) {
+      DYNAMIC_SECTION("HIP device " << device.index() << ": " << device.name()) {
+        SECTION("float xtd::rsqrt(float)") {
+          validate<float, float, xtd::rsqrt, mpfr_rsqrtf>(device, ulps_single);
+        }
+
+        SECTION("double xtd::rsqrt(double)") {
+          validate<double, double, xtd::rsqrt, mpfr_rsqrt>(device, ulps_double);
+        }
+
+        SECTION("double xtd::rsqrt(int)") {
+          validate<double, int, xtd::rsqrt, mpfr_rsqrt>(device, ulps_double);
+        }
+
+        SECTION("float xtd::rsqrtf(float)") {
+          validate<float, float, xtd::rsqrtf, mpfr_rsqrtf>(device, ulps_single);
+        }
+
+        SECTION("float xtd::rsqrtf(double)") {
+          validate<float, double, xtd::rsqrtf, mpfr_rsqrtf>(device, ulps_single);
+        }
+
+        SECTION("float xtd::rsqrtf(int)") {
+          validate<float, int, xtd::rsqrtf, mpfr_rsqrtf>(device, ulps_single);
+        }
+      }
+    }
+  }
+}

test/src/math/rsqrt/rsqrt_t.sycl.cc

@@ -0,0 +1,55 @@
+/*
+ * Copyright 2026 European Organization for Nuclear Research (CERN)
+ * Authors: Andrea Bocci <andrea.bocci@cern.ch>, Aurora Perego <aurora.perego@cern.ch>
+ * SPDX-License-Identifier: MPL-2.0
+ */
+
+// Catch2 headers
+#define CATCH_CONFIG_NO_POSIX_SIGNALS
+#include <catch.hpp>
+
+// xtd headers
+#include "xtd/math/rsqrt.h"
+
+// test headers
+#include "common/sycl/device.h"
+#include "common/sycl/platform.h"
+#include "common/sycl/validate.h"
+#include "mpfr_rsqrt.h"
+
+constexpr int ulps_single = 1;
+constexpr int ulps_double = 1;
+
+TEST_CASE("xtd::rsqrt", "[rsqrt][sycl]") {
+  for (const auto &platform : test::sycl::platforms()) {
+    DYNAMIC_SECTION("SYCL platform " << platform.index() << ": " << platform.name()) {
+      for (const auto &device : platform.devices()) {
+        DYNAMIC_SECTION("SYCL device " << platform.index() << '.' << device.index() << ": " << device.name()) {
+          SECTION("float xtd::rsqrt(float)") {
+            validate<float, float, xtd::rsqrt, mpfr_rsqrtf>(platform, device, ulps_single);
+          }
+
+          SECTION("double xtd::rsqrt(double)") {
+            validate<double, double, xtd::rsqrt, mpfr_rsqrt>(platform, device, ulps_double);
+          }
+
+          SECTION("double xtd::rsqrt(int)") {
+            validate<double, int, xtd::rsqrt, mpfr_rsqrt>(platform, device, ulps_double);
+          }
+
+          SECTION("float xtd::rsqrtf(float)") {
+            validate<float, float, xtd::rsqrtf, mpfr_rsqrtf>(platform, device, ulps_single);
+          }
+
+          SECTION("float xtd::rsqrtf(double)") {
+            validate<float, double, xtd::rsqrtf, mpfr_rsqrtf>(platform, device, ulps_single);
+          }
+
+          SECTION("float xtd::rsqrtf(int)") {
+            validate<float, int, xtd::rsqrtf, mpfr_rsqrtf>(platform, device, ulps_single);
+          }
+        }
+      }
+    }
+  }
+}