seq.hpp

Go to the documentation of this file.

// SPDX-License-Identifier: BSD-2-Clause
 
#ifndef pRC_CORE_RANDOM_SEQ_H
#define pRC_CORE_RANDOM_SEQ_H
 
#include <prc/core/container/array_stack.hpp>
#include <prc/core/random/type_traits.hpp>
 
namespace pRC
{
    template<Size S>
    class SeedSequence
    {
    private:
        static_assert(IsSame<Seed, std::uint32_t>());
 
    public:
        using result_type = Seed;
 
        static constexpr auto size()
        {
            return S;
        }
 
    public:
        template<class... Xs, If<All<IsConvertible<Xs, Seed>...>> = 0>
        constexpr SeedSequence(Xs &&...seeds)
            : mSeed(forward<Xs>(seeds)...)
        {
        }
 
        template<Size N>
        constexpr auto generate() const
        {
            StackArray<Seed, N> out;
            for(Size i = 0; i < N; ++i)
            {
                out[i] = 0x8b8b8b8bU;
            }
 
            constexpr Size t = (N >= 623) ? 11
                : (N >= 68)               ? 7
                : (N >= 39)               ? 5
                : (N >= 7)                ? 3
                                          : (N - 1) / 2;
            constexpr Size p = (N - t) / 2;
            constexpr Size q = p + t;
            constexpr Size m = max(S + 1, N);
 
            auto const T = [&out, &p](auto const k)
            {
                Seed const arg =
                    out[k % N] ^ out[(k + p) % N] ^ out[(k - 1) % N];
                return arg ^ (arg >> 27);
            };
 
            for(Size k = 0; k < m; ++k)
            {
                Seed const r1 = 1664525U * T(k);
                Seed const r2 = r1 +
                    [&]()
                {
                    if(k == 0)
                    {
                        return S;
                    }
                    else if(k <= S)
                    {
                        return k % N + mSeed[k - 1];
                    }
                    else
                    {
                        return k % N;
                    }
                }();
 
                out[(k + p) % N] += r1;
                out[(k + q) % N] += r2;
                out[k % N] = r2;
            }
 
            for(Size k = m; k < m + N; ++k)
            {
                Seed const r3 = 1566083941U * T(k);
                Seed const r4 = r3 - k % N;
 
                out[(k + p) % N] ^= r3;
                out[(k + q) % N] ^= r4;
                out[k % N] = r4;
            }
 
            return out;
        }
 
        constexpr auto &param() const
        {
            return mSeed;
        }
 
    private:
        BEGIN_IGNORE_DIAGNOSTIC_GCC("-Wpedantic")
        StackArray<Seed, S> mSeed;
        END_IGNORE_DIAGNOSTIC_GCC
    };
 
    template<class... Ts>
    SeedSequence(Ts &&...) -> SeedSequence<sizeof...(Ts)>;
}
#endif // pRC_CORE_RANDOM_SEQ_H