canonical.hpp

Go to the documentation of this file.

// SPDX-License-Identifier: BSD-2-Clause
 
#ifndef pRC_CORE_RANDOM_CANONICAL_H
#define pRC_CORE_RANDOM_CANONICAL_H
 
#include <prc/core/basic/common.hpp>
#include <prc/core/basic/functions/ceil_div.hpp>
#include <prc/core/basic/functions/ilog.hpp>
#include <prc/core/basic/limits.hpp>
#include <prc/core/random/type_traits.hpp>
#include <prc/core/value/type_traits.hpp>
 
namespace pRC
{
    template<class T, class RNG, If<IsRandomEngine<RNG>> = 0,
        If<IsFloat<T>> = 0>
    static inline constexpr T generateCanonical(RNG &rng)
    {
        constexpr auto b = NumericLimits<T>::digits();
 
        static_assert(RNG::max() - RNG::min() < NumericLimits<Size>::max(),
            "Unsigned integer overflow.");
 
        constexpr Size R = RNG::max() - RNG::min() + Size(1);
 
        constexpr Size logR = iLog<2>(R);
        constexpr Size k = max(ceilDiv(b, logR), Size(1));
 
        T base = R;
        T basic = rng() - RNG::min();
 
        for(Index i = 1; i < k; ++i, base *= T(R))
        {
            basic += T(rng() - RNG::min()) * base;
        }
        basic /= base;
 
        if(basic >= T(1))
        {
            return T(T(1) - NumericLimits<T>::epsilon() / T(2));
        }
        return basic;
    }
}
#endif // pRC_CORE_RANDOM_CANONICAL_H