diff --git a/include/xtd/math.h b/include/xtd/math.h index 468b5f9..23ea81d 100644 --- a/include/xtd/math.h +++ b/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" diff --git a/include/xtd/math/rsqrt.h b/include/xtd/math/rsqrt.h new file mode 100644 index 0000000..00472ef --- /dev/null +++ b/include/xtd/math/rsqrt.h @@ -0,0 +1,68 @@ +/* + * Copyright 2026 European Organization for Nuclear Research (CERN) + * Authors: Andrea Bocci , Aurora Perego , Simone Balducci + * SPDX-License-Identifier: MPL-2.0 + */ + +#pragma once + +#include +#include + +#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(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(arg)); + } + XTD_DEVICE_FUNCTION inline constexpr float rsqrtf(std::integral auto arg) { + return xtd::rsqrt(static_cast(arg)); + } + +} // namespace xtd diff --git a/test/src/math/rsqrt/mpfr_rsqrt.h b/test/src/math/rsqrt/mpfr_rsqrt.h new file mode 100644 index 0000000..04f4bfc --- /dev/null +++ b/test/src/math/rsqrt/mpfr_rsqrt.h @@ -0,0 +1,26 @@ +/* + * Copyright 2026 European Organization for Nuclear Research (CERN) + * Authors: Andrea Bocci + * SPDX-License-Identifier: MPL-2.0 + */ + +// mpfr::real headers +#include + +// xtd headers +#include + +// test headers +#include "common/mpfr.h" + +inline float mpfr_rsqrtf(xtd::arithmetic auto arg) { + float result; + (1. / mpfr::sqrt(static_cast(static_cast(arg)))).conv(result); + return result; +} + +inline double mpfr_rsqrt(xtd::arithmetic auto arg) { + double result; + (1. / mpfr::sqrt(static_cast(static_cast(arg)))).conv(result); + return result; +} diff --git a/test/src/math/rsqrt/rsqrt_t.cc b/test/src/math/rsqrt/rsqrt_t.cc new file mode 100644 index 0000000..4e1028d --- /dev/null +++ b/test/src/math/rsqrt/rsqrt_t.cc @@ -0,0 +1,52 @@ +/* + * Copyright 2026 European Organization for Nuclear Research (CERN) + * Authors: Andrea Bocci , Aurora Perego + * SPDX-License-Identifier: MPL-2.0 + */ + +// C++ standard headers +#include +#include + +// Catch2 headers +#include + +// 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(device, ulps_single); + } + + SECTION("double xtd::rsqrt(double)") { + validate(device, ulps_double); + } + + SECTION("double xtd::rsqrt(int)") { + validate(device, ulps_double); + } + + SECTION("float xtd::rsqrtf(float)") { + validate(device, ulps_single); + } + + SECTION("float xtd::rsqrtf(double)") { + validate(device, ulps_single); + } + + SECTION("float xtd::rsqrtf(int)") { + validate(device, ulps_single); + } + } +} diff --git a/test/src/math/rsqrt/rsqrt_t.cu b/test/src/math/rsqrt/rsqrt_t.cu new file mode 100644 index 0000000..fda5c56 --- /dev/null +++ b/test/src/math/rsqrt/rsqrt_t.cu @@ -0,0 +1,53 @@ +/* + * Copyright 2026 European Organization for Nuclear Research (CERN) + * Authors: Andrea Bocci , Aurora Perego + * SPDX-License-Identifier: MPL-2.0 + */ + +// Catch2 headers +#define CATCH_CONFIG_NO_POSIX_SIGNALS +#include + +// 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(device, ulps_single); + } + + SECTION("double xtd::rsqrt(double)") { + validate(device, ulps_double); + } + + SECTION("double xtd::rsqrt(int)") { + validate(device, ulps_double); + } + + SECTION("float xtd::rsqrtf(float)") { + validate(device, ulps_single); + } + + SECTION("float xtd::rsqrtf(double)") { + validate(device, ulps_single); + } + + SECTION("float xtd::rsqrtf(int)") { + validate(device, ulps_single); + } + } + } + } +} diff --git a/test/src/math/rsqrt/rsqrt_t.hip.cc b/test/src/math/rsqrt/rsqrt_t.hip.cc new file mode 100644 index 0000000..c71f9ad --- /dev/null +++ b/test/src/math/rsqrt/rsqrt_t.hip.cc @@ -0,0 +1,53 @@ +/* + * Copyright 2026 European Organization for Nuclear Research (CERN) + * Authors: Andrea Bocci , Aurora Perego + * SPDX-License-Identifier: MPL-2.0 + */ + +// Catch2 headers +#define CATCH_CONFIG_NO_POSIX_SIGNALS +#include + +// 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(device, ulps_single); + } + + SECTION("double xtd::rsqrt(double)") { + validate(device, ulps_double); + } + + SECTION("double xtd::rsqrt(int)") { + validate(device, ulps_double); + } + + SECTION("float xtd::rsqrtf(float)") { + validate(device, ulps_single); + } + + SECTION("float xtd::rsqrtf(double)") { + validate(device, ulps_single); + } + + SECTION("float xtd::rsqrtf(int)") { + validate(device, ulps_single); + } + } + } + } +} diff --git a/test/src/math/rsqrt/rsqrt_t.sycl.cc b/test/src/math/rsqrt/rsqrt_t.sycl.cc new file mode 100644 index 0000000..08ddb3c --- /dev/null +++ b/test/src/math/rsqrt/rsqrt_t.sycl.cc @@ -0,0 +1,55 @@ +/* + * Copyright 2026 European Organization for Nuclear Research (CERN) + * Authors: Andrea Bocci , Aurora Perego + * SPDX-License-Identifier: MPL-2.0 + */ + +// Catch2 headers +#define CATCH_CONFIG_NO_POSIX_SIGNALS +#include + +// 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(platform, device, ulps_single); + } + + SECTION("double xtd::rsqrt(double)") { + validate(platform, device, ulps_double); + } + + SECTION("double xtd::rsqrt(int)") { + validate(platform, device, ulps_double); + } + + SECTION("float xtd::rsqrtf(float)") { + validate(platform, device, ulps_single); + } + + SECTION("float xtd::rsqrtf(double)") { + validate(platform, device, ulps_single); + } + + SECTION("float xtd::rsqrtf(int)") { + validate(platform, device, ulps_single); + } + } + } + } + } +}