diff --git a/libcudacxx/include/cuda/std/__random/fisher_f_distribution.h b/libcudacxx/include/cuda/std/__random/fisher_f_distribution.h new file mode 100644 index 00000000000..b2dd093fca8 --- /dev/null +++ b/libcudacxx/include/cuda/std/__random/fisher_f_distribution.h @@ -0,0 +1,196 @@ +//===----------------------------------------------------------------------===// +// +// Part of the LLVM Project, 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) 2025 NVIDIA CORPORATION & AFFILIATES. +// +//===----------------------------------------------------------------------===// + +#ifndef _CUDA_STD___FISHER_F_DISTRIBUTION_H +#define _CUDA_STD___FISHER_F_DISTRIBUTION_H + +#include + +#if defined(_CCCL_IMPLICIT_SYSTEM_HEADER_GCC) +# pragma GCC system_header +#elif defined(_CCCL_IMPLICIT_SYSTEM_HEADER_CLANG) +# pragma clang system_header +#elif defined(_CCCL_IMPLICIT_SYSTEM_HEADER_MSVC) +# pragma system_header +#endif // no system header + +#include +#include +#include + +#if !_CCCL_COMPILER(NVRTC) +# include +#endif // !_CCCL_COMPILER(NVRTC) + +#include + +_CCCL_BEGIN_NAMESPACE_CUDA_STD + +template +class fisher_f_distribution +{ + static_assert(__libcpp_random_is_valid_realtype<_RealType>, "RealType must be a supported floating-point type"); + +public: + // types + using result_type = _RealType; + + class param_type + { + private: + result_type __m_ = result_type{1}; + result_type __n_ = result_type{1}; + + public: + using distribution_type = fisher_f_distribution; + + _CCCL_API constexpr explicit param_type(result_type __m = result_type{1}, result_type __n = result_type{1}) noexcept + : __m_{__m} + , __n_{__n} + {} + + [[nodiscard]] _CCCL_API constexpr result_type m() const noexcept + { + return __m_; + } + [[nodiscard]] _CCCL_API constexpr result_type n() const noexcept + { + return __n_; + } + + [[nodiscard]] friend _CCCL_API constexpr bool operator==(const param_type& __x, const param_type& __y) noexcept + { + return __x.__m_ == __y.__m_ && __x.__n_ == __y.__n_; + } +#if _CCCL_STD_VER <= 2017 + [[nodiscard]] friend _CCCL_API constexpr bool operator!=(const param_type& __x, const param_type& __y) noexcept + { + return !(__x == __y); + } +#endif // _CCCL_STD_VER <= 2017 + }; + +private: + param_type __p_{}; + +public: + // constructor and reset functions + constexpr fisher_f_distribution() = default; + + _CCCL_API constexpr explicit fisher_f_distribution(result_type __m, result_type __n = result_type{1}) + : __p_{param_type{__m, __n}} + {} + _CCCL_API constexpr explicit fisher_f_distribution(const param_type& __p) + : __p_{__p} + {} + _CCCL_API constexpr void reset() noexcept {} + + // generating functions + template + [[nodiscard]] _CCCL_API result_type operator()(_URng& __g) + { + return (*this)(__g, __p_); + } + template + [[nodiscard]] _CCCL_API result_type operator()(_URng& __g, const param_type& __p) + { + static_assert(__cccl_random_is_valid_urng<_URng>); + gamma_distribution __gdm{__p.m() * result_type{.5}}; + gamma_distribution __gdn{__p.n() * result_type{.5}}; + return __p.n() * __gdm(__g) / (__p.m() * __gdn(__g)); + } + + // property functions + [[nodiscard]] _CCCL_API constexpr result_type m() const noexcept + { + return __p_.m(); + } + [[nodiscard]] _CCCL_API constexpr result_type n() const noexcept + { + return __p_.n(); + } + + [[nodiscard]] _CCCL_API constexpr param_type param() const noexcept + { + return __p_; + } + _CCCL_API _CCCL_CONSTEXPR_CXX20 void param(const param_type& __p) noexcept + { + __p_ = __p; + } + + [[nodiscard]] _CCCL_API static constexpr result_type min() noexcept + { + return result_type{0}; + } + [[nodiscard]] _CCCL_API static constexpr result_type max() noexcept + { + return numeric_limits::infinity(); + } + + [[nodiscard]] friend _CCCL_API constexpr bool + operator==(const fisher_f_distribution& __x, const fisher_f_distribution& __y) noexcept + { + return __x.__p_ == __y.__p_; + } +#if _CCCL_STD_VER <= 2017 + [[nodiscard]] friend _CCCL_API constexpr bool + operator!=(const fisher_f_distribution& __x, const fisher_f_distribution& __y) noexcept + { + return !(__x == __y); + } +#endif // _CCCL_STD_VER <= 2017 + +#if !_CCCL_COMPILER(NVRTC) + template + friend ::std::basic_ostream<_CharT, _Traits>& + operator<<(::std::basic_ostream<_CharT, _Traits>& __os, const fisher_f_distribution& __x) + { + using ostream_type = ::std::basic_ostream<_CharT, _Traits>; + using ios_base = typename ostream_type::ios_base; + const typename ios_base::fmtflags __flags = __os.flags(); + const _CharT __fill = __os.fill(); + const ::std::streamsize __precision = __os.precision(); + __os.flags(ios_base::dec | ios_base::left | ios_base::scientific); + _CharT __sp = __os.widen(' '); + __os.fill(__sp); + __os.precision(numeric_limits::max_digits10); + __os << __x.m() << __sp << __x.n(); + __os.flags(__flags); + __os.fill(__fill); + __os.precision(__precision); + return __os; + } + + template + friend ::std::basic_istream<_CharT, _Traits>& + operator>>(::std::basic_istream<_CharT, _Traits>& __is, fisher_f_distribution& __x) + { + using istream_type = ::std::basic_istream<_CharT, _Traits>; + using ios_base = typename istream_type::ios_base; + const typename ios_base::fmtflags __flags = __is.flags(); + __is.flags(ios_base::dec | ios_base::skipws); + result_type __m; + result_type __n; + __is >> __m >> __n; + if (!__is.fail()) + { + __x.param(param_type{__m, __n}); + } + __is.flags(__flags); + return __is; + } +#endif // !_CCCL_COMPILER(NVRTC) +}; + +_CCCL_END_NAMESPACE_CUDA_STD + +#include + +#endif // _CUDA_STD___FISHER_F_DISTRIBUTION_H diff --git a/libcudacxx/include/cuda/std/__random_ b/libcudacxx/include/cuda/std/__random_ index b252a5fae02..934af395c6a 100644 --- a/libcudacxx/include/cuda/std/__random_ +++ b/libcudacxx/include/cuda/std/__random_ @@ -26,6 +26,7 @@ #include #include #include +#include #include #include #include diff --git a/libcudacxx/test/libcudacxx/std/random/distribution/fisher_f.pass.cpp b/libcudacxx/test/libcudacxx/std/random/distribution/fisher_f.pass.cpp new file mode 100644 index 00000000000..95a1d8c1543 --- /dev/null +++ b/libcudacxx/test/libcudacxx/std/random/distribution/fisher_f.pass.cpp @@ -0,0 +1,61 @@ +//===----------------------------------------------------------------------===// +// +// Part of the LLVM Project, 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) 2025 NVIDIA CORPORATION & AFFILIATES. +// +//===----------------------------------------------------------------------===// +// +// REQUIRES: long_tests + +// + +// template +// class fisher_f_distribution + +#include +#include +#include + +#include "random_utilities/stats_functions.h" +#include "random_utilities/test_distribution.h" +#include "test_macros.h" + +template +struct fisher_f_cdf +{ + using P = typename cuda::std::fisher_f_distribution::param_type; + + __host__ __device__ double operator()(double x, const P& p) const + { + // CDF: F(x; m, n) = I_{mx/(mx+n)}(m/2, n/2) + // where I is the regularized incomplete beta function + if (x <= 0) + { + return 0.0; + } + double m_val = p.m(); + double n_val = p.n(); + double z = (m_val * x) / (m_val * x + n_val); + return incomplete_beta(m_val / 2.0, n_val / 2.0, z); + } +}; + +template +__host__ __device__ void test() +{ + [[maybe_unused]] const bool test_constexpr = false; + using D = cuda::std::fisher_f_distribution; + using P = typename D::param_type; + using G = cuda::std::philox4x64; + cuda::std::array params = {P(1, 1), P(5, 2), P(10, 10), P(2, 5), P(20, 30)}; + test_distribution(params, fisher_f_cdf{}); +} + +int main(int, char**) +{ + test(); + test(); + return 0; +}