distribution/detail/GenerateCanonical32.hh

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