random/detail/GenerateCanonical32.hh

0001 //----------------------------------*-C++-*----------------------------------//
0002 // Copyright 2022-2024 UT-Battelle, LLC, and other Celeritas developers.
0003 // See the top-level COPYRIGHT file for details.
0004 // SPDX-License-Identifier: (Apache-2.0 OR MIT)
0005 //---------------------------------------------------------------------------//
0006 //! \file celeritas/random/detail/GenerateCanonical32.hh
0007 //---------------------------------------------------------------------------//
0008 #pragma once
0009
0010 #include "corecel/Types.hh"
0011
0012 namespace celeritas
0013 {
0014 namespace detail
0015 {
0016 //---------------------------------------------------------------------------//
0017 /*!
0018  * Generate random numbers in [0, 1) from a 32-bit generator.
0019  */
0020 template<class RealType = ::celeritas::real_type>
0021 class GenerateCanonical32;
0022
0023 template<>
0024 class GenerateCanonical32<float>
0025 {
0026   public:
0027     //!@{
0028     //! \name Type aliases
0029     using result_type = float;
0030     //!@}
0031
0032   public:
0033     //! Sample a random number with floating point precision
0034     template<class Generator>
0035     inline CELER_FUNCTION result_type operator()(Generator& rng);
0036 };
0037
0038 template<>
0039 class GenerateCanonical32<double>
0040 {
0041   public:
0042     //!@{
0043     //! \name Type aliases
0044     using result_type = double;
0045     //!@}
0046
0047   public:
0048     //! Sample a random number with floating point precision
0049     template<class Generator>
0050     inline CELER_FUNCTION result_type operator()(Generator& rng);
0051 };
0052
0053 //---------------------------------------------------------------------------//
0054 /*!
0055  * Generate a 32-bit float from one 32-bit sample.
0056  *
0057  * Returns a float on [0, 1) assuming a generator range of [0, 2^32).
0058  */
0059 template<class Generator>
0060 CELER_FUNCTION float GenerateCanonical32<float>::operator()(Generator& rng)
0061 {
0062     static_assert(Generator::max() == 0xffffffffu,
0063                   "Generator must return 32-bit sample");
0064
0065     constexpr float norm = 2.32830643654e-10f;  // 1 / 2**32
0066     return norm * rng();
0067 }
0068
0069 //---------------------------------------------------------------------------//
0070 /*!
0071  * Generate a 64-bit double from two 32-bit samples.
0072  *
0073  * Returns a double on [0, 1).
0074  */
0075 template<class Generator>
0076 CELER_FUNCTION double GenerateCanonical32<double>::operator()(Generator& rng)
0077 {
0078     static_assert(Generator::max() == 0xffffffffu,
0079                   "Generator must return 32-bit sample");
0080     static_assert(sizeof(ull_int) == 8, "Expected 64-bit UL");
0081
0082     unsigned int upper = rng();
0083     unsigned int lower = rng();
0084
0085     // Convert the two 32-byte samples to a 53-bit-precision double by shifting
0086     // the 'upper' and combining the lower bits of the upper with the higher
0087     // bits of the lower
0088     constexpr double norm = 1.1102230246251565e-16;  // 1 / 2^53
0089     return norm
0090            * static_cast<double>((static_cast<ull_int>(upper) << (53ul - 32ul))
0091                                  ^ static_cast<ull_int>(lower));
0092 }
0093
0094 //---------------------------------------------------------------------------//
0095 }  // namespace detail
0096 }  // namespace celeritas